[Maxima] Long-float variant giving more precision to numerical computations.
dtc-maxima at scieneer.com
Tue Feb 26 20:24:13 CST 2008
Richard Fateman wrote:
> If you want to test the accuracy of a Bessel function subroutine, you
> should test it near some difficult points, not at 1.0. Probably a
> difficult place is near a zero of the Bessel function, especially
> a zero that is far from the origin.
Yes, good point. I expect Ray was seeing gross errors due to initialization
problems in the library and I was not attempting to verify the accuracy of
the slatec library. However here are some examples close to zero:
> I would not conclude that Slatec is almost accurate enough for double-double;
> in fact, if the library is really well designed, it should be just barely
> good enough for the maximum precision which was anticipated for its use.
> any higher accuracy is wasted effort!
The library does adjust the number of coefficients used based on the float accuracy
and the number of coefficients is determined at runtime. However the number of
coefficients available is limited. The Fortran source does state the claimed
accuracy and it is close to double-double float precision, for example: dbesj0.f
C Series for BJ0 on the interval 0. to 1.60000E+01
C with weighted error 4.39E-32
C log weighted error 31.36
C significant figures required 31.21
C decimal places required 32.00
Obviously the Slatec library is not going to be adequate for the quad-double
or the precise CLISP long-float representations, but it does appear to offer
better accuracy than double-float and is a quick path to more accurate
results. It would the good to have big-float versions someday.
Here are some timing results for:
(time (dotimes (i 1000000) (slatec:dbesj0 0.001l0)))
double-float (64 bit, CMUCL): 0.49 seconds
long-float (80 bit, SCL): 0.45 seconds
double-double-float (CMUCL): 13.9 seconds
long-float (128 bit, CLISP): 18.18 seconds
The SCL long-float (80 bit) code gives a little more precision the double-float
without a significant performance hit and it is well within the accuracy of the
Slatec library. If more precision is needed than the CMUCL double-double-float
or the CLISP long-float code is an option but is much slower and the accuracy
of the Slatec library may be reached.
More information about the Maxima