Uni20
Asynchronous tensor and runtime library.
Loading...
Searching...
No Matches
Numerical precision validation

Numerical support is a claim about an operation, scalar precision, and execution backend. Compiling a scalar or passing tests for another backend is insufficient. Use the shared cases in tests/numerics when extending scalar-generic algorithms. Existing algorithm-specific regression suites remain necessary for layouts, status reporting, restarts, and exceptional inputs.

Shared cases

precision_cases.hpp describes real and complex cases for float32, float64, float80, float128, and mpreal at 128 and 256 bits. Both MPFR cases use the same C++ scalar type. Values, tensor defaults, epsilon queries, and error comparisons are constructed in that case's precision. MPFR exact arithmetic has separate tests: passing an exact rational identity does not demonstrate finite-precision behavior.

The shared contexts construct scalars with scalar_at_precision<R>(n, p) and tensors with a trailing p for fixed and runtime types alike. The GenericPrecisionQueries and GenericPrecisionOperations probes cover this common API, including fixed-type mismatch rejection before output mutation. The operation probe uses default dispatch; separate forced-provider probes remain responsible for provider certification.

All twelve cases remain in the coverage matrix when optional dependencies are disabled. precision_registry.hpp declares which operation/scalar/backend combinations are configured and supported. Only those combinations become executable Google Tests. Unsupported, unavailable, and mathematically inapplicable combinations appear in the coverage report, without creating skipped tests or counting as passes.

The expectations are independent of production concepts: accidentally removing a supported implementation must fail compilation or its numerical test instead of shortening the tested type list. NumericalCoverage.RegisteredProbesMatchDeclaredMatrix also checks that every expected probe is registered exactly once and exports the full declared matrix into the test XML. Update the registry and operation-specific expectations when adding support. A fully configured build should run all its registered tests; reserve runtime skips for conditions discovered during execution, such as a multi-NUMA-node test on a machine with only one visible NUMA node.

Float128 means a 113-bit significand, not MPLAPACK's fallback alias to a narrower long double. That alias mode is reported as unavailable for this case. Float80 similarly requires Uni20's native 64-bit-significand format detection.

Dense decomposition probes force LapackBackend for real/complex QR, LQ, self-adjoint eigensystems, SVD and general eigenanalysis. They cover float32/64, configured binary80/binary128 and MPFR at 128/256 bits. Binary80 complex decomposition cells require a native long-double provider. Reconstruction and orthogonality use scalar loops independent of provider GEMM. Spectral gaps detect narrowing, and a separate 512-bit analytic oracle checks improved QR/eigenvalue/SVD accuracy.

Evidence required for numerical support

The least-squares probes force CPU and LAPACK independently, covering inconsistent rectangular fits, complex residual orthogonality, exact-rational solution references, retained low bits, extended exponent range and a 512-bit accuracy oracle when MPFR is configured. CPU probes include finite MPFR/MPC arithmetic without requiring the MPLAPACK provider.

ProviderSvdLeastSquares and ProviderLuCondition independently force LapackBackend across the configured real/complex scalar cases. They test minimum-norm solutions, epsilon-based rank decisions, original norms, conditioning, retained precision and extended exponents, with independent 512-bit rational oracles for accuracy. Unsupported binary80 complex providers and disabled MPLAPACK cells remain visible in the coverage report.

  1. Correctness: run the same deterministic mathematical problem across the case list, using exact analytic data or an independent higher-precision oracle. Include nonzero imaginary components in complex tests.
  2. Retained precision: use inputs and outputs whose distinguishing increments disappear at a lower precision. Tests must reject narrowing in input loading, workspaces, tolerances, kernels, and result construction where those paths are exercised.
  3. Improved accuracy: use a controlled problem where additional precision lowers a measured error. Do not demand monotonic improvement on arbitrary iterative problems or compare two implementations sharing the same oracle.
  4. Exponent range: test this separately from significand precision. A value outside double's range says nothing about retention of float80's low bits.
  5. Runtime precision: run the same algorithm at multiple MPFR precisions and inspect the output precision. Exact and unspecified states need their own boundary tests; they cannot supply an implicit epsilon.
  6. Execution coverage: force the backend when testing a provider. A successful fallback does not certify that provider. Async, device, view, and distributed paths require additional tests as they gain support.

Avoid EXPECT_NEAR for scalar-generic numerical assertions: its arguments are double. The shared expect_error_at_most computes the error in the tested real type and formats only diagnostics. Construct decimals from text or exact integer ratios; casting a double literal does not create higher-precision data.

Choosing a numerical assertion

Use exact equality when the contract demands it, including provider-wiring checks against the same correctly rounded operation. Use an explicit absolute/relative bound for residuals, accumulated errors and identities near zero. EXPECT_FLOATING_EQ and ASSERT_FLOATING_EQ from <uni20/common/gtest.hpp> instead count representable steps between results; their default tolerance is four ULPs. They support native float32/64/80/128, native complex values componentwise, and mpreal when MPFR is enabled. complex<mpreal> ULP assertions are not yet supported.

For mpreal, the comparison contract is:

  • Two approximate operands must have the same precision, even if their values are equal. A mismatch fails with a diagnostic; convert a reference explicitly with .at(p) when that is intended.
  • An exact operand is rounded to the approximate operand's precision, nearest with ties to even, before counting steps. This includes MPFR's ordinary underflow/overflow behavior. Two exact operands require exact equality; no tolerance gives unequal rationals a floating-point grid.
  • The grid uses that precision and MPFR's current exponent range, without changing either. MPFR has no subnormals: zero and the smallest positive representable value are one step apart. Both signed zeros occupy the same position. Finite operands outside the current exponent range fail explicitly.
  • NaNs and unset values fail. Equal infinities pass at matching precision; other comparisons involving infinity fail. Negative tolerances fail.
  • Diagnostic distances saturate at max<long long> in magnitude. The actual comparison uses the full integer distance, so saturation cannot turn a failure into a pass. Exact inequality and incomparable pairs have no distance and use the same sentinel, with an explanatory assertion message.
auto actual = uni20::mpreal(1, p) / 3;
auto reference = uni20::mpreal("1/3", uni20::Precision::bits(512));
EXPECT_EQ(actual.precision(), p);
EXPECT_FLOATING_EQ(actual, reference.at(p), 1);
static constexpr Precision bits(bit_count_type count)
Select a positive binary significand size.
Definition core/precision.hpp:21
static constexpr Precision exact() noexcept
No approximation has been selected; rational arithmetic remains exact.
Definition core/precision.hpp:17
Owning exact rational or finite-precision MPFR real, with eager arithmetic.
Definition mpreal.hpp:74
#define EXPECT_FLOATING_EQ(a, b,...)
Definition gtest.hpp:88
Definition data_table.cpp:8

Operand expressions and the tolerance are evaluated once. Failure messages show round-trip MPFR values, exact/unset states, working precisions and the reason an invalid pair cannot be compared. Comparison preserves the operands and MPFR precision defaults, exponent range and exception flags. The shared check::FloatingULP<mpreal> comparator also serves CHECK_FLOATING_EQ and PRECONDITION_FLOATING_EQ; the detailed ULP diagnostic is supplied by the GoogleTest assertions. ULP proximity does not check the requested output precision or establish the independence of a reference, so keep those checks separate.

First operation families

Scalar math dispatch uses the same real/complex precision cases. A near-one square-root problem and a magnitude retaining an increment beyond the preceding precision detect narrowing in the callable interface. Core tests separately exercise ADL lookup, unsupported types and exact/explicit-precision semantics; these mechanism tests do not replace the numerical precision cases.

The uni20_numerical_precision_tests target currently covers:

Family Evidence
Scalar epsilon and arithmetic Runtime epsilon, retained increment, complex reciprocal accuracy
Real elementary math Small-argument expm1/log1p with cancellation negative controls; inverse identities retaining increments; error improvement against independent 512-bit series
Real numerical utilities Named decompositions and remainders retaining increments; fused cancellation; rounding ties, signed zeros, NaN/infinity and grid endpoints; integer power/root accuracy against rational and independent Newton references
Native/MPFR Gamma and error functions Accuracy improvement against sqrt(pi) and an independent 512-bit error-function series
Native/MPFR Beta providers Rational identities without MPFR; increased accuracy against an independent pi identity
Typed constants Type retention and improved accuracy of pi, log(2), Euler gamma and Catalan, including binary128 constants and an independent Catalan series
MPFR special functions Gamma/digamma/Beta recurrences, upper-Gamma identity, Bessel Wronskian, zeta/pi identity; accuracy against Ei, Li2, Bessel J and Airy series and an AGM iteration
Native/MPFR pi-scaled functions Exact phase at large arguments; zeros, poles, domains and subnormals; four-grid-spacing checks against 512-bit MPFR; precision improvement from analytic identities and an atan series; direct portable-path validation
Native/MPFR upper incomplete Gamma Integer/half-integer identities, negative shapes, shape-zero limits, transition joins, huge-parameter cancellation and final range boundaries; wider references with explicit cancellation control and precision gains
Native/MPFR Airy Ai Differential equation and value-at-zero identities, zero and region-join neighbors, full-range negative phase, positive subnormal/zero boundaries; cross-checked asymptotic references and precision gains
Native zeta / MPFR Exact trivial zeros, pole neighbors, Gamma/Stirling joins, negative overflow boundaries, rational Bernoulli and pi identities, original-input error bounds and precision gains
Native Ei / MPFR Public probes at all real precisions, cancellation near the positive zero, continued-fraction and asymptotic joins, normal/subnormal/zero/overflow boundaries, independent small-argument series and upper-Gamma identity
Native/MPFR real dilogarithm Whole real line, subnormals, continuation for x>1, nonzero-root neighbors, independent pi/log values and defining series; original-input MPFR512 references with 1024-bit worst-case checks
Native/MPFR digamma Absolute error near zeros, relative accuracy away from zeros, reflection/pole neighborhoods, full exponent range and independent harmonic identities
Native Beta fallback / MPFR Separate portable and public probes; Gamma-ratio cancellation, extreme parameter ratios, original-input error bounds, range limits and independent integer products; large-argument log-Gamma oracle with additional cancellation guard bits
Typed factorial Every native finite order, exact-integer references, overflow boundaries and checked result precision
Native/MPC complex theta jets Value and first two derivatives; conjugation and quasi-periods; independent differentiated Jacobi products at 1280 bits and precision gains
Native/MPFR theta constants Scaled rectangular constants, Jacobi products and a Gamma-function reference; modular joins, truncation budgets, full native range and precision gains
Native/MPFR AGM Extreme operand ratios, subnormal inputs, symmetry, equal-input identities and independent elliptic-integral series
Native/MPFR reciprocal functions Original-input phase, signed zeros, poles and hyperbolic tails beyond sinh/cosh overflow; full-range oracle scans and independent series references
Native/MPFR compound Original increment and integer exponent retention; signed64 limits; range and subnormal neighborhoods; independent exact rational references; public and portable paths
Native/MPFR exponential extensions Small-argument cancellation; logarithmic domain and range boundaries; public and direct portable probes; recorded maximum errors/worst inputs; independent radical and log-series accuracy references
MPFR elementary extensions Small-argument cancellation controls, large pi-scaled argument reduction, inverse and reciprocal identities, compound rounding without loss in 1+x
CPU matrix one-norm and tensor reductions Retained increment, conjugating inner product, squared-norm identity
CPU and provider GEMM Analytic real/complex products with an increment near each precision's resolution
CPU and provider square solve Nearly singular dyadic system with analytic solution, independent scalar residual
Solve accuracy Well-conditioned system with exact rational solution, compared with a 512-bit oracle
CPU and provider reusable LU Small-gap system, multiple RHS columns and repeated solves; rational solution against a 512-bit oracle
CPU/provider float32-factor refinement All configured real/complex working scalars, improved rational-solution accuracy, MPFR residuals at twice working precision; does not imply provider coverage for the working scalar
CPU and provider log determinant Near-one magnitude with an independent analytic logarithm bound; nontrivial complex phase; extended exponent range
Extended-range solve System scaled by 2^-1200, independent of the small-gap test
Projected tridiagonal eigensystem Analytic closely spaced eigenvalues
Lanczos Closely spaced eigenvalues and analytic componentwise eigenvector residuals
Arnoldi Closely spaced real and complex eigenvalues
Hermitian exponential action Diagonal analytic exponential, retaining small output increments

The shared gap is 2^(8-p) for a p-bit significand. Tolerances scale with the case's epsilon and are smaller than the gap. Narrowing controls verify that the gap disappears in a lower native format for the higher-precision cases. Log-determinant range probes additionally factor a non-diagonal matrix scaled by 2^±1200 and 2^±12000: its exact LU factors remain representable, but unscaled complex division can overflow or underflow internally. This exercises elimination as well as the logarithmic reduction.

For the non-dyadic reciprocal and solve fixtures, the oracle transfers each computed binary value exactly to 512-bit MPFR. It compares the error with both the case's epsilon and the error from rounding the analytic solution at the preceding precision (12, 24, 53, 64, 113, or 128 bits as appropriate). The result must be at least sixteen times more accurate than that lower-precision value. The 12-bit case is a comparison baseline, not another supported Uni20 scalar. Without MPFR, these oracle-dependent checks are explicitly unavailable; the analytic increment, range, and residual checks still run.

This matrix generalizes the former binary128-only one-norm, reduction, extended-range solve, tridiagonal, Lanczos, complex Arnoldi, and Hermitian exponential probes. Specialized binary128 Schur, triangular Arnoldi, matrix exponential, SVD, and QR/LQ regressions remain in their existing suites.

Running and reporting

Use the normal project configuration, adding MPFR/MPC and MPLAPACK when testing those providers. See testing, MPFR scalars, and MPLAPACK setup for configuration details.

cmake --build build_codex/<configuration> --target uni20_numerical_precision_tests
ctest --test-dir build_codex/<configuration> --output-on-failure -R '^Numerical'

For an artifact that records actual execution results, run the executable with Google Test XML output, then render that XML:

build_codex/<configuration>/tests/numerics/uni20_numerical_precision_tests \
--gtest_output=xml:build_codex/<configuration>/precision.xml
python3 scripts/report-numerical-precision.py \
build_codex/<configuration>/precision.xml \
--output build_codex/<configuration>/precision.md

The report joins the declared matrix to actual results and distinguishes passed, failed, unsupported, unavailable, and not_applicable. A missing, unexecuted, or skipped expected probe is an error, as is an undeclared probe. Use an unfiltered execution: a filtered XML file cannot stand in for complete coverage. Missing optional dependencies never count as passes. The coverage report retains the unsupported combinations even in a build with no skipped tests. CI runs the matrix in ordinary configurations and saves XML/Markdown artifacts in the MPLAPACK job.

Certification boundaries and remaining work

This is the first shared test matrix, not complete certification of every numerical subsystem. In the current implementation:

  • Float80 CPU arithmetic, reductions, GEMM, solve, reusable LU and log determinants have shared probes. The optional binary80 provider runs the same GEMM, solve, LU/log-determinant and dense-decomposition probes. Its complex LU/solve and decomposition cells are explicitly unsupported when the provider uses distinct _Float64x, whose generic complex division can silently overflow; default dispatch uses the CPU implementation, but that fallback does not count as provider validation.
  • MPFR/MPC CPU matrix norms and reductions have shared probes. Their CPU GEMM kernel accepts exact mode only; finite-precision GEMM uses MPLAPACK. Square solves, reusable LU and log determinants have independent CPU and MPLAPACK probes. CPU probes run with MPLAPACK disabled too, inspecting factor, solution and determinant precision at both 128 and 256 bits. Runtime LU integration tests also cover mixed element precisions, exact inputs converted to finite precision, row swaps, multiple RHS columns, factor reuse across providers, failure statuses, views and concurrent Async solves at different working precisions.
  • Float80 and MPFR/MPC projected LAPACK wrappers and end-to-end Krylov solvers remain unsupported here. Helper-level runtime epsilon tests do not establish Lanczos, Arnoldi, or exponential-action support.
  • Shared coverage must expand to further matrix functions, iterative stopping/restart criteria, async precision propagation, and other numerical subsystems before their scalar support is certified. Existing specialized tests supply additional evidence but must not be mistaken for completed matrix cells.

New numerical operations should join these cases as they are implemented. The dedicated MPFR-only special/extension probes run at 128 and 256 bits. These aggregate probes are inapplicable to native types, which have separate family probes. Missing platform providers remain explicitly unavailable; enabling MPFR never implements a native operation. Pi-scaled trigonometric functions have separate PiScaledIdentities and PiScaledAccuracy probes at every configured real precision; PiScaledPortableAccuracy directly exercises the native fallback, independently of whether a libc provider is selected. The native Gamma/error-function probes run at all configured real precisions; their high-precision oracle requires MPFR. Beta has standard-provider probes whose native binary128 cells remain unsupported. Dedicated Beta probes cover its binary128 and portable paths. ZetaIdentities and ZetaAccuracy cover every configured real type, including reflection near trivial zeros and final overflow boundaries. ExpintIdentities and ExpintAccuracy cover every configured real type and record bound fractions. Constant probes also cover configured binary128. DilogIdentities and DilogAccuracy cover all configured native reals and MPFR128/256, recording the largest fraction of the documented error bound and the corresponding input. UpperGammaIdentities and UpperGammaAccuracy cover finite shape/endpoint regions, endpoint conventions, large-parameter cancellation and final range boundaries. Their error scale is specified in the upper Gamma acceptance table. AiryIdentities and AiryAccuracy cover zero neighborhoods, region joins, positive underflow and original-input phase at huge negative arguments. The Airy acceptance table uses an amplitude-scaled absolute measure on the negative axis; its reference precision grows with the argument exponent to resolve the phase. Native identity probes remain available in MPFR-disabled configurations; independent high-precision accuracy probes require the MPFR oracle.

Keep operation-specific unsupported reasons explicit in the registry. When a combination gains support, enable the same mathematical probe for it.

Exponential extensions have separate ExponentialExtensionIdentities, ExponentialExtensionAccuracy and ExponentialExtensionPortableAccuracy probes. Accuracy probes record per-function maximum grid errors and their inputs as XML properties. See the acceptance table for the measure, targets and separation of observed maxima from guarantees.

Compound has CompoundIdentities, CompoundAccuracy and CompoundPortableAccuracy probes. The accuracy probes record the maximum reference-grid error, original input and integer exponent. Native and MPFR paths use the same represented-input reference; the native portable path is also tested directly. Selected 512-bit oracle results are repeated at 1024 bits. See native compound acceptance.

ReciprocalFunctionIdentities and ReciprocalFunctionAccuracy cover reciprocal trigonometric/hyperbolic functions across configured real types. Native Catalan is included in MathConstantsAccuracy, with correct-rounding checks and an independent Euler-transformed series. The shared oracle has a dedicated regression for negative zero, infinities and NaNs. See the acceptance table.

AgmIdentities and AgmAccuracy validate scaled arithmetic-geometric means through the shared cases. Accuracy includes an MPFR512 oracle, a 1024-bit check of the worst sampled pair and a separate elliptic-integral power series. See the AGM acceptance criteria.

FactorialIdentities and FactorialAccuracy use independently accumulated exact integers. Every finite native order is tested, plus the first overflowing one; MPFR cases use a precision-dependent finite order range. Correct rounding observed in a scan is distinct from the one-spacing-unit acceptance criterion. See typed factorials.

DigammaIdentities and DigammaAccuracy cover the native and MPFR cases. The accuracy measure is absolute near zeros and relative away from them, with an explicit logarithmic scale for reflection cancellation. The shared probes include zero/pole neighbors, full native exponent ranges, independent integer and half-integer identities, and a 1024-bit check of the worst input. See native digamma acceptance.