TR1 Math Special Functions

[Ed Smith-Rowland <3dw4rd at verizon dot net>]

Here is another drop.

It is completely reworked relative to the last one.


Ed

Howdy!

This is a patch for the TR1 Mathematical Special Functions as described in
http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2005/n1836.pdf

All of the function now have useful implementations that I hope are
formatted and commented mostly correctly ;-).

This has been bootstrapped and regtested on powerpc-apple-darwin7.9.0 and on
x86_64-unknown-linux-gnu.  There are about fouri or five new fails on functions
for large arguments as described below.  These will be repaired before delivery
of 4.3!  Of course there are about 150-200 odd new passes ;-).

Design Decisions
----------------
1. In order to avoid duplication I wrote the functions as wrappers around
   template functions.  Various small and large number parameters are derived from
   std::numeric_limits parameters.

2. I have used many resources in this effort for the actual algorithms:
   * The Gnu Scientific Library: http://www.gnu.org/software/gsl/
   * The Handbook of Mathematical Functions ed. by Milton Abramowitz and Irene Stegun
     Dover Publications
   * Numerical Recipes in C by W. H. Press, S. A. Teukolsky, W. T. Vetterling,
     B. P. Flannery, Cambridge University Press (1992), 2nd ed
   * Mathematical Methods for Physicists by George Arfken, 2nd ed. Academic Press
   * Asymptotics and Special Functions by F. W. Olver, Academic Pres 1974

   I believe the algorithms implemented represent a consensus on these functions.

   In general character, these functions use series for small argument and continued
   fractions for larger arguments.  The orthogonal polynomials (Hermite, Laguerre, Legendre)
   use recursion relations by order.  The Bessel functions also use recursion in places.

   I'm very sure that the numerati will have lots to say about the algorithms though :-).

3. For the inevitable numeric constants like pi ans such I used a class analogous
   to numeric_limits that cast long double literals.  IMHO I think the C++ standards
   committee (LWG) could do worse than to consider such a replacement to math.h M_PI, etc.

4. I included a gamma function so that systems without C99 could run the special functions
   that require it.  It works pretty well but I might make it's range of validity wider
   so people could *really* use it all over.

Goals before delivery of gcc-4.3
--------------------------------
1. The only function I feel really bad about is the hypergeometric function.
   We need a better implementation there.

2. The confluent hypergeometric function is a little wobbly in places too
   (although it is in better shape than the hypergeometric function).

3. The Bessel functions need asymptotics on argument and order.  Some of these
   are coded but need to be integrated with good switchover points.

4. The associated and regular Laguerre polynomials fall apart at larger orders.

5. More testcases - especially ones from practice (Hydrogen wavefunctions,
   electromagnetic wave scatttering from cylinders and spheres, ...)

6. I would like to strengthen the gamma function implementation.  I'm using
   it as a stopgap plug for systems without C99 math.  But it would be nice
   to templatize the Lanczos algorithm and make it first class.

Longer Term Goals
-----------------
1. There are some in the C++ standards committee who note that the choice of functions
   lean towards the physical sciences and engineering.  Some suggest that functions
   used bt statisticians and bilologists would broaden the appeal of this part of the library.

2. There are functions that would be useful (some even have partial implementations in
   the library) such as Airy functions, Fresnel function, sine and cosine integrals,
   Jacobian elliptic functions.  These would be good extensions too.

3. The Gnu Scientific Library has a sophisticated error propagation system.  Returning the error
   of the functions is not part of the TR1 specification.  However, it might be useful for developers
   in the underlying template functions.  We could also offer this as an extension.

4. Many math libraries offer versions of special functions that return a vector of values
   for a given argument over a range of orders.  This is most true for orthogonal polynomials.
   This is because one often needs a range of orders for polynomial expansions.  Also, these
   polynomials are computed by iteration on order so to get a high order you compute all the
   intermediate ranges anyway.  This would be a nice extension.  Or perhaps we could make
   the ones we have into a functor to amortize subsequent calls somehow.


Have fun,
Ed Smith-Rowland

Attachment: tr1_math44.diff.bz2
Description: Binary data