ULP Accuracy Measurement Tool

Standalone tool for measuring the accuracy of Eigen‘s vectorized math functions in units of ULP (Unit in the Last Place). Compares Eigen’s SIMD implementations against either MPFR (128-bit high-precision reference) or the standard C++ math library.

Building

From the Eigen build directory:

cd build
cmake ..
cmake --build . --target ulp_accuracy

If MPFR and GMP are installed, the build automatically enables MPFR support (EIGEN_HAS_MPFR). Without them, only --ref=std is available.

Installing MPFR (Debian/Ubuntu)

sudo apt install libmpfr-dev libgmp-dev

Usage

./test/ulp_accuracy [options]

Options:
  --func=NAME    Function to test (required unless --list)
  --lo=VAL       Start of range (default: -inf)
  --hi=VAL       End of range (default: +inf)
  --double       Test double precision (default: float)
  --step=EPS     Sampling step: advance by (1+EPS)*nextafter(x)
                 (default: 0 = exhaustive; useful for double, e.g. 1e-6)
  --threads=N    Number of threads (default: all cores)
  --batch=N      Batch size for Eigen eval (default: 4096)
  --ref=MODE     Reference: 'std' (default) or 'mpfr'
  --hist_width=N Histogram half-width in ULPs (default: 10)
  --list         List available functions

Examples

List all supported functions:

./test/ulp_accuracy --list

Exhaustive float test of sin against std (tests all ~4.28 billion finite floats):

./test/ulp_accuracy --func=sin

Float test against MPFR (more accurate reference, but slower):

./test/ulp_accuracy --func=sin --ref=mpfr

Double precision test with geometric sampling (exhaustive is impractical for double):

./test/ulp_accuracy --func=exp --double --step=1e-6

Test a specific range:

./test/ulp_accuracy --func=sin --lo=0 --hi=6.2832

Output

The tool prints:

• Test configuration: function, range, reference mode, thread count
• Max |ULP error|: worst-case absolute ULP error with the offending input value
• Mean |ULP error|: average absolute ULP error across all tested values
• Signed ULP histogram: distribution of signed errors showing bias direction

Example output:

Function: sin (float)
Range: [-inf, inf]
Representable values in range: 4278190082
Reference: MPFR (128-bit)
Threads: 32
Batch size: 4096

Results:
  Values tested: 4278190081
  Time: 529.04 seconds (8.1 Mvalues/s)
  Max |ULP error|: 2
    at x = -1.5413464e+38 (Eigen=-0.482218683, ref=-0.482218742)
  Mean |ULP error|: 0.0874

Signed ULP error histogram [-10, +10]:
  -2   :        51988 (  0.001%)
  -1   :    186805349 (  4.366%)
  0    :   3904475407 ( 91.265%)
  1    :    186805349 (  4.366%)
  2    :        51988 (  0.001%)

How it works

1. Range splitting: The input range is divided evenly across threads by splitting the linear ULP space.

2. Batched evaluation: Each thread fills batches of input values, evaluates them through Eigen's vectorized path (using Eigen::Array operations), and computes reference values one at a time.

3. ULP computation: IEEE 754 bit patterns are mapped to a linear integer scale where adjacent representable values are adjacent integers. The signed ULP error is the difference between Eigen's result and the reference on this scale. Special cases (NaN, infinity mismatches) report infinite error.

4. Result reduction: Per-thread statistics (max error, mean error, histogram) are merged after all threads complete.

Supported functions

Category	Functions
Trigonometric	sin, cos, tan, asin, acos, atan
Hyperbolic	sinh, cosh, tanh, asinh, acosh, atanh
Exponential/Log	exp, exp2, expm1, log, log1p, log10, log2
Error/Gamma	erf, erfc, lgamma
Other	logistic, sqrt, cbrt, rsqrt

File organization

• ulp_accuracy.cpp — Main tool: ULP computation, worker threads, CLI, result printing
• mpfr_reference.h — MPFR reference function wrappers and scalar conversion helpers

Performance tips

• Float exhaustive sweeps test ~4.28 billion values. With --ref=std this takes ~50 seconds per function; with --ref=mpfr it takes ~500 seconds (10x slower).
• For double precision, exhaustive testing is impractical. Use --step=1e-6 to sample ~2.88 billion values geometrically.
• Thread count defaults to all available cores. MPFR is the bottleneck (single MPFR call per value per thread), so more cores help significantly.