Special Functions
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 |
| incomplete-gamma | Boost |
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ε) |
incomplete-gamma
See boost incomplete gamma documentation for notes and error rates.
lower
| Data Set | Boost (MS C++) | Special-Functions |
|---|---|---|
| small values | Max = 1.54ε (Mean = 0.439ε) | Max = 3.00ε (Mean = 0.516ε) |
| medium values | Max = 35.1ε (Mean = 6.98ε) | Max = 10.00ε (Mean = 0.294ε) |
| large values | Max = 243ε (Mean = 20.2ε) | Max = 20ε (Mean = 0.613ε) |
| integer and half-integer | Max = 13ε (Mean = 2.97ε) | Max = 3ε (Mean = 0.189ε) |
upper
| Data Set | Boost (MS C++) | Special-Functions |
|---|---|---|
| small values | Max = 2.26ε (Mean = 0.74ε) | Max = 2.23ε (Mean = 0.511ε) |
| medium values | Max = 23.7ε (Mean = 4ε) | Max = 9.00ε (Mean = 0.266ε) |
| large values | Max = 469ε (Mean = 31.5ε) | Max = 20.5ε (Mean = 0.621ε) |
| integer and half-integer | Max = 8.72ε (Mean = 1.48ε) | Max = 4.00ε (Mean = 0.174ε) |
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.
Feedback
Was this page helpful?
Glad to hear it! Please tell us how we can improve.
Sorry to hear that. Please tell us how we can improve.