Deviations from the reference

A parity fixture normally holds the value NumPy or SciPy returned, and Tessero must match it at the default tolerance (relative 1e-13 for element-wise functions, 1e-12 for calls; random draws exactly). This page lists every exception and the rule behind it. Each rule is applied by the fixture generators (tools/parity/fixtures*.py), never by editing a fixture, and the counts below are generated from the fixture files by tools/parity/report.py.

The reference environment#

The fixtures are generated by tools/parity/gen-fixtures.sh, and the generators refuse to write outside it. The script disables NumPy's AVX-512 dispatch targets for float64 loops (NPY_DISABLE_CPU_FEATURES="X86_V4 AVX512_ICL AVX512_SPR") and fixes OpenBLAS's kernel and thread count. On AVX-512 machines NumPy evaluates exp, log, power, expm1, log1p, sinh and arcsin with its own vector kernels, which differ from the C library by an ulp on some inputs. The reference would then depend on the machine that generated it. With those targets disabled, NumPy's results equal the C library's; 200 000 inputs per function showed no difference. Each fixture file records this environment and NumPy's per-function dispatch under "env" ([ADR 0012(https://github.com/DMJ-CV-91913/tessero/blob/main/docs/architecture/adr/0012-pinned-reference-environment.md)).

Rules#

Guards: SciPy crashes or never returns. On some inputs SciPy itself segfaults, aborts or loops forever. Examples are the Mathieu functions at infinite arguments, the spheroidal wave functions for large parameters, and cdflib's non-central t sum. tools/parity/crash_scan.py evaluates each row in a forked process and records those rows in tools/parity/crashes-special.json. The function's spec entry carries a guard: (or capped:), Tessero returns NaN there, and the fixture expects NaN.

Undefined rows: SciPy's value is undefined. Some rows read memory SciPy's code never wrote. There are two ways to find them. tools/parity/ub_scan.py runs SciPy's original code under AddressSanitizer. The generator also re-evaluates every reference_ub function in fresh interpreters: once in reverse row order (a different heap history), and three times under glibc's MALLOC_PERTURB_, which fills every new allocation with a byte pattern (two garbage patterns and zero). A row whose value, or the sign of whose zero, changes in any of them depends on uninitialised memory. Both kinds of row are recorded as "undef" and not compared. Tessero's copies of those work arrays are zero-initialised (THIRD_PARTY_NOTICES.md), so its own value on such a row is deterministic; the full fixture suite passes with MALLOC_PERTURB_ set, which tools/parity/gate.sh checks.

U1: invalid shape parameters. For parameters that fail _argcheck, SciPy's public methods return nan (its badvalue). A few private methods (entropy, moment) raise instead when evaluated on the empty set of valid parameters, for example 0.0 ** negative, an index into an empty array, or BiasedUrn's own range check. The fixture expects nan and keeps SciPy's exception under "upstream".

U2: an array result that differs from the element-by-element results. A public method applied to an array must equal the method applied to each element. norminvgauss.sf and isf truncate their output through zip() on arrays. The fixture holds the element-by-element values and keeps SciPy's array under "upstream".

Per-call tolerances derived from the data. epps_singleton_2samp computes g′ pinv(C) g, where C is an estimated covariance matrix. A last-bit difference in C (Tessero's SVD and LAPACK's gesdd sum in different orders) is amplified by cond(C). The generator computes cond(C) from each call's own data and records rtol = max(1e-12, cond(C)·ε), the first-order perturbation bound of the pseudo-inverse, with the reason.

NumPy computes inner, vdot, tensordot and the convolution behind polymul on floats with BLAS dot products, which sum in the BLAS kernel's order; Tessero sums in index order. Both lie within γₙ·Σ|aₖbₖ| of the exact dot product (γₙ = nε/(1 − nε)), so the generator records rtol = 2γₙ·Σ|aₖbₖ|/|result| for a call where that bound exceeds the default, with the reason (tools/parity/fixtures_np_numeric.py).

Resource limits: Tessero raises where SciPy would try. The failure is always an exception, never an abort, a crash or a hang in a PHP worker. These cases are covered by tests/Unit/AdversarialInputTest.php on both backends, under AddressSanitizer as well:

• Arrays sized by a parameter (the support of a discrete distribution, np.arange(n + 1), histogram bins, L-moment orders): Tessero raises MemoryError, or "Maximum allowed size exceeded", when the size exceeds the memory budget (or the machine's RAM when no budget is set). NumPy may instead attempt the allocation.
• A moment order beyond 2^31: Tessero raises. SciPy's private _munp takes a Python int there.
• An infinite moment order in stats.moment: Tessero raises. SciPy never returns.
• Polynomial derivative or integral orders of 1000 or more (np.polyder, np.polyint): Tessero raises. NumPy recurses once per order and raises RecursionError there.

NumPy's own undefined behaviour: Tessero raises or returns a defined value. No fixture pins these, since NumPy's result is not reproducible; tests/Unit/AdversarialInputTest.php covers Tessero's side.

• np.bincount with a value of 2^63 − 1 returns an empty array in NumPy 2.4 and corrupts its heap (the length max + 1 overflows). Tessero raises.
• np.trace and np.diagonal with offset = -2^31 segfault NumPy. NumPy raises OverflowError for offsets outside a C int, and Tessero raises for those as well, but returns the defined (empty) diagonal at −2^31.
• np.tri, np.tril, np.triu and the *_indices functions with a float k whose magnitude is 2^63 or more and a non-empty result: NumPy casts -k to int64, which is undefined. Tessero raises, with the message of NumPy's OverflowError. Below that magnitude, Tessero follows NumPy's arithmetic exactly, including the wrap-around of an integer k near −2^63 and the empty result for k = ±1e20.

Uninitialised NumPy output. numpy.empty_like and numpy.pad(mode='empty') leave values unspecified: those calls compare only the shape and dtype (compare mode shape). Tessero fills them with zeros.

Integer inputs. Where NumPy keeps an integer dtype (ptp, the nan-reductions, cumulative sums and products, partition, correlate/convolve, unique and the set operations, trace, and quantile's discrete methods), Tessero computes in int64 with NumPy's wrap-around, so results are exact above 2^53 as well, and the fixtures compare integers exactly.

Python values and PHP arrays. A PHP list becomes one array before the kernel sees it, so a list that mixes ints and floats is float64 as a whole. NumPy's put, putmask, place, fill_diagonal and put_along_axis convert such a list element by element, so an int above 2^53 in a mixed list keeps its exact value in NumPy and is rounded in Tessero. Pass an int64 NDArray (or a list of ints only) to keep it exact.

NumPy's integer arithmetic in quantiles. For integer samples NumPy's interpolating quantile methods subtract neighbouring order statistics in the sample's integer type, which wraps around when they differ by more than the type's range (int32 samples about 2^31 apart), and return a value outside the sample's range. Tessero interpolates in float64 and returns the value NumPy returns for the same sample as float64. The discrete methods return a sample and agree exactly.

Empty or colliding float inputs, where NumPy's result is incidental. nanquantile of an empty slice with an array q and an axis returns NumPy's nanmean shape, without the q axis. Tessero keeps the q axis, so the shape is the one NumPy gives for a non-empty input. intersect1d and setxor1d of an int64 and a float64 array compare in float64 in both. After unique, NumPy can still treat two int64 values above 2^53 that round to the same float as a match, but Tessero merges them first. unique_values returns 0.0 where NumPy's hash order keeps -0.0.

Known, not guarded. Where SciPy's own algorithm is slow for huge arguments, Tessero runs the same algorithm: orthogonal polynomials with n around 10^9 loop n times, and anderson_ksamp, siegelslopes and percentileofscore are quadratic in the sample size. Callers who pass untrusted sizes should bound them. See the open issues.

Known: undefined behaviour inside SciPy's own code. On ±inf or very large arguments, the vendored special function code (xsf's specfun and cephes, Boost.Math) and NumPy's samplers perform signed integer overflow that UBSan reports: for example lpmv, yn, hyperu, pbdv and smirnovi near 2^31, and the binomial sampler near 2^63. This is the same code SciPy runs, and the sanitizer sweep finds no memory errors on those paths. Tessero leaves them as SciPy has them, so the results stay identical, and records them here; a fix belongs upstream. One small divergence remains: SciPy's betabinom.rvs raises a dtype-cast TypeError for n near 2^63, while Tessero draws.

Current counts#

Edit this page on GitHub