Recommended textbook on function approximation on computer

Question

There is a cephes math library on the Internet to provide accuracy computation of math function, e.g. sin,cos,tan,erf, gamma, lgamma, etc.

It implements the function in an intelligent way. For example, it uses different asymptotic series on the different domain when computing gamma function and switches to Stirling's approximation when the argument is too large.

Basically, the library uses all means of methods, e.g. Asymptotic series, Numerical Integration, Newton's method, Pade approximation, continued fraction.

If I want to fully understand what it is going on and study how to deduce such formula and implement such numerical algorithms, what are the recommended textbook?

For example, I want to compute an efficient and accurate specialized lgamma_exp(x) = lgamma(exp(x)) with argument in exp(x) rather than in x. There is no way except computing exp(x) first, which might overflow on the computer.

Answer 1

This is an answer to your example question.

For large values of the argument of a function, you can rely on asymptotic formulas. In the case of Gamma, you can work with the Stirling approximation

$$\log(\Gamma(e^x))\approx\log\sqrt{2\pi}+(x-1)e^x-\frac x2.$$

the relative error will be on the order of $$\frac{e^{-x}}{12}.$$

Answer 2

You should seriously consider this book:

Muller, J.-M., Brisebarre, N., de Dinechin, F., Jeannerod, C.-P., Lefèvre, V., Melquiond, G., Revol, N., Stehlé, D., Torres, S. : "Handbook of floating point arithmetic". Springer 2010.

This is a link to Springer's page for the book.

I do not remember if it covers the error-function, but it will explain many of the techniques used to write numerical libraries and it does it very well.

Answer 3

Also consider another book by Jean-Michel Muller: Elementary Functions: Algorithms and Implementations (3rd edition, 2016). Link to publisher. In fact, searching for papers etc by anyone in the groups at CNRS/ENS Lyon and/or the ARENAIRE project (e.g., at the HAL archives) would probably be beneficial.

Answer 4

In addition to the books by J.-M. Muller and co-authors already mentioned in other answers, I provided some additional pointers to relevant literature in an answer to a related question.

In the specific case of lgamma_exp, extensive literature research does not seem necessary. For the basic mathematical formulas we can rely on NIST's online resource Digital Library of Mathematical Functions (DLMF). Required knowledge about IEEE-754 floating-point arithmetic can be gleaned from the following seminal paper:

David Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic", ACM Computing Surveys, Vol. 23, No. 1, March 1993, pp. 5-48 (online)

In addition, we need to know that polynomial minimax approximations are an adequate way to approximate most simple special functions, and that tools like Mathematics, Maple, and the free Sollya tool can generate these.

From DLMF §5.7 we learn of the series $\frac{1}{\Gamma(z)} = z + \gamma z^{2} + \cdots$, which tells us that near the origin, $\Gamma(x) \approx \frac{1}{x}$. Since $\exp(38) \approx 3.2\cdot10^{16}$, for IEEE double precision, $\frac{1}{x}$ is accurate to full precision for $x < -38$, therefore lgamma_exp(x) = -x for $x < -38$. From the Stirling approximation in DLMF §5.11.1 we likewise see that for sufficiently large $x$ we have $\Gamma(x) \approx (x-1)\exp(x)$, which is accurate to full double precision for $x >38$.

The only other interval requiring attention is the vicinity of the zeros of $\log\Gamma(x))$ in the positive half plane at $x=1$ and $x=2$. In finite precision floating-point precision, for all $x < \epsilon$, $\exp(\epsilon) = 1$, so we cannot use the naive computation lgamma (exp (x)) if we desire accurate results. Here we can use polynomial minimax approximations $\mathrm{P}_{1}(\exp(x)-1)$ and $\mathrm{P}_{2}(\exp(x)-2)$. Since I am lazy, I establish the bounds of suitable intervals for these approximations numerically, using an arbitrary-precision library as the reference to determine where the error of the naive computation exceeds three ulps. I determined as suitable intervals: $[-0.188, 0.405465)$ for $\mathrm{P}_{1}$ and $[0.405465, 1.1]$ for $\mathrm{P}_{2}$.

The computations of $\exp(x)-1$ and $\exp(x)-2$ suffer from subtractive cancellation. For $\exp(x)-1$ this is easily addressed by use of the expm1 standard math function available in most computing environments. For the other interval, we can compute $\exp(x) - 2$ as 2 * expm1 (x - log(2)), where log(2) is represented to quadruple precision by a pair of double-precision constants, so the computation is accurate when $x \approx \log(2)$.

I generated the polynomial minimax approximations used in the ISO-C code below specifically for this answer using a proprietary implementation of the Remez algorithm. The polynomials are evaluated using a second-order Horner scheme for higher instruction-level parallelism which improves performance on many platforms. In preliminary testing against an arbitrary-precision library (R.P. Brent's MP library: functional but outdated; not recommended for new development), the maximum error found so far is about 3.5 ulp when the code is compiled with the Intel C/C++ compiler on an x86_64 system with strict IEEE-754 conformance (/fp:strict).

double my_lgamma_exp (double x)
{
    const double log2_hi = 6.9314718055994529e-1;
    const double log2_lo = 2.3190468138462996e-17;
    double r, s, t;
    if (fabs (x) > 38) {
        if (x < 0) {
            r = 0 - x;
        } else {
            r = (x - 1) * exp (x);
        }
    } else if ((x >= -0.188) && (x <= 1.1)) {
        if (x < 0.405465) {
            /* minimax polynomial around exp(x)=1; 2nd-degree Horner scheme */
            x = expm1 (x); // exp(x) - 1
            s = x * x;            
            r =       - 4.7091965873693646e-3;
            t =         1.9912834549094705e-2;
            r = r * s - 4.1005903342038519e-2;
            t = t * s + 5.7612615780147633e-2;
            r = r * s - 6.6582808312903363e-2;
            t = t * s + 7.1835727184955733e-2;
            r = r * s - 7.7041643077651434e-2;
            t = t * s + 8.3353766668420856e-2;
            r = r * s - 9.0949364918343004e-2;
            t = t * s + 1.0009890150433524e-1;
            r = r * s - 1.1133433459692958e-1;
            t = t * s + 1.2550968527263201e-1;
            r = r * s - 1.4404989656400469e-1;
            t = t * s + 1.6955717682057603e-1;
            r = r * s - 2.0738555102576703e-1;
            t = t * s + 2.7058080842860699e-1;
            r = r * x + t;
            r = r * x - 4.0068563438654731e-1;
            r = r * x + 8.2246703342411209e-1;
            r = r * x - 5.7721566490153287e-1;
            r = r * x;
        } else {
            /* minimax polynomial around exp(x)=2; 2nd-degree Horner scheme */
            x = 2 * expm1 ((x - log2_hi) - log2_lo); // exp(x) - 2
            s = x * x;
            r =         9.7009080515552103e-9;
            t =       - 6.6363290060543801e-8;
            r = r * s + 2.1575134587124466e-7;
            t = t * s - 4.8610795617393086e-7;
            r = r * s + 9.7870572069649216e-7;
            t = t * s - 2.0308648833487155e-6;
            r = r * s + 4.3609845401391367e-6;
            t = t * s - 9.4380453602222448e-6;
            r = r * s + 2.0510699022734775e-5;
            t = t * s - 4.4927200874090269e-5;
            r = r * s + 9.9457052285033636e-5;
            t = t * s - 2.2315458092948354e-4;
            r = r * s + 5.0966955797798817e-4;
            t = t * s - 1.1927539271577999e-3;
            r = r * s + 2.8905103294554332e-3;
            t = t * s - 7.3855510280415159e-3;
            r = r * s + 2.0580808427809637e-2;
            t = t * s - 6.7352301053207803e-2;
            r = r * x + t;
            r = r * x + 3.2246703342411304e-1;
            r = r * x + 4.2278433509846719e-1;
            r = r * x;
        }
    } else {
        r = lgamma (exp (x));
    }
    return r;
}