How to compute a special function with near machine-precision accuracy

Abstract

Based on my experience as developer and maintainer of some numerical open-source libraries (libcerf, libkww, libformfactor), I will explain key concepts for writing code that computes a special function or integral with high accuracy and high speed.- Choose different numerical algorithms for different argument regions.- Don't be afraid of divergent series or ill-conditioned recursions. - Confine Chebyshev fits to small subregions. - Use code instrumentation and bisection to ensure continuity where the algorithm changes. - Beware of literature that is only concerned with truncation. Near machine precision, cancellation is the bigger problem. Visualization may reveal the difference. - Never rely on non-standard facilities. "Long double" makes no sense if it is not longer than "double"? Tell Apple. - Generate test references and hard-coded coefficients with high-precision scripts (e.g. mpmath based). - Be graceful with relative accuracy measures near zeros and for real or imaginary parts of complex numbers.- Don't trust any performance measure you haven't tweaked yourself. Take caching into account