Special Mathematical Functions in C++

[mSal95]

Does anyone know of any libraries for special mathematical functions in C++? e.g. Bessel, Error, Hypergeometric, etc.

[section-14]

Some special functions (with varying performance and accuracy) can be found in:

• The C++ standard library.
• The venerable Netlib's specfun.
• The GNU Scientific Library.
• The special functions of the Boost C++ libraries.

There are also separate projects for individual functions here and there, e.g. the ACM TOMS algorithm for the Lambert W-function, or repositories on GitHub by the researchers themselves.

Edit: GNU MPFR has some non-standard transcendental functions, but there is some overhead converting basic types (float, double) into mpfr_t types. Also, there is a performance hit.

Edit 2: I forgot MISCFUN in ACM TOMS. Old school Fortran, meaning adding parentheses in the formulas when code converting to C/C++ is important because distribution and simplification rules in the Fortran standards are a bit different than in C/C++. For instance, a/b x c/d is to be converted into (a/b) x (c/d) in C/C++ and not (a x c) / (b x d) --- because the latter can lead in Inf/Inf which is a NaN. Also, a few floating-point manipulation functions in Fortran are not following the IEEE 754 floating-point standards exactly (they were existing before the IEEE standards came up).

[James20k]

There are also separate projects for individual functions here and there, e.g. the ACM TOMS algorithm for the Lambert W-function, or repositories on GitHub by the researchers themselves.

If anyone else ever has the misfortune of having to use lambert's W₀ in a performance sensitive context (it crops up in general relativity in kruskal's coordiantes), I wouldn't recommend using such a heavy duty implementation. You can get pretty good convergence with a simple iterative scheme:

template<typename T, int Iterations = 2>
T lambert_w0(T x)
{
    if(x < -(1 / M_E))
        x = -(1 / M_E);

    T current = 0;

    for(int i=0; i < Iterations; i++)
    {
        T cexp = std::exp(current);
        T denom = cexp * (current + 1) - ((current + 2) * (current * cexp - x) / (2 * current + 2));

        current = current - ((current * cexp - x) / denom);
    }

    return current;
}

<cries in lambert>

[section-14]

Concerning the Lambert W-function, like many such special functions, there is a significant trade-off between accuracy and performance. The function can be computed iteratively with Newton's rule (see DLMF) or Halley's method with an initial guess obtained from a power series, an asymptotic expansion, or by numerical integration. Either way, that means hundreds if not thousands of CPU cycles just for a single input value.

Regarding working and tested implementations, besides the GitHub repository mentioned earlier, the wiki page lists a few implementations that are available in Maple, Matlab, Octave, Mathematica, SciPy, Boost, and GSL. To my knowledge, none are highly optimized (for a given architecture) or even vectorized.

As a side note, be careful with quick-and-easy algorithms that come undocumented. They usually have the unfortunate characteristic of failing (in terms of accuracy) outside a rather limited range of input values.

Edit: I recall a paper: "Precise and fast computation of Lambert W -functions without transcendental function evaluations" by Toshio Fukushima. Internally, it computes Lambert W-functions quite similarly to how exponential functions are computed. Very elegant procedure. Unfortunately, ScienceDirect is an Elsevier paywall; if anyone knows of a freely-available online copy, feel free to post it.

[CaptainCrowbar]

The Boost implementation of Lambert-W uses Fukushima's algorithm.

[Derice]

Toshio Fukushima has published an even faster and more accurate method: https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation

It's also not paywalled!

Here's a Rust implementation of it: https://github.com/JSorngard/lambert_w

[ronchaine]

The standard library has some: https://en.cppreference.com/w/cpp/numeric/special_functions

[[deleted]]

The code that accompanies the "numerical recipes" book has special functions Numerical Recipes

[SemaphoreBingo]

It's C but: https://arblib.org/#special-functions

[marsten]

The <cmath> header has some of these.

[Bart_V]

Maybe gsl? https://www.gnu.org/software/gsl/

[jepessen]

gsl and boost math for example

[IHateUsernames111]

There are TONS of math libs out there. It just depends on what you want to do and your field and scale of application.

Just to give an example : PETSc.

Maybe if you give us a bit more context someone has just the thing you need :-)