# Error analysis of system mathematical functions

This is a relatively demanding test for the accuracy of mathematical system functions like, sin, exp, lgamma, etc. It consists of a group of programs, one for each function, which can be run on a target machine and will produce a short report about the precision of the function in that machine.

There are two main goals for producing this software: first to understand how reliable the numerical evaluations of the system you use are, and secondly to provide some encouragement to the developers to improve some areas with problems.

Currently the programs are in C, but it is not difficult to extend them to other languages. It is also possible to use various floating point systems. The tests are geared to test double precision functions. Neither float nor long doubles are considered.

The function is basically tested over various inputs, for which we have precomputed the result to higher precision. These inputs are difficult, in various senses. First, a list of known problematic values are tried. Secondly, values which will evaluate to problematic values are also tried. Finally, ranges of values are explored randomly, but for each value, nasty neighbours are computed. A nasty value is one for which the exact result will be extremely close to 0.5 ulp. This is done with a cute application of the LLL lattice reduction algorithm. By testing with these numbers we often find the weakest approximations.

The typical summary result for a function looks like:

function tests max err in ulps argument (max err) value (max err)
acos 13040.57432 -0.8750105301999988592.63625389482078987

where the function tested is `acos` (arc cosine), the test had 1304 points, and the worst error was 0.57432 ulps for the argument -0.875... which gave 2.636... as a result. Clicking on the function name will show the actual output of the program for acos, which contains some more information. In particular, the 20 worst errors are shown. This is very helpful to find why these errors may be happening. (See for example j0 where all the problems happen very close to the zeros.).

For a result which is stored in a double precision variable (with no extra bits or extended precision), it is impossible to achieve a result better than 0.5 ulp (for most reasonable functions). In practical terms, a maximum error of 0.51 ulp is a superb result, a maximum of 0.6 ulp is a very good result, a maximum of 1 ulp is getting to be a bit weak.

# Results available

Currently results are available for:

hardware operating system compiler
Intel Pentium or AMD PC Linux `gcc -O9 -lm`
Sun - Sparc Solaris `cc -xO4 -lm`
DEC - Alpha OSF1 `cc -O4 -lm`
SGI Octane Irix64 `cc -O3 -lm`

# Model of computation

The model of error in the computation is the standard one for numerical analysis. The inputs to the functions (floating point numbers) are considered exact rationals. The result of the function computation is compared against the value of the function at the exact rational. The error is measured in ulp (units in the last place), that is the values of the results are multiplied by a power of the base such that an error in the last significant place is an error of 1.

Some functions in some systems suffer from very large ulp errors for arguments which produce values close to zero. This is still a serious flaw in our view, but a bit more forgivable by others.

We have included some special functions for reasons explained below.

name mathematical reason
INV 1/x Testing hardware division
Pix Pi x Testing hardware multiplication
pow2_x 2x the `pow` function has two arguments, test one
powx_275 x2.75 the `pow` function has two arguments, test the other
Gamma exp(lgamma(x)) the lgamma hides too much error

The computations are done in two ways, first they are done and stored in a double (and we try to confuse the compiler so that it cannot optimize it away). Secondly they are done together with the subtraction of the largest part of the result. The second computation will reveal if the accumulator of the computer has more precision than the regular doubles. This is reported in the full results of each function.

# Surprises

There are some pleasant surprises and unpleasant ones. There are also some hardware characteristics, which will inevitably show up in the results.

The zeros of functions cause havoc to most systems. This is well known to be a problem and it is easy to solve by a special approximation around the zero. This may need a careful argument reduction, for example for trigonometrics. The zeros of lgamma on the negative range and the zeros of the Bessel functions caused havoc uniformly.

``Chapeau'' (as the French would say) for solaris, without any extra precision it gets a superb collection of results.

I was surprised to find that some hardware does not round correctly for division and multiplication. This is the case of the Pentium/AMD, which get very small ulp > 0.5 both for multiplication and division. It was pointed out to me by Rolf Strebel and Brad Lucier that this problem is really due to double rounding, first at 80 bits then at 64. It can be corrected by setting the FPU control word to FPU_DOUBLE using the _FPU_SETCW macro. With this setting, the results are much worse for the functions, while correct for mult/div/sqrt.

The results for Pentium/AMD and Linux, outside the pitfalls, show a very uniform good behaviour, which can be attributed to the extra bits in the accumulator.

# How it is all done

If you want to try these tests on your favourite machine, you can download the data and compile and run all the programs under the directory `C`.

If your system has `make`, then your life will be much easier, just add the entry at the top of the makefile in the same style as the others, name your machine appropriately (say xxx) and run:

```    setenv ARCH xxx
make clean
make xxx/summary.html
```
and your summary file will be completely built.

In all cases you should check that your system is properly described in the file `float.M`. A superb tool for finding about your system is enquire. It is very portable, will run smoothly and will give you lots of information about your hardware and C compiler. The file float.M can be reconstructed from the values obtained from enquire. At this point you will need to have Maple to run the program that generates all the new test files for your particular hardware. In this case, `make` will still suffice (but you will have to be more patient).

If you want to include a new function, this is relatively trivial to do. Follow the examples of any of the `functname.M`. The main function to be run is `GenerateTestFunction` which has five arguments: Maple's name of the function, language (C so far), C name of the function, set of problematic values, set of test ranges (you should give at least one range). Then add your function to the `results` part of the makefile and run it. That's all.

The programs to be executed are a sandwich of three parts: a header (`C/header.h`) file, a set of constants generated by the Maple program `GenerateTestFunction` and a trailer (`C/trailer.h`) part which contains the executing program. The header and the trailer are included in the Maple generated program, which is then self-contained. To learn more details just read the comments of the programs.

To verify a single value, there is a Maple function called `Check` in the file Check.M. This is useful when you suspect about a particular result and would like to know exactly how are the argument and results represented. Check will print the arguments and values as given and as they will (should) be represented in the particular hardware. It also does a backwards error analysis.

The technique to find nasty values to evaluate using the LLL lattice reduction algorithm is described in this note.

# Contributions, copyright, distribution, disclaimers

We welcome contributions, i.e. results of running the programs in other hardware, fixes, extensions, additional comments, or any contribution you feel is necessary to make. Please mail them to me.

All the material that is under this website is unrestrictedly available to the world. Copyright © Gaston H. Gonnet.

This program/data is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

Neither Gaston Gonnet, nor the ETH (Eidgenössiche Technische Hochschule Zürich) nor any of its members assert the accuracy or validity of the results presented here, nor are responsible for them in any thinkable way. The results were obtained using our best efforts from our existing computers, but we cannot be responsible for any errors or omissions. In case of doubt, readers are encouraged to reproduce the results obtained here by their own means.

The entire directory can be downloaded from a compressed tar file.

# ----

with best wishes, and hoping that this may be useful for you, Gaston Gonnet