P.J. Plauger Dinkumware, Ltd. pjp@dinkumware.com
Four of the "special math" function families added with TR1 are
cylindrical Bessel functions. Each has a parameter nu that
describes its order, as in:
// regular modified cylindrical Bessel functions: double cyl_bessel_i(double nu, double x); float cyl_bessel_if(float nu, float x); long double cyl_bessel_il(long double nu, long double x); // cylindrical Bessel functions (of the first kind): double cyl_bessel_j(double nu, double x); float cyl_bessel_jf(float nu, float x); long double cyl_bessel_jl(long double nu, long double x); // irregular modified cylindrical Bessel functions: double cyl_bessel_k(double nu, double x); float cyl_bessel_kf(float nu, float x); long double cyl_bessel_kl(long double nu, long double x); // cylindrical Neumann functions; // cylindrical Bessel functions (of the second kind): double cyl_neumann(double nu, double x); float cyl_neumannf(float nu, float x); long double cyl_neumannl(long double nu, long double x);
J. Marraffino and M. Paterno have proposed in N1665 that we
define these functions only for nonnegative nu and
x. That's a help to implementors, but it doesn't
go far enough. We at Dinkumware have put considerable effort
into writing top quality versions of all the special math
functions. Much of that effort has gone into the cylindrical
Bessel functions. With a heroic effort, we can now generate
accurate results even for the largest possible values of
x. But for large values of nu
(greater than around 100) we have yet to develop adequate
techniques. Nor have we found any existing implementations
that have done better.
The problem stems from the fact that, for other than very
small or very large values of x, computational
techniques involve some form of recurrence on nu,
in steps of 1 from smaller values. For large nu,
these techniques either lose accuracy or take unacceptably
long time. It is possible, in principle, to express a
Bessel in terms of Airy functions whose arguments combine
nu and x, but our testing shows
that existing attempts to do so result in very large errors.
Nor have we yet found a way to do better.
A related problem exists with a number of other functions. These take unsigned integer arguments, as in:
// associated Laguerre polynomials: double assoc_laguerre(unsigned n, unsigned m, double x); float assoc_laguerref(unsigned n, unsigned m, float x); long double assoc_laguerrel(unsigned n, unsigned m, long double x); // associated Legendre functions: double assoc_legendre(unsigned l, unsigned m, double x); float assoc_legendref(unsigned l, unsigned m, float x); long double assoc_legendrel(unsigned l, unsigned m, long double x); // Hermite polynomials: double hermite(unsigned n, double x); float hermitef(unsigned n, float x); long double hermitel(unsigned n, long double x); // Laguerre polynomials: double laguerre(unsigned n, double x); float laguerref(unsigned n, float x); long double laguerrel(unsigned n, long double x); // Legendre polynomials: double legendre(unsigned l, double x); float legendref(unsigned l, float x); long double legendrel(unsigned l, long double x); // spherical Bessel functions (of the first kind): double sph_bessel(unsigned n, double x); float sph_besself(unsigned n, float x); long double sph_bessell(unsigned n, long double x); // spherical associated Legendre functions: double sph_legendre(unsigned l, unsigned m, double theta) float sph_legendref(unsigned l, unsigned m, float theta) long double sph_legendrel(unsigned l, unsigned m, long double theta) // spherical Neumann functions; // spherical Bessel functions (of the second kind): double sph_neumann(unsigned n, double x); float sph_neumannf(unsigned n, float x); long double sph_neumannl(unsigned n, long double x);
Once again, the computational time increases linearly with the size of the unsigned integer parameters. Typically, accuracy also deteriorates for large enough values.
We thus propose, as a matter of practical necessity, an upper limit on the value of any of these parameters that an implementation is required to support. It may accept larger values, and even produce good results,, but it is not required to do so. We suggest that is upper limit be 128 in all cases, for simplicity and uniformity.