From: Raymond Toy <toy.raymond@gm...>  20080413 14:03:14

Juan Jose GarciaRipoll wrote: > The following does not imply anything. I will look into this, but I > presume one needs special functions to handle the high precision case, > as in CMUCL. I also wonder how many digits of accuracy this > systematically brings in other cases. > > The problem is in the calculation of log(sqrt(x^2+y^2)). If x is close to 1 and y is small or vice versa, sqrt(x^2+y^2) is inaccurate. The solution is to change the computation to log(sqrt(x^2+y^2)) = 1/2*log(x^2+y^2) = 1/2*log(1+((x1)*(x+1)+y^2) But don't use the log function to evaluate this. Use log1p(z) = log(1+z). Most C libraries have a log1p function already, but if not, you can just use the series formula since z is small. Ray 