Long-float variant giving more precision to numerical computations.

Hello Richard,

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:

(slatec:dbesj0 0.001w0)
0.999999750000015624999565972229008
Expected:
0.999999750000015624999565972229003906182...

(slatec:dbesj0 1w-12)
0.9999999999999999999999997499999904923978
Expected:
0.999999999999999999999999750000000000000000...

> 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.

Regards
Douglas Crosher