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(a)math.cas.cz; gpc(a)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 ------