Frank wrote:
Mathematically it all seems alright. I didn't try to prove the numeric stability, but my tests show that it's better than that of the previous version. Just one thing:
The "if c <= MinReal" branch (needed to avoid some overflows) is executed when |y|=1 and x^2 <= MinReal. In this case,
From the code, the condition for this case is:
| 1 - |y| | <= |x| and x^2 <= MinReal
To conclude |y| = 1, you need that EpsReal^2 >= MinReal, don't you? While this might be true on every FP system, would it cause problems in your code if not?
Frank
The "c <= MinReal" test is completely wrong :-( It should be "c <= 5 / MaxReal".
But back to your question. The code really assumes EpsReal^2 > 5/MaxReal. If you feel that this may be a problem, something more complicated is needed:
if c + Sqr (b) <= 5 / MaxReal then b := -Ln (Abs (Cmplx (Re (z), b)) / (1 + Abs (Im (z)))) / 2 {*} else if Abs (b) > Abs (Re (z)) then b := Log1P (4 * (Abs (Im (z)) / b) / (b + Re (z) * (Re (z) / b))) / 4 else b := Log1P (4 * (Abs (Im (z)) / Re (z)) / (b * (b / Re (z)) + Re (z))) / 4;
This still assumes 1/EpsReal < 1.6*MaxReal. I can't imagine a floating-point implementation violating this; anyway, in such a case, the line with {*} above should be replaced by
b := -Ln (Abs (Cmplx (Re (z), b)) / Abs (Cmplx (Re (z), 1 + Abs (Im (z))))) / 2
This code is essentially untested (i.e. it works on i586/Linux, but that's not a target with EpsReal^2 <= 5/MaxReal).
Emil Jerabek