Hello list,
the last time I mentioned this bug (or misbehavior maybe?) I didn't get any
feedback. Anyway, I hope I get some feedback now on this patch, which
unfortunately makes complex number division more complex (;-) but fixes
mistaken results in case of denominator overflow. Any ideas for optimizing
this algorithm are welcome.
You may reproduce with:
writeLn(cstr((-8E+231 - 1E+230*i) / (-1E+144 + 4E+199*i)));
without the patch it returns:
-0.000000000000000E+000
but with the patch:
-2.500000000000000E+030+2.000000000000000E+032i
Dimitris
--- rtl/inc/ucomplex.pp.orig3 2006-08-04 04:31:33.000000000 +0300
+++ rtl/inc/ucomplex.pp 2006-11-07 21:28:56.000000000 +0200
@@ -315,14 +315,34 @@
inline;
{$endif TEST_INLINE}
{ division : z := znum / zden }
+ { The following algorithm is used to properly handle
+ denominator overflow:
+
+ | a + b(d/c) c - a(d/c)
+ | ---------- + ---------- I if |d| < |c|
+ a + b I | c + d(d/c) a + d(d/c)
+ ------- = |
+ c + d I | b + a(c/d) -a+ b(c/d)
+ | ---------- + ---------- I if |d| >= |c|
+ | d + c(c/d) d + c(c/d)
+ }
var
- denom : real;
+ tmp, denom : real;
begin
- { WARNING: denom may overflow here, spoiling the result! }
- with zden do denom := (re * re) + (im * im);
- { generates a fpu exception if denom=0 as for reals }
- z.re := ((znum.re * zden.re) + (znum.im * zden.im)) / denom;
- z.im := ((znum.im * zden.re) - (znum.re * zden.im)) / denom;
+ if ( abs(zden.re) > abs(zden.im) ) then
+ begin
+ tmp := zden.im / zden.re;
+ denom := zden.re + zden.im * tmp;
+ z.re := (znum.re + znum.im * tmp) / denom;
+ z.im := (znum.im - znum.re * tmp) / denom;
+ end
+ else
+ begin
+ tmp := zden.re / zden.im;
+ denom := zden.im + zden.re * tmp;
+ z.re := (znum.im + znum.re * tmp) / denom;
+ z.im := (-znum.re + znum.im * tmp) / denom;
+ end;
end;
operator / (znum : complex; r : real) z : complex;
_______________________________________________
fpc-devel maillist - [email protected]
http://lists.freepascal.org/mailman/listinfo/fpc-devel
