Skip to content

fix(simd,vml): simd_exp_f32 and vstanh within 1 ULP; tanh keeps ±0 and subnormals - #369

Merged
AdaWorldAPI merged 5 commits into
masterfrom
fix/vstanh-small-arg
Oct 24, 2026
Merged

AdaWorldAPI merged 5 commits into
masterfrom
fix/vstanh-small-arg

Conversation

@AdaWorldAPI

@AdaWorldAPI AdaWorldAPI commented Oct 24, 2026 •

Copy link
Copy Markdown
Owner

Why

ONNX Tanh is meant to map directly onto ndarray::hpc::vml::vstanh, so vstanh has to be accurate. Measured against the correctly rounded result (f64 tanh/exp rounded once to f32), it was not:

before after
vstanh, |x| in 1e-38 … 1e-10 786,163,491 ULP; every value returned 0 0 ULP
vstanh, |x| in 1e-10 … 1e-4 868,220,004 ULP; 61,186 of 131,072 returned 0 0 ULP
vstanh, |x| in 1e-4 … 0.1 14,256 ULP 1 ULP
vstanh, |x| in 0.1 … 4 111 ULP 1 ULP
vstanh(-0.0) +0.0 -0.0
vstanh on subnormals 0 the input itself
simd_exp_f32, [-87, 88.5] 39 ULP (the doc claimed ~2) 1 ULP

Causes

  1. vstanh cancelled. It computed 2·sigmoid(2x) − 1, and below |x| ≈ 1e-7 the sigmoid rounds to exactly 0.5. The kernel now works on a = |x|:

    • a + a·z·P(z) with z = a² for a < 0.625 (the Cephes tanhf polynomial);
    • 1 − 2/(exp(2a) + 1) otherwise;
    • negated where x < 0, and x returned where x == 0.

    It uses only F32x16 compares and select, because the NEON U32x16 has no &/|. The tail and the strided path run through the same 16-lane kernel, so a value gets the same bits wherever it sits in the array.

  2. simd_exp_f32 used a degree-5 Taylor series. Its r⁶/720 term is about 2.4e-6 relative at |r| = ln2/2. It now uses the Cephes expf minimax 1 + r + r²·P(r), still with plain * and +, so the bit-for-bit scalar reference keeps every backend identical. Every user of simd_exp_f32 (vsexp, softmax, GELU, the model layers) gets the more accurate exp too.

No public API changes and no new public functions.

Evidence

  • Tests first. Commits 0e176e1 (the tanh gate) and 9650bc5 (the exp gate) failed on the old kernels before the fix landed.

    • Tanh gate: ≤ 2 ULP over 8 bands, both signs, reaching the tail. Special values are exact in chunk and tail, and tail == chunk bit for bit.
    • Exp gate: ≤ 2 ULP.
    • 4f6c132 adds four more tests: every length from 0 to 65 agrees with a single lane; strided and 2-D views match the contiguous result; the neighbourhood of the 0.625 branch point is monotone and odd; saturation and extremes behave.
  • Exhaustive sweep (examples/vstanh_exhaustive.rs, all 2³² f32 inputs, AVX-512):

    old kernel new kernel
    max ULP 868,220,004 1, in every magnitude band
    nonzero inputs returned as 0 1,728,053,246 0
    odd-symmetry violations 975,101,106 0
    outside [−1, 1] / non-monotone — 0 / 0

    NaN is preserved. The oracle (f64 tanh rounded once) agrees with mpmath's correctly rounded value on 200,000 of 200,000 sampled inputs, and on every worst case. Exactly one input lies near a rounding boundary in f64; it also agrees with mpmath.

  • Throughput, same run: vstanh 3.06 ns/elem, the old kernel 2.46 ns/elem, scalar f32::tanh 11.66 ns/elem.

  • Disable runs. Each was applied with an asserted anchor, and each turned its tests red:

    Disabled Red tests
    small-argument branch ULP sweep, special values
    sign restore 4 tests
    only the x == 0 select special values (-0.0 → +0.0)
    libm in the tail tail == chunk; every length
    old Taylor exp exp ULP, bit-for-bit reference, tanh sweep
  • Backends. The same ULP probe gives identical numbers on four builds:

    • AVX-512 (x86-64-v4);
    • AVX2 (x86-64-v3);
    • the x86-64 baseline;
    • aarch64 NEON under qemu-aarch64-static.

    That is exp ≤ 1 ULP on 4,194,304 points, and tanh ≤ 1 ULP on 2·2¹⁸ points per band.

  • wasm. scripts/wasm-parity.sh passes. Its harness restated the old exp polynomial as its bit-for-bit reference and now restates the new one.

  • Native gates. cargo test --release --lib: 2590 passed, 0 failed. cargo clippy --lib --tests -- -D warnings and cargo fmt --check are clean.

Not covered

  • vstanh accuracy on wasm. The harness checks exp there, not tanh.
  • The nightly core::simd arm.
  • The exhaustive sweep on the non-AVX-512 tiers. Those were only sampled.

Downstream: the error decomposition in AdaWorldAPI/lance-graph#1757 shows the tanh stage adding about 1e-8 to a Fisher-Z round trip, against about 1e-3 for i8 quantization.


Generated by Claude Code

…/inf/NaN, tail == chunk

Fails on the current kernel: 2·sigmoid(2x) − 1 cancels for |x| below
~1e-7 (785,499,064 ULP at x = -9.5e-11), and -0.0 returns +0.0.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
Fails today: the degree-5 Taylor polynomial is 39 ULP off at reduced
argument |r| = ln2/2 (r^6/720 ≈ 2.4e-6 relative), against the doc
comment's '~2 ULP'.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
simd_exp_f32: replace the degree-5 Taylor polynomial (39 ULP at reduced
argument |r| = ln2/2) with the Cephes expf minimax 1 + r + r²·P(r).
Plain * and + as before, so the scalar reference still matches bit for
bit on every backend. Measured: max 1 ULP on 4,194,304 points over
[-87, 88.5], against exp in f64 rounded once (AVX-512 and AVX2).

vstanh: 2·sigmoid(2x) − 1 cancelled for small |x| (every |x| below
~1e-7 returned 0, up to 868,220,004 ULP) and turned -0.0 into +0.0.
The 16-lane kernel now works on |x|:
- |x| < 0.625: a + a·z·P(z), z = a² (Cephes tanhf polynomial);
- otherwise 1 − 2/(exp(2a) + 1);
then copies the sign bit of x back. The tail and the strided path go
through the same kernel, so a value gets the same bits anywhere in the
array. Measured: max 1 ULP in every band from 1e-38 to 100, both signs;
±0, subnormals, ±inf and NaN exact.

Users of simd_exp_f32 (vsexp, softmax, GELU, layer kernels) get the more
accurate exp too. Full lib suite: 2590 passed, 0 failed.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
…p reference

The NEON U32x16 has no & / |, so the sign-bit copy did not compile for
aarch64. tanh now negates where x < 0 and returns x where x == 0 (keeps
-0.0), using only F32x16 compares and select, which every backend has.
Measured under qemu-aarch64 (NEON): exp max 1 ULP, tanh max 1 ULP in every
band, the same as AVX-512 and AVX2.

wasm-simd-parity restated the old Taylor polynomial as its bit-for-bit
exp reference; it now restates the Cephes coefficients. wasm parity: OK.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
@coderabbitai

coderabbitai Bot commented Oct 24, 2026 •

Copy link
Copy Markdown
Contributor

Important

  • 🔍 Trigger review

This repository does not receive automatic reviews because it has fewer than 10 stars.

Run configuration
  • Source mode: Auto
  • Repository: Unavailable
  • Branch: Unavailable
  • File: No configuration file from this repository applied
  • Commit: Unavailable
  • Layers: Organization UI
  • Review profile: CHILL
  • Plan: Advanced
  • Run ID: 6d36b628-8aea-46c1-94b5-8caccc28ea62
  • Autofix · Keep fixing CodeRabbit findings and required CI, and resolving merge conflicts

Warning

Billing warning: we have not been able to collect payment for this subscription for more than 72 hours. Please update the payment method or pay any pending invoices in Billing to avoid service interruption.

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands.

…saturation

examples/vstanh_exhaustive.rs evaluates every f32 (each non-negative bit
pattern with its negation) against tanh in f64 rounded once, and runs the
previous 2*sigmoid(2x)-1 kernel as a baseline arm. Measured on AVX-512:
- repaired: max 1 ULP over all 2^32 inputs; 0 zeroed nonzero inputs;
  0 bitwise odd-symmetry violations; 0 outside [-1, 1]; 0 monotonicity
  violations; every NaN stays NaN; 1 reference flip-risk input.
- previous kernel: max 868,220,927 ULP; 1,728,053,246 nonzero inputs
  returned 0; 975,101,106 odd-symmetry violations.
- 3.06 ns/elem (old kernel 2.46, scalar f32::tanh 11.66).
The flip-risk input and every band's worst input agree with mpmath at 60
digits, and the f64 reference equals mpmath correct rounding on 200,000
log-uniform inputs (0 differ).

New unit tests: every length 0..65 around the 16-lane chunk matches the
one-element result; 2-D, transposed and strided views match contiguous;
4001 consecutive floats across the 0.625 branch threshold, both signs,
within 2 ULP, monotone and odd; MAX/-MAX and large inputs saturate to ±1.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
@chatgpt-codex-connector

chatgpt-codex-connector Bot commented Oct 24, 2026 •

Copy link
Copy Markdown

Codex Review Summary

This comment shows the latest Codex review activity on this pull request.

Review Status Commit Review trigger
📝 Code Review ✅ Completed 2026-10-24T08:33:29.730994Z 4f6c132 Draft marked ready
ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review" or "@codex security review".

Codex reacts with 👀 while any review is running, comments if it has suggestions, and reacts with 👍 once all reviews finish with no findings.

@AdaWorldAPI
AdaWorldAPI merged commit 7ef88e4 into master Oct 24, 2026
30 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants