Evaluating special functions

FLINT's special functions take balls as input (arb_t, acb_t, or power series with ball coefficients) and return balls that contain the exact value of the function at every point of the input. Every algorithm therefore comes with an error bound: for a truncated series, for the propagation of the input radius, and for the rounding errors, which ball arithmetic tracks automatically. Most of the code goes back to Arb [Joh2017; Joh2014a].

Principles#

Enclosures, not approximations. A result [m \pm r] is correct if it contains the true value; it is useful if r is small. The functions are written so that r is close to the best possible for exact inputs at the target precision, but they never return an incorrect enclosure to obtain a smaller one. When no finite enclosure can be given — a pole in the input, an input that straddles a branch cut, a missing algorithm, or too low a precision for the size of the arguments — the result is indeterminate (infinite radius). The semantics and the cost model of ball arithmetic are explained in Ball arithmetic and Ball arithmetic.

Truncation bounds. Every series is truncated with an explicit bound for the remainder: tail bounds by geometric series for convergent hypergeometric series (Convergent series), Olver-type bounds for asymptotic expansions (Asymptotic expansions, The Stirling series), Cauchy majorants for Taylor series of solutions of differential equations (The corner case), and Euler–Maclaurin remainder bounds for zeta functions (Zeta functions, polylogarithms and L-functions). The bounds are computed with the mag_t type, which rounds upward, so that bounding costs little compared to the evaluation.

Working precision. Ball arithmetic propagates rounding errors correctly but does not prevent cancellation: summing \sum (-x)^k/k! for x = 100 at 64 bits gives a correct but useless ball. The functions therefore estimate the cancellation of the method they use, usually in double precision, and add that many guard bits; conversely they lower the working precision to the accuracy of the input, since more precision cannot improve the result.

Midpoint evaluation and propagated error. For an input ball x = [m \pm r], evaluating an algorithm directly in ball arithmetic often overestimates the result (each operation widens the ball, and correlations are lost). Many functions instead evaluate at the exact midpoint m and add r \cdot \sup_{|t - m| \le r} |f'(t)|, using a cheap bound for the derivative: e^{-z^2} for the error function, \psi(x) for the gamma function, Bessel and Airy function bounds, Legendre polynomial derivative bounds. This also allows the evaluation at the midpoint to use a lower precision.

Choosing among expansions. Most functions have several representations — a convergent series at the origin, an asymptotic expansion at infinity, transformations of the argument, integral representations — with very different costs depending on the argument, the parameters and the precision. FLINT chooses among them with cheap floating-point estimates of the number of terms and of the cancellation, and some functions retry with another method if the result is inaccurate.

Argument reduction and functional equations. Recurrences and symmetries move the argument to a region where an expansion converges quickly: \Gamma(z) = \Gamma(z+r)/(z)_r and the reflection formula, the linear fractional transformations of {}_2F_1 and the dilogarithm, the modular group for theta functions and modular forms, \exp(x) = \exp(x/2^r)^{2^r} and table lookups for the elementary functions.

Summation algorithms. The cost of a series is dominated by the multiplications involving full-precision numbers. Rectangular splitting reduces their number to O(\sqrt N) for N terms when the term ratios have small coefficients; binary splitting computes series with short arguments and rational coefficients in quasi-linear time; the bit-burst algorithm extends this to arbitrary arguments for holonomic functions (Summation algorithms).

Caching. Constants (\pi, \log 2, Euler's \gamma, \Gamma(1/3), ...) and Bernoulli numbers are computed once at the highest precision requested so far and cached per thread or globally (Mathematical constants).

Power series. Functions of power series arguments (_series) and Taylor expansions at a point (_jet) give derivatives with rigorous bounds. They are computed either from the differential equation of the function (Bessel, Airy, Coulomb functions), from the derivative and power series integration (error functions, exponential integrals), or by running the scalar algorithm with power series arithmetic (gamma function, hypergeometric series with power series parameters, which is also how limits at removable singularities are computed).

Map of functions and methods#

The power series versions in arb_poly and acb_poly (arb_poly_gamma_series, acb_poly_zeta_series, ...) use the same algorithms with power series arithmetic.

Interfaces#

The functions are available at three levels.

• Ball arithmetic. The modules in the table, with arb_t and acb_t arguments, are the primary interface. The real functions of arb_hypgeom mostly wrap the complex ones and return an indeterminate result when the value is not real.
• Generic rings. gr_special defines gr_gamma, gr_erf, gr_hypgeom_2f1, gr_bessel_j, gr_zeta and many more for every ring; the contexts of real and complex balls call the functions above, the floating-point contexts of nfloat and mp_real have their own implementations of the elementary functions, and exact contexts such as ca return exact values where they can (for example \Gamma(1/2) = \sqrt\pi). Some power series algorithms exist only in generic form, such as _gr_poly_bessel_j_series.
• Floating-point wrappers. arb_fpwrap evaluates the functions for double and complex double arguments with correct rounding or a guaranteed relative accuracy, increasing the precision until the ball determines the result.

Use in FLINT#

The module pages give the conventions (precision, branch cuts, regularization, optional outputs); the guide Using acb_hypgeom explains them for the hypergeometric functions with examples.

@misc{flint-doc-algorithms-special-functions-overview,
  author       = {{The FLINT team}},
  title        = {{Evaluating special functions}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/special-functions/overview.html},
  note         = {Generated 2026-10-20}
}