# How do I develop numerical routines for the evaluation of my own special functions?

This question has been cross-posted to ComputationalScience.SE here.

When performing computational work, I often come across a univariate function, defined in terms of an integral or differential equation, which I would like to rapidly evaluate (say, millions of times per second) over a specified interval to a given precision (say, one part in $10^{10}$). For example, the function $$f(\alpha) = \int_{k=0}^\infty \frac{e^{-\alpha^2 k^2}}{k+1}\ \mathrm{d}k$$ over the interval $\alpha \in (0,10)$ came up in a recent project. Now it happens that this integral can be evaluated in terms of standard special functions (in particular, $\operatorname{Ei}(z)$ and $\operatorname{erfi}(z)$), but suppose we had a much more complicated function for which no such evaluation was known. Is there a systematic technique I can apply to develop my own numerical routines for the evaluation of such functions?

I am sure plenty of techniques must be out there, as fast algorithms seem to exist for basically all of the common special functions. However, I emphasize that the sort of technique I am looking for should not rely on the function having some particular structure (e.g. recurrence relations like $\Gamma(n+1) = n\Gamma(n)$ or reflection formulas like $\Gamma(z) \Gamma(1-z) = \pi \csc(\pi z)$). Ideally, such a technique would work for just about any (sufficiently well-behaved) function I come across.

You can take for granted that I do have some slow method of evaluating the desired function (e.g. direct numerical integration) to any precision, and that I'm willing to do a lot of pre-processing work with the slow method in order to develop a fast method.

• Well-behaved functions should allow nice approximations by polynomials (piecewise). However for the desired precision this may be challenging. On the other hand, the surrounding algorithm might allow you to work with less precision most of the time, e.g., to solve $f(x)=b$ by binary search you will not need much precision as long as you are not yet close. – Hagen von Eitzen Nov 23 '15 at 19:48
• @HagenvonEitzen That is true, but the applications I work in typically demand high precision. For example, the project I mentioned above uses $f(\alpha)$ to compute the entries of some near-singular matrix $A$. The eigenvalues of $A$ are very sensitive to changes in the entries, so requiring accuracy to the nearest $10^{-10}$ is not an unreasonable demand. – David Zhang Nov 24 '15 at 1:16
• Based on my long-time professional experience in implementing software for efficient evaluation of mathematical functions in double precision or better, I doubt there is a magic bullet. In practical terms I would examine potential solutions roughly in this order: Polynomial approximation, possibly coupled with a transformation of arguments to approximately "linearize" them; rational approximation; asymptotic approximation; recurrence relation; continued fraction expansion. For more complicated functions, often two different approaches are needed and it may take months to find a good solution. – njuffa Nov 24 '15 at 7:21
• @njuffa I suspected that might be the case. Could you expand on that list of potential solutions in an answer below? I'd be glad to accept it. – David Zhang Nov 24 '15 at 22:54
• A significant question you have not answered here is whether there is any online-offline aspects to your problems. For instance, it could happen that you need to be able to evaluate your $f(\alpha)$ for any $\alpha \in (0,10)$ online to an accuracy of $10^{-10}$ in 1 ms. But you have, say, 24 hours of runtime available for offline setup of this online method. In this case you can probably construct a piecewise polynomial approximation offline and just plug into it when you need to evaluate online. This is a common situation in modern numerical problems. – Ian Nov 24 '15 at 22:56

Unfortunately, there is no single approach that will lead to robust, accurate, and high-performance implementations across the large universe of special functions. Often, two or more methods must be used for different parts of the input domain, and the necessary research and implementation work may take weeks for elementary functions and months for higher transcendental functions.

Since considerable mathematical as well as programming skill is required to produce high-quality implementations, my first recommendation would be to utilize and leverage existing mathematical libraries as much as feasible. These could be commercial libraries, such as the NAG numerical libraries or the IMSL numerical libraries by RogueWave, or open source libraries such as the GNU Scientific Library (GSL) or the math and numerics portions of the Boost library. You may also find relevant source code in online repositories such as Netlib's collected algorithms from ACM TOMS.

In practical terms, on modern SIMD-enhanced processors, the extensive use of tables is no longer advisable and (piecewise) polynomial approximations are usually the most attractive. The reason table-based methods have fallen into disfavor on high-performance processor architectures is that over the past decade the performance of functional units (measured in FLOPS) has grown much faster than the performance of memory sub-systems (measured in GB/sec). The reasoning in the following paper matches my own professional experience:

Marat Dukhan and Richard Vuduc, "Methods for high-throughput computation of elementary functions". In Parallel Processing and Applied Mathematics, pp. 86-95. Springer, 2014. (slides)

In terms of performance, the polynomial approximations benefit from the fused multiply-add operation (FMA) present in modern processor hardware (both CPUs and GPUs). This operation also helps reduce round-off errors while offering some protection against subtractive cancellation. For smallest error and best efficiency, one would want to use minimax approximation.

Commonly used tools such as Maple and Mathematica have built-in facilities to generate these. While they generate approximations that are (very close to) optimal in the mathematical sense, they do not typically account for the error incurred by representation of coefficients and evaluating operations in limited floating-point precision. The Sollya tool offers this functionality through its fpminimax command. Lastly, you can write your own approximation code which would likely be based on the Remez algorithm.

For some functions, polynomial approximations are not really practical, as too many terms would be required to reach IEEE-754 double precision. In those cases, one can chose one of two strategies.

The first strategy is to cleverly transform the input arguments, using basic arithmetic and simple elementary functions, such that the resulting function is "well-behaved" with respect to polynomial approximation. Typically such transformation tends to "linearize" the function to be approximated. A good instructional example of this approach is the computation of erfc in the following paper:

M. M. Shepherd and J. G. Laframboise, "Chebyshev Approximation of $$(1 + 2x)\exp(x^2)\operatorname{erfc} x$$ in $$0 \leqslant x < \infty$$". Mathematics of Computation, Vol. 36, No. 153 (Jan., 1981), pp. 249-253 (online)

The second approach is to use the ratio of two polynomials, that is, a rational approximation, for example in the form of a Padé approximant. The tools already mentioned can help with that; there is also a copious amount of literature published on rational approximation, which in general is a harder problem than polynomial approximation.

For special functions (as opposed to elementary functions), straightforward polynomial and rational approximations are often inaccurate and/or inefficient. They require the application of more advanced mathematical concepts such as asymptotic expansions, recurrence relations, and continued fraction expansions. Even if the use of these solves the problem mathematically, there may still be numerical issues, for example when evaluating continued fractions in the forward direction. Not surprisingly, entire books have been written on the computer evaluation of certain functions, such as Bessel functions and Mathieu functions.

In the following, I am giving a quick overview of useful literature, starting with coverage of the mathematical foundations, moving to methods suitable for elementary functions and simple special functions such as erfc and tgamma, and finally advanced methods for special functions that are harder to compute both in terms of performance and accuracy. Obviously, this can only scratch the surface, much relevant material on particular functions can be found in individual papers, for example in journals and proceedings from AMS, SIAM, ACM, and the IEEE.

Much of the literature has not yet caught up to modern hardware and software environments, in particular the presence of the FMA operation and SIMD architectures. In terms of robust computer codes for the evaluation of mathematical functions, one could wish for closer co-operation between mathematics and science on one hand, and computer science and computer engineering on the other hand. Among the works below, those by Markstein and Muller are the most advanced in this regard.

Milton Abramowitz and Irene A. Stegun (eds.), "Handbook of Mathematical Functions. With Formulas, Graphs, and Mathematical Tables". New York, NY: Dover 1972 (online version)

Frank Olver, et. al. (eds.), "NIST Handbook of Mathematical Functions". New York, NY: Cambridge University Press 2010 (online version)

A. Erdelyi, et. al., "Higher Transcendental Functions." Vol. 1-3. New York, NY: McGraw-Hill 1955

Oskar Perron, "Die Lehre von den Kettenbrüchen, 3rd ed." Vol. 1+2. Stuttgart (Germany): Teubner 1954, 1957

John F. Hart, "Computer Approximations". Malabar, FL: Krieger Publishing 1978

William J. Cody and William Waite, "Software Manual for the Elementary Functions". Englewood Cliffs, NJ: Prentice-Hall 1980

Peter Markstein, "IA-64 and Elementary Functions". Upper Saddle River, NJ: Prentice-Hall 2000

Jean-Michel Muller, "Elementary Functions. Algorithms and Implementation 3rd. ed.". Birkhäuser 2016

Nelson H. F. Beebe, "The Mathematical-Function Computation Handbook". Springer 2017

Jean-Michel Muller, et. al., "Handbook of Floating-Point Arithmetic 2nd ed.". Birkhäuser 2018

Nico M. Temme, "Special Functions. An Introduction to the Classical Functions of Mathematical Physics". New York, NY: Wiley 1996

Amparo Gil, Javier Segura, and Nico M. Temme, "Numerical Methods for Special Functions". SIAM 2007

Frank W. J. Olver, "Asymptotics and Special Functions". Natick, MA: A K Peters 1997

Jet Wimp, "Computation with Recurrence Relations". Boston, MA: Pitman 1984

A. N. Khovanskii, "The application of continued fractions and their generalizations to problems in approximation theory". Groningen (Netherlands): Noordhoff 1963

A. Cuyt, et. al., "Handbook of Continued Fractions for Special Functions". Springer 2008

• Wow. This makes me glad I'm sticking, at least for the moment, to the world of pure analysis and leaving the pesky computations to you amazing folks :) – Alan Nov 25 '15 at 7:45
• In many cases, those "pesky computations" would benefit from a fresh perspective. As Cleve Moler (of MATLAB fame) has noted, for various special functions we still fall back to Fortran codes (original or ported) written in the 1970s, while both mathematics and computer architecture have advanced since that time. – njuffa Nov 25 '15 at 8:21
• @njuffa Excellent answer! Thanks in particular for the long list of references. It is clear to me that this field is far deeper than I had first imagined. If I may ask, where have you worked to obtain professional experience in this field? I'm interested in this sort of work, and if there are career prospects here, then all the better. – David Zhang Nov 25 '15 at 20:48
• My degree is in CS, but as far as this particular expertise goes I am self-taught going back to the 1980s. I worked on floating-point units for x86 processors (which incorporate instructions for some transcendental functions) and more recently I was responsible for the standard C/C++ math library of a widely used parallel programming environment. Most people who do this kind of work are mathematicians, and it is usually in a larger context of providing all kind of numerical and performance libraries (think BLAS, LAPACK, FFT, image & signal processing). It is a small field, not that many jobs. – njuffa Nov 25 '15 at 21:15