Float64 accuracy of std.math on CPU, measured against the system libm (Mojo 1.0.0)

Mojo Standard Library
, ,
1

Float64 accuracy of std.math on CPU, measured against the system libm (Mojo 1.0.0)

Hello,

A few months ago, on a 1.0 beta, I ran into a small numerical discrepancy in exp/log while gradient-checking some derivatives I had worked out by hand. I never chased it down properly. Recently I asked an AI assistant (Claude Opus 5) to look into it seriously on 1.0.0. It wrote the measurement harness, ran a ~250k-point survey in its own sandbox, and produced the numbers below.

To be clear about what is mine and what isn’t:

  • Section 1 I ran myself, on my own machine. It is a single self-contained file, and my output was byte-identical to the assistant’s, on a different machine and a different glibc patch level.
  • Everything else (the 250k-point survey, the tables in sections 2 and 3, the issue search in section 4) comes from the assistant’s runs. I have not independently reproduced any of it.

I don’t know the constraints behind the current implementations, so some or all of this may be known, intended, or measured the wrong way. Apologies if the post is out of place; I only want to help or inform.

TL;DR

functions
:white_check_mark: identical to libm, 0 ULP, every point sin cos tan asin acos atan atan2 log10 expm1 sqrt cbrt hypot
:warning: relative difference ~1e-10 exp exp2 log log2 cosh x ** 0.5
:warning: relative difference ~1e-8 / ~1e-6 x ** y / log1p
:cross_mark: different result class than libm exp → 0.0, sinh → nan, cosh → inf, log of a subnormal

1. Self-contained reproducer

One file, no dependencies. Both columns are computed in the same process: the left one through std.math, the right one through a direct FFI call into the libm the process is already linked against. The three control rows at the top exist to show the comparison itself works.

# std.math vs the system libm, Float64, CPU.  Run with:  mojo run libm_check.mojo
# "bitwise" compares the two IEEE-754 bit patterns, so NaN compares equal to NaN.

from std.math import exp, log, log1p, sin, cos, sinh, cosh, sqrt
from std.ffi import external_call


def pad(s: String, n: Int) -> String:
    var o = s
    while o.byte_length() < n:
        o += " "
    return o


def row(label: String, mojo: Float64, libm: Float64):
    var flag = String("DIFF")
    if mojo.to_bits() == libm.to_bits():
        flag = String("match")
    print(pad(label, 22), pad(String(mojo), 25), pad(String(libm), 25), flag)


def main():
    var bar = String("")
    for _ in range(84):
        bar += "-"

    print(pad(String("expression"), 22), pad(String("std.math"), 25),
          pad(String("system libm"), 25), "bitwise")
    print(bar)

    # controls: expected to agree
    row("sin(1.0)", sin(Float64(1.0)), external_call["sin", Float64](Float64(1.0)))
    row("cos(1e22)", cos(Float64(1e22)), external_call["cos", Float64](Float64(1e22)))
    row("sqrt(2.0)", sqrt(Float64(2.0)), external_call["sqrt", Float64](Float64(2.0)))
    print(bar)

    # rounding differences
    row("exp(8.0)", exp(Float64(8.0)), external_call["exp", Float64](Float64(8.0)))
    row("2.0 ** 0.5", Float64(2.0) ** Float64(0.5),
        external_call["pow", Float64](Float64(2.0), Float64(0.5)))
    row("log1p(-0.25)", log1p(Float64(-0.25)),
        external_call["log1p", Float64](Float64(-0.25)))
    print(bar)

    # different result class
    row("exp(-708.745)", exp(Float64(-708.745)),
        external_call["exp", Float64](Float64(-708.745)))
    row("exp(-745.0)", exp(Float64(-745.0)),
        external_call["exp", Float64](Float64(-745.0)))
    row("sinh(709.783)", sinh(Float64(709.783)),
        external_call["sinh", Float64](Float64(709.783)))
    row("cosh(709.4365)", cosh(Float64(709.4365)),
        external_call["cosh", Float64](Float64(709.4365)))
    row("log(5e-324)", log(Float64(5e-324)),
        external_call["log", Float64](Float64(5e-324)))

Output on my machine (Mojo 1.0.0, Ubuntu 24.04.3 under WSL2, glibc 2.39):

expression             std.math                  system libm               bitwise
------------------------------------------------------------------------------------
sin(1.0)               0.8414709848078965        0.8414709848078965        match
cos(1e22)              0.523214785395139         0.523214785395139         match
sqrt(2.0)              1.4142135623730951        1.4142135623730951        match
------------------------------------------------------------------------------------
exp(8.0)               2980.9579870518287        2980.9579870417283        DIFF
2.0 ** 0.5             1.4142135623734946        1.4142135623730951        DIFF
log1p(-0.25)           -0.2876821432670699       -0.2876820724517809       DIFF
------------------------------------------------------------------------------------
exp(-708.745)          0.0                       1.570208859696427e-308    DIFF
exp(-745.0)            0.0                       5e-324                    DIFF
sinh(709.783)          nan                       8.991046692770538e+307    DIFF
cosh(709.4365)         inf                       6.358097963416392e+307    DIFF
log(5e-324)            -709.0895657128241        -744.4400719213812        DIFF

Reading the last block:

  • exp returns exactly zero from x ≤ -708.745 down. libm returns real subnormals to x ≈ -745.13.
  • sinh returns NaN from x ≈ 709.783, where libm still returns a finite ~9e307.
  • cosh overflows to inf from x ≈ 709.4365; libm stays finite to x ≈ 710.4758.
  • Subnormal inputs to log/log2 look flushed to zero: log(5e-324) is off by 35.35 in absolute terms, log2 gives -1023 instead of -1074, and 5e-324 ** 0.5 is off by a factor of ~5e7.

The sinh and cosh thresholds are close together, so both may go through a shared exp-based path. Neither of us looked at the implementation.

2. Rounding differences

Normal inputs and normal outputs only, so the cases above are excluded.

function % matching libm max ULP max rel. diff
exp 0.19% 2.6e6 2.9e-10
exp2 0.84% 1.8e6 2.1e-10
log 3.96% 5.5e6 6.2e-10
log2 4.26% 4.0e6 6.2e-10
cosh 0.14% 2.6e6 2.9e-10
x ** 0.5 0.10% 4.5e6 6.4e-10
x ** y 0.05% 8.8e7 1.2e-8
log1p 75.39% 6.9e9 1.1e-6
sinh 89.68% 2 3.3e-16
tanh 98.16% 2 2.2e-16
  • sinh and tanh are fine on rounding. sinh’s only issue is the NaN above.
  • log1p is the outlier: accurate for tiny arguments, not in the ordinary range.
  • x ** 0.5 differs from sqrt(x) even where the exact result is representable, as the reproducer shows.

3. How the survey was run

  • Mojo 1.0.0, x86-64 Linux, Float64, scalar, no vectorisation.
  • ~250k points: seeded pseudo-random samples over many decades, plus hand-picked ones (zero, subnormals, overflow/underflow thresholds, near the poles of tan, sin(1e22)). sqrt also got 220k raw bit patterns covering every exponent, 20k of them subnormal.
  • Inputs and outputs carried as exact IEEE-754 bit patterns. ULP distance = difference of the binary representations under the usual monotone ordering, never computed in floating point.
  • Reference: the system libm (glibc), reached through CPython’s math for the scoring, and directly through FFI in the reproducer above. Both agree.

4. Already on record

  • #1584 accuracy loss in math.exp, open since Dec 2023
  • #1303 accuracy loss in math.log, open since Nov 2023
  • #1235 accuracy loss in math.exp2, open since Nov 2023
  • #6082 fixed tanh for float64 via an expm1-based formula, merged March 2026 — consistent with the measurements above: expm1 matched libm everywhere tested, tanh is at 2 ULP

No issue turned up on the three result-class cases or on log1p, but the search may well have missed them.

There’s also an existing thread here on MPFR-based correctness testing for Mojo, which is a far more rigorous approach: MPFR as a correctly rounded reference with per-rounding-mode checks, as used by CORE-MATH and LLVM-libc. Better foundation than this, if anyone wants to take it further.

5. Caveats

  • glibc is not a proof of correct rounding. The table in section 2 measures divergence from the system libm, not incorrectness in the IEEE sense. It would need an MPFR re-run before anything there could be called an error bound. The result-class cases in section 1 don’t depend on the reference.
  • Two machines, but effectively one platform: both x86-64 Linux with glibc 2.39. No data for macOS, arm64, or accelerators, and the GPU picture likely differs.
  • The sample is not exhaustive. 0 ULP above means “nothing found”, not “proven correct”.
  • It’s entirely possible this reflects a deliberate speed/accuracy trade-off I’m not aware of.

6. Workaround, in case anyone is in the same spot

Calling libm directly works and costs little:

from std.ffi import external_call

def exp_c(x: Float64) -> Float64:
    return external_call["exp", Float64](x)

Roughly +19% per transcendental call in scalar code, and identical to libm on the 8k points checked for exp, log, cosh and pow, subnormals and boundaries included.

Why this level of accuracy matters in my case (click to expand)

Float64 on CPU, scalar and deterministic. Points and frames on the unit sphere (normalisation, dot products, rotations by an angle), softmax / log-sum-exp over scores, cross-entropy-style objectives on the resulting log-probabilities.

Two things made elementary-function accuracy load-bearing rather than cosmetic:

  • Gradient checks. Analytic gradients are derived and implemented by hand, then validated against central finite differences (L(x+h) - L(x-h)) / (2h). What certifies the derivative isn’t the size of the error but its rate: it should fall as O(h²) over a couple of decades before flattening onto the cancellation floor, roughly eps * |L| / h. With correctly rounded functions and h ≈ 1e-5 that floor is near 1e-11 and the ramp is visible. A relative error of 1e-10 in exp or log raises the floor enough that the O(h²) regime disappears, and the check can no longer separate a wrong derivative from library noise. That is what put me onto this in the first place, months ago.
  • Subnormal probabilities. exp(s - max s) reaches the subnormal range routinely once a distribution becomes confident. A subnormal probability is fine; one flushed to exactly zero is not, because log(0) then propagates into the gradient.

Thanks for reading. The survey harness is small (a Mojo program emitting input/output bit patterns, a Python script scoring them) and I can share it if it’s of any use.

Where section 1 was run (mine):

Mojo:    1.0.0 (ed45d567), installed with `uv pip install mojo==1.0.0`
OS:      Ubuntu 24.04.3 LTS on WSL2, x86-64
libc:    glibc 2.39 (2.39-0ubuntu8.8)

Where sections 2 and 3 were run (the assistant’s sandbox):

Mojo:    1.0.0 (ed45d567), installed from PyPI
OS:      Ubuntu 24.04.4 LTS, x86-64
libc:    glibc 2.39 (2.39-0ubuntu8.7)
Python:  3.12.3 (only used to score results against the system libm)
2

Please create an issue in modular github repository: GitHub - modular/modular: The Modular Platform (includes MAX & Mojo) · GitHub

3

Filed as [BUG]: Float64 std.math on CPU: exp flushes subnormal results to zero, sinh returns NaN, cosh overflows early, log1p loses ~1e-6 · Issue #7068 · modular/modular · GitHub.
Thanks.

I kept it to the cases that aren’t already tracked in #1584 / #1303 / #1235.