You can put a lot of work into writing math functions and still get embarrassed when they go out into the field. It's very hard to maintain something resembling full precision over the entire set of valid input values (the domain). It's equally hard to ensure that the function avoids spurious intermediate overflows and underflows over the entire set of representable output values (the range). The trick lies in testing the math functions, often repeatedly, as you develop them.
If you think writing the functions is hard, try your hand at writing tests. A good test calculates "true" values with exquisite accuracy. It can afford to run much slower than production code, or be much larger, but it still has to get really correct results. The test must also be sure to check the function at a sufficient number of representative values to discover any systematic errors. Typically, a math function accepts an astronomical number of input values, so exhaustive testing must be replaced by judicious sampling. Finally, the test must poke at end points, special cases, or regions of known difficulty to be sure the math function doesn't cut corners.
Fortunately, a good set of math tests has been around for a long time. William J. Cody, Jr. and William Waite [1] wrote the "elefunt" tests for the elementary math functions many years ago. They worked in Fortran, with a rather different mix of floating-point architectures. But the work they did was so comprehensive and so elegant that it has stood the test of time. You will even find a translation of the elefunt tests in the Pascal validation suite used widely in the 1980s.
I translated elefunt into Standard C a while back, in the process of refining the math functions in my Standard C and draft Standard C++ libraries. When Tim Prince offered to develop the full set of quad-precision math functions presented here, I happily made my version available to him for testing. He has clearly used them to good advantage.
P.J. Plauger