Precision guide

Most Nautilus numerical functions operate in f32 (IEEE 754 single precision), providing approximately 6-7 significant decimal digits. Nautilus.Special functions accept f32 or f64 inputs. Nautilus.Rolling uses f64 series.

The table reports f32 errors. Improvements from f64 vary by function. Six functions (gamma, log_gamma, beta, lbeta, ellipk, and ellipe) are limited by f32 rounding rather than by their own coefficients, so at f64 their errors are much smaller. The other approximations improve by less, depending on the function and argument. Measure the argument range your calculation uses.

Treat the table as an f32 guide. It has no figure for beta or lbeta, and bessel_y1 has an absolute error of about 1.2e-4 at either dtype just below its large-x seam at x = 7.5, well away from any zero.

The measured dtype differences include:

• erf is capped by Abramowitz & Stegun 7.1.26's own 1.5e-7 error. At x = 0.5 an f64 erf measures 1.385e-7 against an f32 erf's 1.861e-7; over the whole range the maxima are 1.3884e-7 and 4.438e-7.
• airy_ai and airy_bi above x = 5 use only the leading asymptotic term and gain nothing from f64 there; the two dtypes agree to three digits. Below x = 5 f64 is far better.
• f32 erfc returns exactly 0 from about x = 3.92 onward; f64 returns 1.546e-8 at x = 4 and remains nonzero to about x = 5.5. Its relative error is still 2.8e-3 at x = 4: subtracting erf(x) from 1 cancels nearly equal values while retaining erf's absolute approximation error. A nonzero f64 tail is not necessarily accurate to three significant digits.

Precision by function family

Family	Error guide (relative unless marked absolute)	Notes
erf	~1e-7	Horner rational approximation
erfc	same absolute error as erf	1 - erf(x), so relative error grows in the upper tail
erfinv	~1e-8	Acklam rational approximation via norminv
gamma	~1e-7	Lanczos (g=7) with reflection
log_gamma	~1e-9	Lanczos (g=7) with reflection
digamma	~1e-7	Recurrence + asymptotic (x >= 6)
trigamma	~1e-6	Recurrence + asymptotic (x >= 6)
ellipk, ellipe	~1e-8	AGM recurrence (quadratic convergence)
bessel_j0, j1	Use absolute comparisons near zeros	Rational polynomial + large-x trig
bessel_y0	Use absolute comparisons near zeros	Rational + log-singularity
bessel_y1	~1.2e-4 absolute at x ≈ 7.4	Large-x branch from x = 7.5; use absolute comparisons near zeros
bessel_i0, i1	f32	Polynomial + asymptotic, crossover 3.75
bessel_k0, k1	f32	Polynomial/log + asymptotic, crossover 2.0
airy_ai, airy_bi	f32 for |x| <= 5	See large-negative-x note below
Distribution CDFs	~1e-5 to 1e-7	Depends on underlying special functions
normal_cdf	~1e-7	Delegates to erf
normal_inv_cdf	~1e-7	Acklam rational via erfinv

Known precision issues

Bessel function zeros

For bessel_j0, bessel_j1, bessel_y0, and bessel_y1, a small absolute error can produce a large relative error near a zero: relative error divides the absolute difference by the reference value's magnitude. At an exact zero, that ratio is undefined. Compare absolute differences there and calibrate tolerances against a reference over the arguments your calculation uses. The bessel_y1 measurement at x ≈ 7.4 describes its approximation-branch boundary, not an error bound near zeros or across the full domain.

Airy functions at large negative x

airy_ai uses a power series for |x| <= 5 and an exponential asymptotic for x > 5. For large negative x, only the power series is available (no asymptotic oscillatory branch is implemented). The function enters an oscillatory regime for x < 0 and the power series degrades for |x| much larger than 5. airy_bi has the same limitation but is less affected because it grows exponentially for positive x where the asymptotic branch covers it.

Cancellation in subtraction-heavy expressions

Any computation involving subtraction of nearly-equal f32 values will suffer catastrophic cancellation. This affects:

• cosine_distance when vectors are nearly parallel (1 - sim near 0)
• variance_vec for data with very small variance relative to the mean
• gamma_cdf for extreme shape/scale ratios

When f32 precision is insufficient in Nautilus.Special, call it at f64 directly. Its functions have an explicit {f32, f64} dtype bound. Read the caveats above first, because several of them are coefficient-limited rather than dtype-limited. The distribution, linear-algebra, and solver APIs described here are f32-only.

Nautilus.Special does not admit f16 or bf16; these calls are type errors.

Comparison to scipy

SciPy typically computes in f64, while most Nautilus functions return f32. Compare values at the precision and inputs appropriate to the calculation; the functions' approximations have different numerical domains.