By Cleve Moler, MathWorks

We've just added the code for the *tetragamma* function to MATLAB; it will be in the next release. I suspect that most of you have never heard of the *tetragamma* function. I hadn't, until my colleagues working on the Statistics and Machine Learning Toolbox needed the *digamma* function to compute maximum likelihood fits to binomial probability distributions.

You may be familiar with the *gamma* function. It's defined by an infinite integral, but it is known best as the continuous function that interpolates the integer factorial function

Γ(*n*) = Γ(*n*-1)!

The *digamma function*, which is also called the *psi* function, is the logarithmic derivative of the gamma function

Ψ(*x*) = *d*[ln Γ(*x*)]/*dx* = Γ'(*x*)/Γ(*x*)

If you take higher order derivatives, you get functions named *trigamma, tetragamma, pentagamma,* and so on.

All of this is explained in chapter 6 of A&S, the definitive reference for these matters, "Handbook of Mathematical Functions," by Milton Abramowitz and Irene Stegun. A&S was first published by the National Bureau of Standards in 1964. A Web- and CD-based successor to A&S, the Digital Library of Mathematical Functions, is under development at the National Institute for Standards and Technology, see http://dlmf.nist.gov.

Much of A&S is devoted to tables that we don't need anymore, except to check the results computed by MATLAB and other mathematical software. But the algorithms we use in modern software for special mathematical functions are derived from techniques used in hand computation over 50 years ago. There are too few people today who know how to write robust software for computing these functions.

A predecessor to Abramowitz and Stegun is an amazing book by Jahnke and Emde, "Tables of Higher Functions", published in Germany in 1909. (Reprinted by Dover in 1945 and by McGraw Hill in 1960.) Jahnke and Emde, and their students, computed tables of dozens of different functions. They also created, again by hand, intricate 2- and 3-dimensional graphs showing the behavior of these functions in the complex plane. I've included two of my favorite graphs here. One shows the complex valued gamma function. The height of the surface is the modulus, or absolute value, and the contour lines are the modulus and phase. Since

Γ(*n*) = Γ(*n*+1) /*n*

you can see that this function has poles at the negative integers. The other graph shows the Hankel function, or complex valued Bessel function. I've made MATLAB versions of these graphs and added color and lighting, but it's not so easy to recreate the labeling and other nuances of the hand-drawn work.

To appreciate the subtleties involved in the computation of special functions, consider a more basic and familiar function, sin(*x*). You might think you could use the power series

sin(^{x}) = *x* - ^{x3}/3! + ^{x5}/5! - ...

After all, this series converges for all *x*. But it turns out to be both inefficient and inaccurate. Here is a toy MATLAB function implementing the power series approach.

function s = sin(x) s = 0; t = x; n = 1; while s+t ~= s; s = s + t; t = -x.^2/((n+1)*(n+2)).*t; n = n + 2; end

Note that there is no `x^n or factorial(n)`

. Each term `t`

is computed by multiplying the previous term by the appropriate factor. The number of terms required is determined by roundoff error. Eventually, `t`

becomes so small that it is numerically negligible compared to the accumulated sum `s`

.

This function works pretty well if *x* is small enough. Take *x* = π/4 = .785398.... In this case, the series starts with

.785398... - .080746... + .002490... - .000037... +...

The first term is the largest. It takes just nine terms to compute sin(*x*) = .707107... to full double precision accuracy.

Now try *x* = 7.5π = 23.5619… . Since *x* is an odd multiple of π/2, sin(*x*) should be -1. But the power series starts out with

23.5619... - 2180.13... + 605165.9... - 799922...

The terms get too big before they start to get small. The biggest term is *x*^{23}/23!, which is larger than 10^{9}. We lose 9 of the 16 digits available with double precision floating point arithmetic. The final computed sum is -0.99999997360709. If that weren't bad enough, our function takes a long time because it uses 47 terms.

Of course, it gets worse. The power series always converges, but as *x* increases, some of the intermediate terms become so large that the computed sum eventually loses all accuracy. The loss of accuracy results from the magnitude of the terms, not their number, even though that increases as well. Fortunately, we can make use of an important property of sin(*x*), its periodicity,

sin(*x*+2π) = sin(*x*)

So, if | *x* | > π, subtract an integer multiple of 2π to bring it into the range where the power series has acceptable numerical properties. This is called *argument reduction*. We can actually combine periodicity with other trig identities to reduce the argument to the interval | *x* | < π/4. On the reduced interval, sin(*x*) is approximated by a polynomial of degree 13 with only odd powers of *x*

sin(*x*) ≈ *x* - *c*_{1}*x*^{3} + *c*_{2}*x*^{5} + … + *c*_{6}*x*^{13}= *p*(*x*)

The six coefficients are close to, but not exactly equal to, the power series coefficients 1/3!, 1/5!, …, 1/13! . They minimize the maximum relative error, *| (sin(x) - p(x))/sin(x) |,* over the interval. Six terms are enough to make this approximation error less than 2^{-52}, which is the roundoff error involved when all the terms, and the sum, are less than one.

The resulting algorithm for approximating sin(*x*), as well as other trig functions, exp(*x*) and

log(*x*), is carefully implemented in *fdlibm*, a "Freely Distributable Math Library," developed at Sun Microsystems by K. C. Ng and others, (see www.netlib.org/fdlibm). We use *fdlibm* in MATLAB for all the machines we support. This ensures uniform behavior, particularly for exceptional conditions.

Intel microprocessor architecture includes instructions for evaluating polynomial approximations to sin(*x*) and cos(*x*) on a reduced interval, but MATLAB doesn't use them because the Microsoft optimizing C compiler does not provide a complete argument reduction.

This brings up a delicate and controversal topic: what is the correct value of the MATLAB expression

`sin(pi)`

We all agree that sin(π) is zero. But the quantity that MATLAB denotes by `pi`

is a floating point number that is not exactly equal to the theoretical mathematical quantity denoted by π. It turns out that, π is about 1/4 of the way between pi and the next larger floating point number, `pi+2*eps`

. In fact,

*x* - `pi`

≈ 1.2246·10^{-16}

Near π, sin(*x*) is very nearly linear, with a slope of -1. So `sin(pi)`

should be about 1.2246e-16. The exact value reveals information about the value of π used internally by the math library for argument reduction. Many people expect `sin(pi)`

to be zero, but then we have to worry about `sin(k*pi)`

for increasingly large integer `k`

.

This approach to computing sin(*x*), except for the details about π, carries over to other mathematical functions. The general outline is:

- Transform the argument to one or more intervals where the function can be accurately approximated.
- Evaluate a polynomial or rational approximation for the transformed argument.
- If necessary, adjust the result to account for the initial argument transformation.

The argument reduction formulae for the log and exponential functions allow the reduction to an interval near *x* = 1.

exp(*k x*) = (exp(*x*))^{k}

log(2^{k} *x*) = *k* log(2) + log(*x*)

The algorithm for the real gamma function shows the complexity of the more specialized algorithms.

If 0 ‹ *x* ≤ `eps`

, Γ(*x*) ≈ 1/*x*.

If `eps`

‹ *x* ‹ 1, use Γ(*x*) = Γ(*x*+1)/*x*.

If 1 ≤ *x* ≤ 2, compute a rational approximation.

If 2 ≤ *x* ≤ 12, use Γ(*x*+1) = *x*Γ(*x*) several times.

If 12 ‹ *x*, use an approximation related to Sterling's asymptotic series,

log(Γ(*x*)) ≈ (*x* - 1/2) log(*x*) - x + 1/2 log(2π) + c_{1}/*x* + c_{2}/*x* ^{2} +...

The decisions that go into these algorithm designs-the choice

of reduction formulae and interval, the nature and derivation of the approximations-involve skills that few mathematicians have mastered. The algorithms that MATLAB uses for gamma functions, Bessel functions, error functions, Airy functions, and the like are

based on Fortran codes written 20 or 30 years ago by W. J. Cody at Argonne National Labs and Don Amos of Sandia National Labs. (You can see code in the MATLAB directory `toolbox/matlab/specfun`

.) Cody and Amos are now retired. I hope we have people around today to fill their shoes.

Published 2002