Maurice Lombardi wrote:
Frank Heckenbach wrote:
Emil Jerabek wrote: If you find other things to improve in the maths algorithms, just go ahead (numerics has never been my favourite area ;-).
While reading this thread, I had a look to math.pas and found an other well known issue. The function
Complex_Abs := SqRt (Sqr (Re (z)) + Sqr (Im (z)))
overflows if either Re or Im is in between sqrt(MaxReal) (~10^154) and MaxReal (~10^308), while the result is (usually) below MaxReal.
Why not apply a little trigonometry.
Theta = arctan(Re/Im); C_Abs = Re * cos(theta);
with a few conditionals to avoid zeroes. Actually, the presence of zeroes greatly simplifies the computation.
You might reserve it for when either Re or Im is greater than sqrt(maxreal)/2 if you have qualms about the trigonometric function accuracy. Note the /2, so the sum doesn't overflow either.
(* only overflows if the result is an overflow value *) IF re = 0.0 THEN cabs := abs(im) ELSE IF im = 0.0 THEN cabs := abs(re) ELSE IF (abs(re) < crit) AND (abs(im < crit) THEN BEGIN (* sqrt of sum of squares, can't overflow *) END ELSE BEGIN (* trig version, can overflow, probably less accurate *) END;