Features
FeaturesWhen requirements are tight, we must carefully examine all potential sources of error. Make sure your math library isn't the weak link in the chain.
Introduction
The purpose of the software described in this article is to remove some of the mystery surrounding the mathematical functions in your runtime library. Most of us consider these functions to be black boxes; that is, they work because the compiler manufacturer says they work. When challenged to show that they work, what can we do? With this article you will find software that is designed to shed light on some of these black boxes.
"Knowledge is power " is the tenet that guided development of these programs. It is common practice for users to take the mathematical functions in a runtime library at face value. After all, the square root function calculates the square root of a number doesn't it? Exactly, too, right? But given accuracy information about a function, how do you verify its authenticity? If you can't demonstrate a function's performance then you must accept its output on faith. If you're like me, that probably makes you feel more than a little uncomfortable.
If you are using mathematical functions in your runtime library you should know how well they work. It is not wise to simply accept what your compiler manufacturer gives you on faith. Libraries are written by humans, and humans make mistakes. Besides, complex accuracy specs are easily misinterpreted by manufacturer and customer alike. Despite a manufacturer's best intentions, you may never be sure what they're giving you unless you test it yourself. If you would like to know how well mathematical functions in your runtime library work, the programs presented in this article will supply much of that knowledge.
Also, if you develop mathematical functions yourself, these test programs will supply the performance statistics that you need to verify that your functions perform correctly, or that your programs work as advertised.
Description of Test Package
There are two levels of test programs in this package. The first level of programs tests the standard-precision mathematical functions in your runtime library. These programs measure performance of each function against a baseline function that produces results at greater precision.
The second level of test programs generates performance figures for the baseline functions.
The ideas of Cody and Waite [1] and Plauger [2] played a large role in the development of this test package. The architecture of the test programs is primarily that of Plauger.
Following Cody and Waite, each test program conducts its testing in two stages. In the first stage, the function under test is called using random arguments over two or more ranges of its input. Except for the square root function, these test programs produce arguments that are uniformly distributed. Logarithmically distributed variates are supplied to the square root function. The test program gathers statistics for maximum magnitude of the relative error (MRE), and the root mean square error (RMS E). MRE approximates worst-case behavior, while RMS indicates overall quality [1, page 13]. For MRE, printouts consist of MRE and its base two logarithm, the argument at which MRE occurred, and an expression of the number of bits lost. For RMS E, printouts consist of RMS E, its base two logarithm, and an expression of the number of bits lost.
The second stage consists of calling the function under test with specific arguments that "might cause or detect trouble" [1, p. 15]. During this stage, each function is called with arguments at the extreme ranges of its domain, arguments that test even or odd symmetry, and those that test out-of-bounds conditions.
Each test program, both first-level and second-level, uses information about the particular arithmetic in use. The first level (i.e. tests of run-time library functions) uses information obtained from standard header file float.h. The first level program calls function GetMachar (Listing 4) to extract this data from float.h and to derive other data not provided in float.h. GetMachar obtains the following data:
Name Description 1. Radix Radix of floating point representation 2. FracDigs Number of base Radix digits in floating point fraction 3. Rounds TRUE if addition rounds 4. EpsNeg Smallest x such that 1.0 - x != 1.0 5. Eps Smallest x such that 1.0 + x != 1.0 6. Max Maximum number 7. Min Smallest unnormalized positive number 8. ExpDigits Number of digits used for exponent 9. MaxExp Maximum exponent 10. MinExp Minimum exponent 11. NumGuard Number of multiplier guard digitsThe programs that test the baseline functions call function GetExtCh to extract characteristics of the arithmetic in the extended-precision library. Both GetMachar and GetExtCh owe their origin to Plauger's code in elefunt, which in turn comes from Cody and Waite.
Tests of Runtime Library Functions
I have written programs that test the following standard functions:
1. asin, acos 2. atan, atan2 3. exp 4. log, log10 5. pow 6. sin, cos 7. sinh, cosh 8. sqrt 9. tan, cot 10. tanhExcept for cot, this list matches the ANSI C set.
I have also written test programs for the following functions:
11. asinh Inverse hyperbolic sine and cosine 12. acosh Inverse hyperbolic cosine 13. atanh Inverse hyperbolic tangent 14. cbrt Cube root function 15. j0, j1, jn Bessel functions of the first kind 16. y0, y1, yn Bessel functions of the second kindFollowing the example set by P.J. Plauger's elefunt [2] , each test program can be compiled in three precisions -- float, double, and long double. As you can surmise from the foregoing lists, there are 16 primary test source files. This means that there are 48 executable functions in all. The naming convention is:
xv<name>P.exewhere <name> has up to five characters and P is one of f, d, or l, signifying float, double, or long double respectively. For example, the processor that verifies the double-precision square root function is xvsqrtd.exe. The xv stands for "extended precision verification."
The methods that I use to test the library functions depart somewhat from the methods of Cody and Waite. Their methods depend largely on self-tests with few outside agents involved. The test programs that I have written compare the output of the library functions with the output produced by a set of extended-length arithmetic functions. As pointed out by Cody and Waite [1, p. 11], this is the ideal method for conducting accuracy tests. The extended-length functions were written by S. L. Mosier [3] .
Source Code Configuration
The following dependency list serves an example of how the test programs are built. The list refers to one of the test programs for the square root function, namely xvsqrtd.exe.
xvsqrtd.exe xvsqrtd.cxp xvsqrt.h xverfun.h qfloat.h machars.h xtstdefs.h dprotos.h getmachx.hThe indentations indicate #include level. As you might suspect, the hierarchy is similar for xvsqrtf.exe and xvsqrtl.exe. The extension cxp stands for "C++ extended precision."
Not all the .h files are header files in the sense of providing prototypes and data definitions. Three are outright program files:
xvsqrt.h - contains main and driving functions xverfun.h - calls driving functions and performs statistical calculations getmachx.h - returns machine characteristics for the arithmetic under testI describe these and the top-level program file, xvsqrtd.cxp, next.
Listing 1 shows program file xvsqrtd.cxp. To show the similarity of the parent program files, I have included the contents of xvsqrtf.cxp and xvsqrtl.cxp. Except for differences in comment lines, the critical difference is in the #define statement at line eight of each file.
The principal user of this definition is header file machars.h, shown in Listing 5. For example, when symbol DBL is defined in xvsqrtd.cxp, lines 20-21 in machars.h are activated. The synonym for TYPE is set to double, and the NAME() macro for double-precision functions is defined. The typedef for TYPE is quickly reflected at lines 33-36, which are part of the structure that will hold machine characteristics for double-precision arithmetic.
All three program files shown in Listing 1 include file xvsqrt.h, which is shown in Listing 2. xvsqrt.h includes a main at line 48. Of immediate interest are the uses of data type qfloat at lines 11, 16 and 30, the TYPE statement at line 18, and the NAME() macro at line 25. Since symbol DBL is #defined in parent program file xvsqrtd.cxp, TYPE is a synonym for double and macro NAME(sqrt) yields sqrt. If the parent program file were xvsqrtl.cxp, LNGDBL would be #defined, TYPE would be a synonym for long double, and NAME(sqrt) would yield sqrtl. Data type qfloat holds extended precision floating-point numbers.
As stated at the beginning of this section, testing is performed in two stages. This is shown clearly in Listing 2. The first stage, when random arguments are thrown at the function under test, extends from line 64 to line 76. Printouts from this stage are produced by subroutine XVerFun, shown in Listing 3. The second stage, during which the function under test is driven with special values begins at line 78. The following special values are used:
Item Special Value 1. 0.0 2. 1.0 3. 1.0 - EpsNeg 4. 1.0 + Eps 5. Min 6. Max 7. Integers from 1 to 1000 8. -1.0where, as described previously, EpsNeg, Eps, Min, and Max are floating point numbers for the arithmetic precision under test:
Name Description EpsNeg Smallest x such that 1.0 - x != 1.0 Eps Smallest x such that 1.0 + x != 1.0 Max Maximum number in the arithmetic Min Smallest unnormalized positive numberItem number 8 is called a boundary value.
Machine characteristics for the arithmetic in use are produced by function GetMachar. This function returns a structure of type MacharStru, declared as a typedef in machars.h (Listing 5) as MACHAR_STRU. The function is called at line 59, Listing 2, and an informative printout is issued at line 61.
Tests of the Baseline Functions
To give you some confidence that the Moshier functions produce decent baseline values, I have written "confidence builders" -- test programs for the functions that correspond to those numbered 1-16 in the lists at the beginning of the previous section. The code in programs corresponding to those numbered 1-10 follows closely the code from Plauger's elefunt. The code in the remaining six follows the code that I wrote to verify the corresponding standard-precision functions.
Like the Cody and Waite programs, these confidence builders rely largely on self tests, although some produce baseline values from Taylor series expansions.
There are 16 extended-test source files. The naming convention is:
xt<name>.exewhere <name> has up to five characters. For example, the program that tests the extended-precision square root function is xtsqrt.exe. The xt stands for "test of extended precision function."
Development Language
All test programs are written in C++. This includes both the verification programs (tests of standard-precision library functions) and the confidence builders. Without C++ the code would be very difficult to read. Moreover, I would not have been able to translate Plauger's elefunt code nearly so easily. In some cases, it was simply a matter of changing the target data type to qfloat, a class defined in C++ wrapper file qfloat.h. This file is largely the brainchild of S. L. Mosier.
There is nothing fancy about the C++ nature of the code. One new data type is introduced, a class called qfloat. An object of this class is used to hold extended-precision floating-point numbers.
Listing 6 shows the code fragment from qfloat.h that defines the class qfloat.
Some Test Results
Of the 16 tests listed previously in the section Tests of Runtime Library Functions, I have selected double-precision and long double precision versions of the tests of exp, pow, and sin/cos for show and tell. I chose these tests because 1) these functions presented some of the greatest challenges to the extended precision programs and 2) runtime versions of type float are not available in the Microsoft C/C++ 7.0 library.
I ran these tests on three C runtime libraries:
1. Microsoft C/C++ 7.0
2. Moshier's double and long double libraries
3. P.J. Plauger's Standard C library
Recaps of the test runs on the functions in these libraries are shown in Figures 1, 2, and 3 respectively.
I also ran corresponding tests on the extended-precision versions. The recaps appear in Figure 4. These tests are meant to show that there is enough precision in the results to represent the true function values to the number of bits required.
The following definitions are useful for understanding the column headings in the figures:
(Log 2) MRE
This column shows the base 2 logarithm of the Maximum Relative Error. That is, 2val yields the MRE. MRE is calculated as follows:
Bits Lost
The calculation for bits lost is:
Val = Bits in Significand + log2 MRE if Val > 0 Bits Lost = Val else Bits Lost = 0(Log 2) RMSE
This column shows the base 2 logarithm of the Root Mean Square Error. That is, 2val yields the RMS E. RMS E is calculated as follows:
Bits Lost
The calculation for bits lost is:
Val = Bits in Significand + log2( RMS E ) if Val > 0 Bits Lost = Val else Bits Lost = 0Test Range
This column contains one or two lines per test that recap Range. If the function was called with two random variables, for example pow(x,y), this column shows two ranges, one each of x and y.
Individual Tests
Microsoft, Figure 1
The test results for Microsoft's double precision functions are the best of the lot. The same cannot be said for Microsoft's long double precision. For example, in the tests of expl, the second and third tests show that there is a problem in this function. To be sure, the ranges of arguments are larger than for double precision and they should be. When run with the same argument ranges, expl produces answers of about the same quality as those for its double precision counterpart, exp.
Moshier, Figure 2
The results for both Moshier's double and long double precision functions are as good or better than those reported in Cody and Waite. This includes the fourth tests of pow and pow1 where the number of lost bits is greatest.
Plauger's Standard C, Figure 3
The results for both Plauger's double and long double precision functions are as good or better than those reported in Cody and Waite. This includes the fourth tests of pow and pow1 where the number of lost bits is greatest. The errors compare favorably with Moshier's but are, in general, slightly higher.
Moshier's Extended-Precision Functions, Figure 4
As stated previously, the tests on these functions are meant to show that the functions behave as expected. That is, they produce answers that are of high enough quality to be used as baseline functions. The errors in these functions are the highest of the standard functions -- those listed in 1-10 of the section Tests of Runtime Library Functions. Nonetheless, these functions retain at least 128 bits of accuracy. o
References
[1] W. J. Cody and W. M. Waite. Software Manual for the Elementary Functions (Prentice-Hall, Inc., Englewood Cliffs, NJ, 1980).
[2] P. J. Plauger. Elefunt Code (unpublished, adapted with permission from Cody & Waite [1980]).
[3] S. L. Moshier. Methods and Programs for Mathematical Functions (Ellis Horwood Ltd., Chichester, 1989).
[4] P. J. Plauger. The Standard C Library (Prentice-Hall, Inc., Englewood Cliffs, NJ, 1992).
K. B. Williams received a B.A. in Mathematics from California State University at Sacramento in 1955. He is a veteran of 40 years in the scientific applications programming trenches and is the founder of King Baker Enterprises, in Stillwater, Oklahoma. He consults in applied numeric techniques. He can be reached at 505 Banyan Way, Melbourne Beach, FL 32956, or via e-mail at Kbwms@aol.com.
