GSL K0/K1

[Patrick Alken <alken@colorado.edu>]

Pavel,

  The Bessel function results look good. If everything is tested you can
merge it all into the master branch. I'm cc'ing gsl-discuss so that this
discussion will start to be archived.

Regarding the error estimates, an attempt was made early on in GSL to
provide error estimates for all the special functions. From section 7.2
of the manual:

----
The field val contains the value and the field err contains an estimate
of the absolute error in the value.
----

I confess I don't fully understand all the error estimation, or even how
accurate the estimates are. I don't know if any studies were done to
verify the accuracy of these error estimates. One of the original GSL
developers (Gerard Jungman) wrote a lot of this code (I've cc'd him on
this email) - maybe he can help shed light on this.

Since you've already done exhaustive comparisons of your new Bessel
functions against arbitrary precision libraries, maybe you could simply
change the error estimate to be:

err = factor * GSL_DBL_EPSILON

where factor are the numbers you've found in your error analysis. This
way the calculation would be very fast (and probably more accurate than
the current error estimation).

We can wait and see if Gerard has any suggestions on this too.

Gerard - since you haven't been following this email chain, Pavel has
found more accurate Chebyshev expansions for the Bessel K0/K1 functions
which have lower errors than the original GSL implementations when
compared with arbitrary precision libraries. Through lots of random
testing he has found upper limits on the errors of each function given
by some positive constant times GSL_DBL_EPSILON. We are wondering if its
safe to change the error estimation of these functions (I don't know the
origin of the current calculations for these error estimates).

Patrick