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 incomplete-gamma Boost Error ratesThe 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-gammaSee 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 InfinityThe 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 dataThe 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.
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4