3 Special Functions

The term “special function” has no formal definition. However, for the purposes of the math library, a special function is one that is not elementary.

Currently, math/special-functions does not export any functions that accept or return complex numbers. Mathematically, some of them could return complex numbers given real numbers, such hurwitz-zeta when given a negative second argument. In these cases, they raise an exn:fail:contract (for an exact argument) or return +nan.0 (for an inexact argument).

Most real functions have more than one type, but they are documented as having only one. The documented type is the most general type, which is used to generate a contract for uses in untyped code. Use :print-type to see all of a function’s types.

Because lambert : Zero -> Zero, Typed Racket proves during typechecking that one of its exact cases is (lambert 0) = 0.

Because the theorem lambert : Flonum -> Flonum is stated as a type and proved by typechecking, Typed Racket’s optimizer can transform the expressions around its use into bare-metal floating-point operations. For example, (+ 2.0 (lambert 3.0)) is transformed into (unsafe-fl+ 2.0 (lambert 3.0)).

The most general type Real -> (U Zero Flonum) is used to generate lambert’s contract when it is used in untyped code. Except for this discussion, this the only type documented for lambert.

3.1 Real Functions

Error is no more than 10 ulps everywhere that has been tested, and is usually no more than 4 ulps.

Error is no more than 11 ulps everywhere that has been tested, and is usually no more than 2 ulps. Error reaches its maximum near negative roots.

Except near negative roots, maximum observed error is 2 ulps, but is usually no more than 1.

Near negative roots, which occur singly between each pair of negative integers, psi0 exhibits catastrophic cancellation from using the reflection formula, meaning that relative error is effectively unbounded. However, maximum observed absolute-error is (* 5 epsilon.0). This is the best we can do for now, because there are currently no reasonably fast algorithms for computing psi0 near negative roots with low relative error.

If you need low relative error near negative roots, use bfpsi0.

From spot checks with m > 0, error appears to be as with psi0: very low except near negative roots. Near negative roots, relative error is apparently unbounded, but absolute error is low.

Mathematically, erfc(x) = 1 - erf(x), but having separate implementations can help maintain accuracy. To compute an expression containing erf, use erf for x near 0.0. For positive x away from 0.0, manipulate (- 1.0 (erfc x)) and its surrounding expressions to avoid the subtraction:

For erf, error is no greater than 2 ulps everywhere that has been tested, and is almost always no greater than 1. For erfc, observed error is no greater than 4 ulps, and is usually no greater than 2.

This function has two real branches. The lambert variant computes the upper branch, and is defined for x >= (- (exp -1)). The lambert- variant computes the lower branch, and is defined for negative x >= (- (exp -1)). The only exact case is (lambert 0) = 0.

The Lambert W function often appears in solutions to equations that contain n log(n), such as those that describe the running time of divide-and-conquer algorithms.

For both branches, error is no more than 2 ulps everywhere tested.

When s is an odd, negative exact integer, (zeta s) computes (bernoulli (- 1 s)), which can be rather slow.

Maximum observed error is 6 ulps, but is usually 3 or less.

When s is an odd, negative exact integer, (eta s) computes (bernoulli (- 1 s)), which can be rather slow.

Maximum observed error is 11 ulps, but is usually 4 or less.

While hurwitz-zeta currently raises an exception for s < 1, it may in the future return real values.

Maximum observed error is 6 ulps, but is usually 2 or less.

If you are doing statistical work, you should probably use gamma-dist instead, which is defined in terms of gamma-inc and is more flexible (e.g. it allows negative x).

If you are doing statistical work, you should probably use beta-dist instead, which is defined in terms of beta-inc and is more flexible (e.g. it allows negative x).

Similar identities should hold as with gamma-inc.

While most areas of this function have error less than 5e-15, when a and b have very dissimilar magnitudes (e.g. 1e-16 and 1e+16), it exhibits catastrophic cancellation. We are working on it.

3.2 Flonum Functions