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 ------
Joe wrote:
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
??? You don't have to guess the coefficients, they are initialized in the code. (If you don't want to rebuild the RTS just to test it (I didn't either), move the "to begin do" stuff into your main program, and import the GPC unit (because of the SplitReal function).)
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(?)
Not sure what sort of "tweaking" you have in mind. Anyway, my Log1P implementation tries to use as few terms of the power series as possible, and it sums them in reverse order. I don't know if this is optimal, but it appeared to work with acceptable precision: the largest error I found was 2 ulp's when compiled without optimization, and 1 ulp when compiled with -O2, both on an i586.
Emil
da Silva, Joe wrote:
Emil Jerabek wrote:
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.
Indeed. So for all checks containing A (1e17) I'm adding similar tests with B=MaxReal (except those whose correct results contain B/2 or Sqr(B) in which cases I'm using MaxReal/2 and SqRt(MaxReal), resp.). Alright? In fact, both the last and the previous version of Complex_ArcTan passes all these tests on IA32 (tested on Linux and DJGPP).
One platform where I'm getting overflows with the previous version (even without the new tests) is alpha-dec-osf4.0b. I'll try the last version there.
Frank
-
da Silva, Joe -
Emil Jerabek -
Frank Heckenbach