Please see below ...
Joe.
-----Original Message----- From: Emil Jerabek [SMTP:ejer5183@artax.karlin.mff.cuni.cz] Sent: Wednesday, February 06, 2002 4:39 AM To: gpc@gnu.de Subject: Re: Complex_Arctan
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;
[Joe da Silva]
Perhaps it's impossible to avoid overflow under all conditions? In this case, if Im(z) is large and Re(z) is small, then "d" will overflow here because "c" is now large too. :-/
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 :-)
[Joe da Silva]
No, it was invented in Australia! See the following! : <G> http://www.theage.com.au/news/state/2001/07/02/FFX0ADFPLOC.html
(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
[Joe da Silva]
Perhaps I should consult the manual first, but what's this "-pedantic" option for? "SqRt", "Odd" and "WriteLn" are all standard identifiers in ISO 7185, so why are they listed as "UCSD extensions"?
Emil Jerabek
da Silva, Joe wrote:
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;
[Joe da Silva]
Perhaps it's impossible to avoid overflow under all conditions? In this case, if Im(z) is large and Re(z) is small, then "d" will overflow here because "c" is now large too. :-/
But that's the `x >= b' case. :-)
It seems that even the wheel was invented in San Diego :-)
[Joe da Silva]
No, it was invented in Australia! See the following! : <G> http://www.theage.com.au/news/state/2001/07/02/FFX0ADFPLOC.html
Nice. I'm surprised, though, that this didn't happen in the US first ... ;-)
Frank