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. One Quick and Dirty trick is adapted from the Numerical Recipes library (version 2), where the function is called pythag:
function Complex_Abs (z: Complex): Double; var absx, absy: double; begin absx:=abs(Re(z)); absy:=abs(Im(z)); if absx>absy then Complex_Abs:=absx*sqrt(1.0+sqr(absy/absx)) else if absy=0 then Complex_Abs:=0.0 else Complex_Abs:=absy*sqrt(1.0+sqr(absx/absy)); end;
But then I've got to reverse Emil's optimization in Complex_Ln, so it won't fail for the same problem. Or there a better way to do complex Ln?
Now I had a look for this problem, both in the djgpp libc and libm libraries, where there is a function hyphot which does the job, and is said to be ieee754 standard, and in the CEPHES long double math library, (which defines cabs), which I use (with an import unit) when I need accuracy, since djgpp has no long double library (only elementary + - / * operations I think). libc overcomes the problem by using the extended accuracy and exponent range of the x87 floating point unit. gpc cannot use directly the same trick, since sqrt is computed in double accuracy only because djgpp does not define sqrtl (this is tested by configure). One way out would be to use an other configure test for hypot, to use directly this libc function (if this is possible).
Yes, I just hope it's implemented in a "good" way, otherwise it would be worse to use this function compared to the code above.
I'm doing both, and it works on my system. (maur8.pas)
Frank