Please see below ...
Joe.
-----Original Message----- From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] Sent: Friday, January 25, 2002 5:23 AM To: gpc@gnu.de Subject: Complex_Arctan
Hi all,
I've had a look at the complex arctangent function in the RTS. It turned out that
- it raises a runtime error or returns an infinite imaginary part on
arguments too close to the singularities +/- i, for example ArcTan (Cmplx (1e6 * MinReal, +/- 1))
- it often returns a real number when it should return something with a
small imaginary part instead (example: ArcTan (Cmplx (1e17, 1e17)) on IA32, YMMV), and in general, the imaginary part of the result is very unreliable when |Im (z)| << |Re (z)|^2 (because it is essentially computed as `Ln (something close to 1)'; example: Im (Arctan (Cmplx (1000, 0.0001))) gives 9.999989726e-11 instead of 9.999990000e-11, again this is machine-dependent)
The following patch to p/rts/math.pas solves these problems (I hope it doesn't introduce others).
function Log1P (x: Real): Real; asmname 'log1p';
function Complex_Arctan (z: Complex): Complex; const Ln2 = 0.6931471805599453094172321214581765680755; var b, c, d: Real; begin 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); b := 1 - Abs(Im(z)); d := (1 + Abs(Im(z))) * 0.5 * b - c;
Also, as a general question, is division more expensive than multiplication, or does the complier's optimization take care of such matters anyway?
if Abs (b) > Abs (Re (z)) then b := Log1P (4 * (Abs (Im (z)) / b) / (b + Re (z) * (Re (z) / b))) / 4 else if c <= MinReal then b := (Ln2 - Ln (Abs (Re (z)))) / 2 else b := Log1P (4 * (Abs (Im (z)) / Re (z)) / (b * (b / Re (z)) + Re (z))) / 4; if Im (z) < 0 then b := -b; Complex_Arctan := Cmplx (ArcTan2 (Re (z), d) / 2, b) end;
Emil Jerabek
da Silva, Joe wrote:
Also, as a general question, is division more expensive than multiplication, or does the complier's optimization take care of such matters anyway?
I think for numbers whose inverse can be represented exactly (such as 2 and 0.5), it does. In a quick check I just made, it would, e.g., convert `a := b / 2' to `a := b * 0.5' and `a := b / 0.5' to `a := b + b' with optimization.
Frank
Joe da Silva wrote:
function Log1P (x: Real): Real; asmname 'log1p';
function Complex_Arctan (z: Complex): Complex; const Ln2 = 0.6931471805599453094172321214581765680755; var b, c, d: Real; begin 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); b := 1 - Abs(Im(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. 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?
Also, as a general question, is division more expensive than multiplication, or does the complier's optimization take care of such matters anyway?
AFAIK, the compiler turns Foo / 2 into Foo * 0.5, when optimizing. (This obviously can't take care of division with non-constant divisor, and in such a case division _is_ usually much more expensive than multiplication, although this is somewhat target-dependent.)
Emil
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
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; ------------------------------------
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.
[pas]% cat exprconst.pas program Foo (Output);
const Bar = 1 / Sqrt (16);
var OK: Boolean value True; S: String (30);
procedure Baz (X: Integer); var Y: Real; E: Integer; begin if Odd (X) then WriteStr (S, Bar) else WriteStr (S, Bar); Val (S, Y, E); if (E <> 0) or (Y <> 0.25) then begin WriteLn ('Failed ', X, ': ', S); OK := False end end;
begin Baz (0); Baz (1); Baz (2); Baz (3); if OK then WriteLn ('OK') end. [pas]% gpc exprconst.pas [pas]% ./a.out Failed 0: 4.609351182172047e-01 [pas]% gpc -O2 exprconst.pas [pas]% ./a.out Failed 0: Inf Failed 2: Inf
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 :-) (Apart from that, shouldn't it produce warnings rather than errors?)
[pas]% gpc -pedantic exprconst.pas exprconst.pas:4: `SqRt' is a UCSD Pascal extension exprconst.pas:7: type/var initialization is an ISO 10206 Extended Pascal extension exprconst.pas:8: UCSD and Borland Pascal want string capacity in brackets exprconst.pas:8: ISO 7185 Pascal does not have the `String' type exprconst.pas: In procedure `Baz': exprconst.pas:15: `Odd' is a UCSD Pascal extension exprconst.pas:15: `WriteStr' is an ISO 10206 Extended Pascal extension exprconst.pas:15: `WriteStr' is an ISO 10206 Extended Pascal extension exprconst.pas:16: `WriteStr' is an ISO 10206 Extended Pascal extension exprconst.pas:16: `WriteStr' is an ISO 10206 Extended Pascal extension exprconst.pas:17: `Val' is a Borland Pascal extension exprconst.pas:17: `Val' is a Borland Pascal extension exprconst.pas:19: `WriteLn' is a UCSD Pascal extension exprconst.pas: In main program: exprconst.pas:29: `WriteLn' is a UCSD Pascal extension
Emil Jerabek
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
Frank Heckenbach wrote:
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
So is mine.
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)?
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.
(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 ;-).
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 :-)
Emil Jerabek
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
I wrote:
Emil Jerabek wrote:
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 was able to fix this particular bug. However, such constants still don't work as a bound of a global array. But at least GPC gives a proper error message in this case (`size of ... isn't constant'). But it works with a local array (emil11c.pas).
Frank
-
da Silva, Joe -
Emil Jerabek -
Frank Heckenbach