Special Functions

Implemented in Common Lisp

The library assumes working with 64 bit double-floats. It will probably work with single-float as well. Whilst we would prefer to implement the complex domain, the majority of the sources do not. Tabled below are the special function implementations and their source. This library has a focus on high accuracy double-float calculations using the latest algorithms.

function source
erf libm
erfc libm
inverse-erf Boost
inverse-erfc Boost
log-gamma libm
gamma Cephes

Error rates

The following table shows the peak and mean errors using Boost test data. Tests run on MS Windows 10 with SBCL 2.0.10. Boost results taken from the Boost error function, inverse error function and log-gamma pages.

erf

Data Set Boost (MS C++) Special-Functions
erf small values Max = 0.841ε (Mean = 0.0687ε) Max = 6.10e-5ε (Mean = 4.58e-7ε)
erf medium values Max = 1ε (Mean = 0.119ε) Max = 1ε (Mean = 0.003ε)
erf large values Max = 0ε (Mean = 0ε) N/A erf range 0 < x < 6

erfc

Data Set Boost (MS C++) Special-Functions
erfc small values Max = 0ε (Mean = 0) Max = 1ε (Mean = 0.00667ε)
erfc medium values Max = 1.65ε (Mean = 0.373ε) Max = 1.71ε (Mean = 0.182ε)
erfc large values Max = 1.14ε (Mean = 0.248ε) Max = 2.31e-15ε (Mean = 8.86e-18ε)

inverse-erf/c

Data Set Boost (MS C++) Special-Functions
inverse-erf Max = 1.09ε (Mean = 0.502ε) Max = 2ε (Mean = 0.434ε)
inverse-erfc Max = 1ε (Mean = 0.491ε) Max = 2ε (Mean = 0.425ε)

log-gamma

Data Set Boost (MS C++) Special-Functions
factorials Max = 0.914ε (Mean = 0.175ε) Max = 2.10ε (Mean = 0.569ε)
near 0 Max = 0.964ε (Mean = 0.462ε) Max = 1.93ε (Mean = 0.662ε)
near 1 Max = 0.867ε (Mean = 0.468ε) Max = 0.50ε (Mean = 0.0183ε)
near 2 Max = 0.591ε (Mean = 0.159ε) Max = 0.0156ε (Mean = 3.83d-4ε)
near -10 Max = 4.22ε (Mean = 1.33ε) Max = 4.83d+5ε (Mean = 3.06d+4ε)
near -55 Max = 0.821ε (Mean = 0.419ε) Max = 8.16d+4ε (Mean = 4.53d+3ε)

The results for log gamma are good near 1 and 2, bettering those of Boost, however are worse (relatively speaking) at values of x > 8. I don’t have an explanation for this, since the libm values match Boost more closely. For example:

(spfn:log-gamma -9.99999237060546875d0) = -3.3208925610275326d0
(libm:lgamma    -9.99999237060546875d0) = -3.3208925610151265d0
Boost test answer                         -3.320892561015125097640948165422843317137

libm:lgamma provides an additional 4 digits of accuracy over spfn:log-gamma when compared to the Boost test answer, despite using identical computations. log-gamma is still within 12 digits of agreement though, and likely good enough for most uses.

gamma

Data Set Boost (MS C++) Special-Functions
factorials Max = 1.85ε (Mean = 0.491ε) Max = 3.79ε (Mean = 0.949ε)
near 0 Max = 1.96ε (Mean = 0.684ε) Max = 2.26ε (Mean = 0.56ε)
near 1 Max = 2ε (Mean = 0.865ε) Max = 2.26ε (Mean = 0.858ε)
near 2 Max = 2ε (Mean = 0.995ε) Max = 2ε (Mean = 0.559ε)
near -10 Max = 1.73ε (Mean = 0.729ε) Max = 0.125ε (Mean = 0.0043ε)
near -55 Max = 1.8ε (Mean = 0.817ε) Max = 0ε (Mean = 0ε)

NaN and Infinity

The lisp specification mentions neither NaN nor infinity, so any proper treatment of these is going to be either implementation specific or using a third party library.

We are using the float-features library. There is also some support for infinity in the extended-reals package of numerical-utilities, but it is not comprehensive. Openlibm and Cephes have definitions, but we don’t want to introduce a large dependency just to get these definitions.

Test data

The test data is based on Boost test data. You can run all the tests using the ASDF test op:

(asdf:test-system :special-functions)

By default the test summary values (the same as in Boost) are printed after each test, along with the key epsilon values.


Last modified 30 June 2021 : Add special-functions documentation (e4b24f3)