Emil Jerabek wrote:
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))))) / 2This code is essentially untested (i.e. it works on i586/Linux, but that's not a target with EpsReal^2 <= 5/MaxReal).
As Knuth said, I've only proved the correctness of the code, but never tested it. ;-)
Thanks for the new code.
Joe da Silva wrote:
function Complex_Arctan (z: Complex): Complex; c := Sqr (Re (z)); b := 1 - Abs (Im (z)); d := ((1 + Abs (Im (z))) * b - c) / 2;
[Joe da Silva]
I wonder, would it be better to use the following instead, to reduce problems with overflow for large values (magnitude) of Re(z) and/or Im(z)? :
c := Re(z) * 0.5 * Re(z); d := (1 + Abs(Im(z))) * 0.5 * b - c;It seems you found the second place in the code which makes assumptions about the FP representation (I wasn't aware of this one): it is safe for d to overflow, _provided_ EpsReal^2 > 16 / MaxReal.
... and provided the FP representation supports Inf (not all do AFAIK).
Your change doesn't solve this, it just reduces the threshold to EpsReal^2 > 4 / MaxReal. (BTW, the original implementation in math.pas also needs a similar condition.)
Frank, do you _really_ think there are machines with such a small MaxReal?
It doesn't really matter what I think, does it? ;-) Actually, I haven't looked at any non-IEEE FP implementation in detail.
I'm applying Joe's patch (since it's certainly a little improvement). If you have more ideas for improvement, let me know.
Frank