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;
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).
libm and CEPHES use basically the same trick to avoid the need for extended exponent range: they extract the exponents, rescale the two arguments (by shifting their exponents) to be in acceptable range, compute and then rescale the result the opposite way. There are also auxiliary issues to ensure lowest bit accuracy described in %DJGPP%\src\libm\math\e_hypot.c (extracted from v2/djlsr203.zip). The libm hypot function could be also retrieved through configure test I think.
Basically what to do depends on what kind of complexity you are willing to incorporate into the compiler itself, as opposed to specialized libraries.
Maurice