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))))) / 2

This 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

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?

When I pointed out these special cases in previous mails, I didn't mean to says we should cater for all kinds of practically impossible situations, I rather meant to ask you if you see any real problems (since my numeric experience is rather limited). E.g., when using the previous version on a machine that really has such strange properties, would it merely produce inaccurate results in some borderline cases (which is probably acceptable on such a strange machine), or produce clearly visible errors (so a user on such a platform would recognize it and report it as a bug, and we'd know about the problem and where it happens), or could it silently produce substantially wrongs results (which would be a real problem)?

I encountered a compiler bug during testing :-( Strictly speaking, the Standard doesn't allow function calls in constant expressions, but I expect the compiler would reject it if it were not supposed to work.

At least that's a known bug (at least known to us; I'll have to check the to-do list if it's listed there). I'll look into it ... (emil11.pas)

I also wondered if the compiler catches this as a violation of the standard when compiled with `--extended-pascal' or `-pedantic'. It compiled happily with the former (when I replaced `Val' with `ReadStr', of course), but the result with `-pedantic' was somewhat surprising. It seems that even the wheel was invented in San Diego :-)

Of course, the errors about `SqRt', `Ord' and `WriteLn' are wrong (happens for all standard functions because of a systematic bug). Fixing it. (emil10.pas)

(Apart from that, shouldn't it produce warnings rather than errors?)

It did until a recent discussion here when people noted that there should be errors with `--extended-pascal' (or other dialect options which don't contain the given feature), so I changed it. One might want to make it warnings with `-pedantic' and errors with `--extended-pascal', but doing this will be some extra work. If you like to do it, go ahead. ;-) I won't do it for such an unimportant thing (IMHO -- since a program that compiles without problems with `-pedantic' must be in the intersection of ISO 7185, UCSD and Borland Pascal and maybe more, which gives a rather restricted language ;-).

Frank

Emil Jerabek wrote:

That depends. I assume it would just produce more or less inaccurate results, if the machine permits to recover gracefully from overflows, i.e. either it supports infinities, or it returns +/- MaxReal in its stead. Otherwise anything may happen.

Could you construct a test case which would fail (so it would be noticed when running the test suite), or are the checks in emil9.pas (as I recently sent to the list) good enough already?

I don't care much about what `-pedantic' does. It just seemed to me that its behavior is inconsistent with its description in the manual :-)

OK, I'll do the easy thing and adjust the manual. ;-) Thanks for pointing this out.

Frank