Please see below ...
Joe.
-----Original Message----- From: Frank Heckenbach [SMTP:frank@g-n-u.de] Sent: Wednesday, February 06, 2002 11:40 AM To: jerabek@math.cas.cz; gpc@gnu.de Subject: Re: Complex_Arctan
Emil Jerabek wrote:
Frank Heckenbach wrote:
Emil Jerabek wrote:
[...]
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
So, here's another version :-) It should never overflow, AFAICS.
const Ln2 = 0.6931471805599453094172321214581765680755; var FunnyNumber: Real;
function Complex_ArcTan (z: Complex): Complex; var x, y, a, b, c, d: Real; begin x := Abs (Re (z)); y := Abs (Im (z)); b := 1 - y; d := (1 + y) / 2; if b < 0 then d := -d; b := Abs (b); if x >= b then begin c := b / x; d := d * c - x / 2; if x <= FunnyNumber then b := Ln2 - Ln (4 * x * Sqrt (1 + Sqr (c)) / (1 + y)) / 2 else b := Log1P (2 * (y / x) / (b * c / 2 + x / 2)) / 4; c := 1; end else begin c := x / b; a := x * c / 2; d := d - a; if b <= FunnyNumber then b := Ln2 - Ln (4 * b * Sqrt (1 + Sqr (c)) / (1 + y)) / 2 else b := Log1P (2 * (y / b) / (b / 2 + a)) / 4 end; d := ArcTan2 (c, d) / 2; if Re (z) < 0 then d := -d; if Im (z) < 0 then b := -b; Complex_ArcTan := Cmplx (d, b) end;
to begin do begin FunnyNumber := 2.1 / Sqrt (MaxReal); end;
Well, this one is a little longer than the previous ones, so should I use it?
[Joe da Silva]
Probably. It looks more reliable than the previous version.
The only concern I still have, is that the test program, emil9, does not test for accuracy with large "arctan" arguments (only up to 1e17). Since "tan" can produce huge results even with just real arguments, an "arctan" function must conversely provide accurate results with huge arguments.
BTW, AFAICT, DJGPP does not have a Log1P function, so I wasn't able to verify the above anyway. I didn't know the values to use for the coefficients "Log1PCoefs[]", so I couldn't try the suggested Log1P code, either. Although Emil presented a power series which I assume represents these coefficients, I don't know whether these need to be "tweaked". If I understand this correctly, it is an infinite series. However, AFAIK, when you apply a finite version of such a series, you need to "tweak" it's values to get good accuracy (no, I don't know how to do this).
I remember experimenting years ago with the power series for "sine" or "cosine" and being amazed at how a small finite "tweaked" power series was much more accurate than the "conventional" power series, no matter how many terms I calculated for it (after a while, the more terms you calculate, the greater the error you introduce). I presume the "log1p" power series also has this effect(?)
------ snip ------