by Gaston H. Gonnet, ETH, Informatik
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 | 1304 | 0.57432 | -0.875010530199998859 | 2.63625389482078987 |
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.
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 |
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.
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.
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.
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.