Accuracy

Accuracy is what this package sells, so it is measured, published and enforced by the test suite, not asserted.

The per-region table is generated: accuracy-table.md. This page explains how it is produced, which algorithm each function uses and where it comes from, what the negative control showed, and where the package is not (yet) as good as it should be.

How the error is measured

Oracle. Every fixture value is computed by mpmath 1.4.1 twice, at 192 and at 320 bits of working precision (with extra precision where the formula cancels, e.g. log Γ(a) + log Γ(b) - log Γ(a+b) for huge arguments, which is evaluated with 1200 extra bits). A point whose two evaluations differ by more than 2^-80 relative is dropped, and the drop is counted in the fixture header: an oracle that silently lost precision must not become the truth the tests are held to. (The first runs dropped 247 logbeta points this way, which led to the extra-precision fix.)

Exact reference, not a rounded one. Each row stores the mpmath value as a double-double hi + lo, hi being the correctly rounded float64. The error of a result y is |y - (hi + lo)| / ulp(hi), so a correctly rounded result scores at most 0.5, and the half-ulp of the reference’s own rounding is not added to anyone’s error.

SciPy on the same points. The generator also stores SciPy 1.18.1’s result for the same float64 arguments, and the Go test measures it exactly as it measures this package. SciPy is a point of comparison, never the reference.

Regions. Each function’s domain is cut into regions — tails, the neighbourhood of zeros and poles, huge and tiny arguments, subnormals, special values (NaN, ±Inf, ±0) — sampled with a fixed seed (log-uniform where the scale matters). About 1000 points per region; 160 000 points over all fixtures. A result that should be NaN or ±Inf and is not (or the reverse) is reported as wrong.

Budgets. TestAccuracy fails if any region’s maximum error exceeds its budget (budgets_test.go: 1.5 times the error measured when the budgets were generated, at least 1 ulp). The suite runs on amd64, arm64, 386, arm, riscv64, loong64, ppc64le and s390x in CI: architectures where Go fuses multiply-adds and ones where it does not must all meet the same budgets.

Regenerate:

uv run --with mpmath==1.4.1 --with scipy==1.18.1 --with numpy==2.5.4 \
    python testdata/gen/gen.py            # fixtures (minutes)
go test -run TestAccuracyTable -accuracy.table=docs/accuracy-table.md
go test -run TestAccuracyBudgets -accuracy.budgets=budgets_test.go

Headline

Over 62 fixtures (419 regions, about 250 000 points), the maximum error of this package is at or below SciPy’s in 418 of 419 regions; the one exception is k1 at small x, 1.05 ulp against SciPy’s 1.05 (both use Boost’s K1 rational; ours is larger in the third digit).

• About half an ulp (correctly rounded or nearly) for the gamma family, erf/erfc/erfcx and their inverses, Ndtr/Ndtri/LogNdtr, zeta and Hurwitz zeta, polygamma, combinatorics, logistic helpers and elliptic integrals — including the tails, the zeros of log Γ at 1 and 2, the poles, the trivial zeros of ζ and the subnormal range, where SciPy’s errors run from tens of ulps to 1e15 ulp or a wrong NaN/Inf.
• Incomplete gamma and its inverses: at most 2.1 ulp (SciPy: up to 8.9e15). Incomplete beta: at most 9 ulp; its inverses at most 20 ulp for parameters >= 1e-3 (SciPy: up to 4.5e15).
• Bessel functions: at most about 5 ulp everywhere, 0.5 ulp at the zeros of J0, J1, Y0, Y1 (SciPy: 1e10–1e12 ulp there) and under 2 ulp for arguments up to 1e300 (SciPy loses every digit).
• Not solved: the beta inverses with a or b far below 1e-3, where the root depends on digits beyond double-double (both this package and SciPy can return a wrong root; the budget is +Inf and documented).

Bold in accuracy-table.md marks a region where ours is at or below SciPy; where both are “wrong” the row is bold too.

Algorithms and sources

Two building blocks run through everything:

• Double-double arithmetic (internal/dd): error-free transformations (Knuth’s TwoSum, Dekker’s Fast2Sum, FMA-based TwoProd) and the double-word algorithms of Joldes, Muller & Popescu (ACM TOMS 44(2), 2017). dd.Exp and dd.Log are table-driven after Tang (ACM TOMS 15(2), 1989 and 16(4), 1990), carried to double-double (about 2^-76 and 2^-72 relative), with 300-bit mpmath tables.
• Fitted polynomials (testdata/gen/fit.py, zz_*_coeffs.go): Chebyshev interpolants of a well-chosen regular part of a function, fitted to mpmath at 200 bits on short pieces, degree chosen so that the interpolation error is below 2^-64, constant term kept in double-double. The function’s zeros, poles and logarithmic singularities are factored out analytically, so the fitted part is smooth and the relative accuracy holds where the function vanishes.

Ported code keeps its attribution in a comment at the function: Cephes (Stephen L. Moshier, BSD via SciPy) and Boost.Math (Boost Software License 1.0).

Negative control

A test suite that cannot fail proves nothing. testdata/negative_control.sh perturbs one coefficient, runs TestAccuracy, and restores the file.

Results (run on this branch, arm64):

The budgets are tight enough that an error of a few ulps in one coefficient of one piece is caught, through every function built on it.

Where SciPy’s result is not a number it should be

Some SciPy “wrong” entries in the table are conventions, not errors; they are listed so the table is not misread:

• gammaln(5e-324) and gammaln(-5e-324) return inf (exact: 744.44).
• erfcinv(5e-324) returns inf (exact: 27.2133).
• betaln(a, b) returns nan for 528 of the 1000 points with a, b in (100, 1e300] (exact values are finite, e.g. betaln(5.08e223, 7.77e195) = -5.05e197).
• Conventions: gammaln(-inf) is -inf in SciPy and +inf here (C99 lgamma); log_factorial(-1) is compared with gammaln(0) = inf where this package returns NaN (n! is undefined for n < 0); polygamma(1, -inf) is inf in SciPy and NaN here.

Known limitations

• Speed. Double-double costs a factor of about 3 single-threaded for most functions and 10–15 for the Bessel functions of real order, Hurwitz zeta and the incomplete beta function (perf.md).
• Digamma at its negative zeros: 31 ulp at 1e-16 relative distance from a zero (SciPy: 2.5e17); limited by double-double itself.
• LogBeta where log B(a, b) ≈ 0 (a, b near 1): 30 ulp; the three log-gammas carry ~2^-56 relative errors that the cancellation exposes (SciPy: 45 000).
• Bessel: Jn/Yn within 5.4 ulp at the zeros of orders >= 2; Jv at negative non-integer order and x < |v| 7.2 ulp (13.2 ulp on loong64, ppc64le and s390x, where Go fuses multiply-adds differently: results are not bit-identical across architectures), Ive at negative order 12 ulp (both from the reflection formula’s cancellation); K1/K1e at small x tie with SciPy (1.05/1.49 ulp, Boost’s rational unchanged); orders above 1e7 only where the series or the large-x expansion applies (NaN elsewhere); J and Y at x = +Inf return the limit 0 where SciPy returns NaN.
• Incomplete beta: up to 5 ulp for a, b in [1, 100) and 9 ulp for a, b in [100, 1e4) near the mean (SciPy: 447 and 3487); BetaIncCInv 20 ulp for a, b in [1e-3, 1) (51 on amd64; SciPy 4.5e15).
• Beta inverses with a or b far below 1e-3: I_x(a, b) is flat to beyond 2^-100 over most of (0, 1), so the float64 root depends on digits past double-double; results can be wrong (so are SciPy’s). Budget +Inf, documented at BetaIncInv.
• Oracle caveat: mpmath’s own betainc returned an identical wrong value at both precisions for b = 1e-200 (the self-check cannot catch an error that does not depend on precision), so the incomplete-beta oracles compute the smaller tail by series or tanh-sinh quadrature instead.
• Complex arguments are not implemented (the nd layer rejects complex dtypes).
• Not yet implemented: Struve, Airy, Kelvin, spherical Bessel, incomplete elliptic integrals, Lambert W.