logaritmisk/trueskill-tt
1
0
Fork 0
Code Issues 1 Pull Requests Packages Projects Releases Activity

erfc is a 1.2e-7 approximation and sets the accuracy floor for the whole crate #41

New Issue
Closed
opened 2026-09-07 13:29:12 +00:00 by logaritmisk · 2 comments
logaritmisk commented 2026-09-07 13:29:12 +00:00
Owner
Copy Link
Copy Source

Deferred from the tail-precision work in 7341669. Not urgent, and deliberately not bundled with that fix — it is a different kind of change with a different risk profile.

What it is

erfc (src/lib.rs) is the Numerical Recipes erfcc routine: a single rational-in-t expression wrapped in one exp, with a documented fractional error of 1.2e-7. Everything numeric in the crate goes through it — cdf, sf, pdf ratios, ppf, compute_margin, the truncation moments, log_evidence, and both prediction paths.

The good news, measured against an incomplete-gamma reference accurate to ~1e-15, is that it is well-behaved: the error is genuinely relative, and it does not degrade in the tail.

       x          erfc(x) [crate]      erfc(x) [reference]      rel err
     0.5   4.79500092276757439e-1   4.79500122186975797e-1      6.24e-8
     2.0   4.67773498933404374e-3   4.67773498104706413e-3      1.77e-9
     8.0  1.12242980240431987e-29  1.12242971729824276e-29      7.58e-8
    20.0 5.39586506702822336e-176 5.39586561160751332e-176      1.01e-7
    26.0 5.66319185484591015e-296 5.66319240885590743e-296      9.78e-8

That property is load-bearing — it is why 7341669 could fix the tail bugs without touching erfc at all. The bugs were cancellation in the callers, not error in erfc.

Why it is still worth revisiting

1. It sets a floor on every identity the crate can assert. The approximation is not exactly antisymmetric, so mathematical identities hold only to ~1e-7:

  • erfc(z) + erfc(-z) differs from 2 by ~3e-8, so sf(z) + cdf(z) differs from 1 by ~3e-8.
  • erfc(0) returns 1.000_000_03 rather than 1, so cdf(mu, mu, sigma) is not exactly 0.5.

Two tests in src/lib.rs (survival_function_matches_the_naive_form_where_that_form_works, survival_and_cdf_partition_the_mass) are asserted at 1e-6 purely because of this. They would be exact to ~1e-15 with a better erfc, and a tighter identity test is a better regression net.

2. It is the residual in the new prediction code. The adaptive quadrature in predict_win_probabilities converges to ~1e-15 against the exact two-team closed form, and the ranking grid converges as O(h^2) — but both floor out at ~2.4e-8, which is erfc, not the integrator. MIN_GRID_POINTS is documented as sized to that floor; a better erfc would let it drop or let the accuracy rise.

What it is not

Worth being explicit, so this does not get overstated:

  • It is not currently causing wrong answers. Measured contribution to a finished prediction probability is ~6e-9. Skill posteriors and rankings are unaffected at any precision anyone cares about.
  • It was not the cause of the tail bugs. Those were cancellation in cavity_evidence and v_w, fixed in 7341669 without changing erfc.

Options

  1. Cody's rational Chebyshev approximation — the standard replacement, ~1e-16, three branches by magnitude, no dependency. Perhaps 60 lines. Probably 2-3x the arithmetic of erfcc.
  2. A dependency (libm, statrs, puruspe). Accurate and maintained, but the crate is deliberately dependency-light and #![forbid(unsafe_code)], so each candidate needs checking against both.
  3. Leave it. Defensible: the error is relative, tail-stable, and demonstrably below anything that changes a rating.

The catch, and why this needs measuring rather than deciding

erfc sits in the inference hot loop — TruncFactor::propagate calls it for every rank-adjacent pair, every event, every iteration, on every sweep. A 2-3x slowdown there is a real cost against a benefit currently worth ~6e-9.

So this should not be done on principle. It should be done only if benches/ shows the cost is small, or if a future feature needs identities to hold tighter than 1e-7. The erfcx added in 7341669 already sidesteps erfc entirely for x >= 2 via a continued fraction, and is accurate to ~1e-15 there — so part of the hot path has quietly improved already.

Acceptance, if picked up

  • erfc accurate to ~1e-15 relative, verified against an independent reference across [0, 27] including the far tail.
  • erfc(0) == 1.0 exactly, and erfc(z) + erfc(-z) == 2.0 to ~1e-15.
  • The two identity tests above tightened from 1e-6 to 1e-12 or better.
  • cargo bench before and after, with the inference-loop delta recorded in the PR — and the change abandoned if it is not worth paying.
Deferred from the tail-precision work in `7341669`. Not urgent, and deliberately **not** bundled with that fix — it is a different kind of change with a different risk profile. ## What it is `erfc` (`src/lib.rs`) is the Numerical Recipes `erfcc` routine: a single rational-in-`t` expression wrapped in one `exp`, with a documented fractional error of **1.2e-7**. Everything numeric in the crate goes through it — `cdf`, `sf`, `pdf` ratios, `ppf`, `compute_margin`, the truncation moments, `log_evidence`, and both prediction paths. The good news, measured against an incomplete-gamma reference accurate to ~1e-15, is that it is **well-behaved**: the error is genuinely *relative*, and it does not degrade in the tail. ``` x erfc(x) [crate] erfc(x) [reference] rel err 0.5 4.79500092276757439e-1 4.79500122186975797e-1 6.24e-8 2.0 4.67773498933404374e-3 4.67773498104706413e-3 1.77e-9 8.0 1.12242980240431987e-29 1.12242971729824276e-29 7.58e-8 20.0 5.39586506702822336e-176 5.39586561160751332e-176 1.01e-7 26.0 5.66319185484591015e-296 5.66319240885590743e-296 9.78e-8 ``` That property is load-bearing — it is why `7341669` could fix the tail bugs without touching `erfc` at all. The bugs were cancellation in the *callers*, not error in `erfc`. ## Why it is still worth revisiting **1. It sets a floor on every identity the crate can assert.** The approximation is not exactly antisymmetric, so mathematical identities hold only to ~1e-7: - `erfc(z) + erfc(-z)` differs from `2` by ~3e-8, so `sf(z) + cdf(z)` differs from `1` by ~3e-8. - `erfc(0)` returns `1.000_000_03` rather than `1`, so `cdf(mu, mu, sigma)` is not exactly `0.5`. Two tests in `src/lib.rs` (`survival_function_matches_the_naive_form_where_that_form_works`, `survival_and_cdf_partition_the_mass`) are asserted at `1e-6` purely because of this. They would be exact to ~1e-15 with a better `erfc`, and a tighter identity test is a better regression net. **2. It is the residual in the new prediction code.** The adaptive quadrature in `predict_win_probabilities` converges to ~1e-15 against the exact two-team closed form, and the ranking grid converges as `O(h^2)` — but both floor out at **~2.4e-8**, which is `erfc`, not the integrator. `MIN_GRID_POINTS` is documented as sized to that floor; a better `erfc` would let it drop or let the accuracy rise. ## What it is *not* Worth being explicit, so this does not get overstated: - **It is not currently causing wrong answers.** Measured contribution to a finished prediction probability is ~6e-9. Skill posteriors and rankings are unaffected at any precision anyone cares about. - **It was not the cause of the tail bugs.** Those were cancellation in `cavity_evidence` and `v_w`, fixed in `7341669` without changing `erfc`. ## Options 1. **Cody's rational Chebyshev approximation** — the standard replacement, ~1e-16, three branches by magnitude, no dependency. Perhaps 60 lines. Probably 2-3x the arithmetic of `erfcc`. 2. **A dependency** (`libm`, `statrs`, `puruspe`). Accurate and maintained, but the crate is deliberately dependency-light and `#![forbid(unsafe_code)]`, so each candidate needs checking against both. 3. **Leave it.** Defensible: the error is relative, tail-stable, and demonstrably below anything that changes a rating. ## The catch, and why this needs measuring rather than deciding `erfc` sits in the **inference hot loop** — `TruncFactor::propagate` calls it for every rank-adjacent pair, every event, every iteration, on every sweep. A 2-3x slowdown there is a real cost against a benefit currently worth ~6e-9. So this should not be done on principle. It should be done only if `benches/` shows the cost is small, or if a future feature needs identities to hold tighter than 1e-7. The `erfcx` added in `7341669` already sidesteps `erfc` entirely for `x >= 2` via a continued fraction, and is accurate to ~1e-15 there — so part of the hot path has quietly improved already. ## Acceptance, if picked up - `erfc` accurate to ~1e-15 relative, verified against an independent reference across `[0, 27]` including the far tail. - `erfc(0) == 1.0` exactly, and `erfc(z) + erfc(-z) == 2.0` to ~1e-15. - The two identity tests above tightened from `1e-6` to `1e-12` or better. - `cargo bench` before and after, with the inference-loop delta recorded in the PR — and the change abandoned if it is not worth paying.
logaritmisk commented 2026-09-07 13:43:20 +00:00
Author
Owner
Copy Link
Copy Source

Independent confirmation and concrete replacement candidates, from a literature/ecosystem survey run separately from the measurements above.

The measurement reproduces

The crate's erfc was ported to Python and compared against 40-digit mpmath over 20,000 points — a different reference and a different method from the incomplete-gamma comparison in the issue body:

  • worst absolute error in Φ: 4.15e-08
  • worst relative error in Φ: 1.05e-07
  • Φ(0) returns 0.5000000150…

Those match the figures above to the digit, from an independent implementation. The 1.2e-7 documented bound is real and is being hit.

One thing the issue body understates

The error compounds across CDF factors. predict_win_probabilities forms a product of N-1 normal CDFs, so at eight competitors the relative error is ~7e-7, not ~1e-7. That does not change the conclusion that nothing user-visible is affected, but it does mean the accuracy floor rises with team count rather than staying flat — worth knowing before anyone tightens a tolerance on the prediction paths.

Replacement candidates, measured

crate accuracy deps unsafe verdict
spec_math 0.1.6 ~1e-16 zero zero best fit — pure-Rust Cephes, MIT/Apache, provides erf/erfc/erf_inv/norm_cdf/norm_cdf_inv; norm_cdf(norm_cdf_inv(p)) round-trips to 2e-15 down to p=1e-300. Risk: 39.9k downloads, 10 stars
libm 0.2.16 ~1e-16 zero none in erf.rs safer bet on maintenance — rust-lang maintained
statrs ~1e-10 24 crates — do not use for erf. erfc(0.5) gives 4.79500122236e-1 where Python, libm, spec_math, puruspe and errorfunctions all agree on 4.79500122187e-1. Its erf_inv is fine; the forward function is the weak link
Cody rational Chebyshev, vendored ~1e-16 zero zero option 1 in the issue body; still viable

f64::erf in std is nightly-only behind #![feature(float_erf)] (rust-lang/rust#136321) — still open with FCP unchecked, so not a route yet.

Does not change the recommendation

This remains deferred, and for the same reason: erfc is in TruncFactor::propagate, so it runs for every rank-adjacent pair, every event, every iteration, every sweep. The acceptance criteria above already require a cargo bench delta and abandoning the change if it is not worth paying. Nothing here argues for doing it sooner — it just means that when someone does pick it up, spec_math and libm are the two to benchmark, and statrs can be skipped.

Independent confirmation and concrete replacement candidates, from a literature/ecosystem survey run separately from the measurements above. ## The measurement reproduces The crate's `erfc` was ported to Python and compared against 40-digit `mpmath` over 20,000 points — a different reference and a different method from the incomplete-gamma comparison in the issue body: - worst **absolute** error in Φ: **4.15e-08** - worst **relative** error in Φ: **1.05e-07** - `Φ(0)` returns `0.5000000150…` Those match the figures above to the digit, from an independent implementation. The 1.2e-7 documented bound is real and is being hit. ## One thing the issue body understates **The error compounds across CDF factors.** `predict_win_probabilities` forms a product of `N-1` normal CDFs, so at eight competitors the relative error is ~**7e-7**, not ~1e-7. That does not change the conclusion that nothing user-visible is affected, but it does mean the accuracy floor rises with team count rather than staying flat — worth knowing before anyone tightens a tolerance on the prediction paths. ## Replacement candidates, measured | crate | accuracy | deps | unsafe | verdict | |---|---|---|---|---| | `spec_math` 0.1.6 | ~1e-16 | zero | zero | best fit — pure-Rust Cephes, MIT/Apache, provides `erf`/`erfc`/`erf_inv`/`norm_cdf`/`norm_cdf_inv`; `norm_cdf(norm_cdf_inv(p))` round-trips to 2e-15 down to p=1e-300. Risk: 39.9k downloads, 10 stars | | `libm` 0.2.16 | ~1e-16 | zero | none in `erf.rs` | safer bet on maintenance — rust-lang maintained | | **`statrs`** | **~1e-10** | 24 crates | — | **do not use for `erf`.** `erfc(0.5)` gives `4.79500122236e-1` where Python, libm, `spec_math`, `puruspe` and `errorfunctions` all agree on `4.79500122187e-1`. Its `erf_inv` is fine; the forward function is the weak link | | Cody rational Chebyshev, vendored | ~1e-16 | zero | zero | option 1 in the issue body; still viable | `f64::erf` in std is nightly-only behind `#![feature(float_erf)]` ([rust-lang/rust#136321](https://github.com/rust-lang/rust/issues/136321)) — still open with FCP unchecked, so not a route yet. ## Does not change the recommendation This remains deferred, and for the same reason: `erfc` is in `TruncFactor::propagate`, so it runs for every rank-adjacent pair, every event, every iteration, every sweep. The acceptance criteria above already require a `cargo bench` delta and abandoning the change if it is not worth paying. Nothing here argues for doing it sooner — it just means that when someone does pick it up, `spec_math` and `libm` are the two to benchmark, and `statrs` can be skipped.
logaritmisk added the numericsperformance labels 2026-09-07 13:52:11 +00:00
logaritmisk commented 2026-09-07 20:25:36 +00:00
Author
Owner
Copy Link
Copy Source

Done in 3dd6593. The answer to "is it worth paying for" turned out to be that there is nothing to pay.

Accuracy

libm (the Rust port of FDLIBM), measured against an independent incomplete-gamma reference over 20 000 points per range:

range previous (NR erfcc) libm
[-3, 0] 7.95e-8 2.15e-14
[0, 0.5] 8.69e-8 4.70e-14
[0.5, 2] 9.38e-8 1.24e-12
[2, 6] 1.04e-7 5.20e-14
[6, 26] 1.07e-7 1.75e-13

And the identities the issue flagged:

  • erfc(0) returns exactly 1.0, where it was 1.000_000_03.
  • abs(erfc(z) + erfc(-z) - 2) is 2.22e-16 — one ULP — where it was 6.00e-8.

Performance: no change detected

benches/batch.rs, against a baseline taken immediately before the swap:

Batch::iteration   change: [-3.3423% -0.2286% +2.2576%] (p = 0.89 > 0.05)
                   No change in performance detected.

The acceptance criteria asked for the benchmark delta and for the change to be abandoned if it was not worth paying. There is no delta to weigh.

Why it is free, and a correction to the issue's reasoning

This issue assumed a 2-3× arithmetic cost. A naive microbenchmark agrees — sweeping x uniformly over [-2.5, 2.5], libm is 1.65× slower. That sweep is the wrong input distribution.

Capturing the arguments inference actually passes to erfc during a converge (508 848 calls):

|x| < 0.5     96.16%
0.5 - 0.84     2.05%
0.84 - 1.25    1.24%
1.25 - 2       0.54%
2 - 6          0.00%
> 6            0.00%

98% fall below 0.84375 — exactly the threshold under which FDLIBM's erfc is a rational polynomial with no exponential at all, while the Numerical Recipes form always evaluates one. On the captured trace libm is the faster of the two: 3.23 ns/call against 3.78.

So the hot-loop concern was real but pointed the wrong way, and it could only be settled by measuring the distribution rather than reasoning about instruction counts.

One thing worth recording as a caution: an ad-hoc Instant-based harness reported a 16% end-to-end speedup, which was an artifact of that harness allocating and leaking memory per run. Criterion's "no change" is the trustworthy number. The microbenchmark got the sign wrong in one direction and the ad-hoc harness got the magnitude wrong in the other; only the project's own benchmark and the captured trace agreed with each other.

What it bought

The issue predicted this would mostly tighten test tolerances. It did more:

  • compute_margin against exact standard-normal quantiles: 8.4e-8 → 1.7e-16. This was the residual left over after the erfc_inv sign fix in 683813e, and it is now at round-off. The draw margin is exact.
  • cdf(mu, mu, sigma) is now exactly 0.5.
  • sf + cdf sums to one within a ULP, from 3e-8 out.
  • erfcx's two branches now agree to round-off across the x = 2 crossover, where they were ~1e-7 apart. That matters beyond tolerances: the log-space evidence path (ln_sf, ln_interval) goes through the continued fraction while the linear path goes through erfc, and they were previously inconsistent with each other by more than either's own error.
  • Ten test tolerances tightened from 1e-6 to 1e-13..1e-15 — a 10^7 to 10^9 tightening of the regression net, which is the real durable win.
  • MIN_GRID_POINTS no longer sits on a cdf floor. Its documentation said refining past 8 192 nodes "buys nothing"; that was true only because of erfc. The error there is now purely the grid's own O(h²), so the constant is an accuracy/cost choice rather than a wall. It is unchanged, but for a different and now-correct reason.

Cost

One dependency: libm 0.2.16 — zero transitive dependencies, rust-lang maintained. statrs was ruled out as advised (its erf is only ~1e-10); spec_math was not needed once libm proved free.

Five numerical goldens moved, by 2.4e-9 to 6e-7 — the magnitude of the removed error. Re-recorded with more digits so future drift stays visible, and verified as movement toward truth rather than re-baselined.

Done in `3dd6593`. The answer to "is it worth paying for" turned out to be that there is nothing to pay. ## Accuracy `libm` (the Rust port of FDLIBM), measured against an independent incomplete-gamma reference over 20 000 points per range: | range | previous (NR `erfcc`) | libm | |---|---|---| | [-3, 0] | 7.95e-8 | **2.15e-14** | | [0, 0.5] | 8.69e-8 | **4.70e-14** | | [0.5, 2] | 9.38e-8 | **1.24e-12** | | [2, 6] | 1.04e-7 | **5.20e-14** | | [6, 26] | 1.07e-7 | **1.75e-13** | And the identities the issue flagged: - `erfc(0)` returns **exactly 1.0**, where it was `1.000_000_03`. - `abs(erfc(z) + erfc(-z) - 2)` is **2.22e-16** — one ULP — where it was 6.00e-8. ## Performance: no change detected `benches/batch.rs`, against a baseline taken immediately before the swap: ``` Batch::iteration change: [-3.3423% -0.2286% +2.2576%] (p = 0.89 > 0.05) No change in performance detected. ``` The acceptance criteria asked for the benchmark delta and for the change to be abandoned if it was not worth paying. There is no delta to weigh. ## Why it is free, and a correction to the issue's reasoning This issue assumed a 2-3× arithmetic cost. A naive microbenchmark agrees — sweeping `x` uniformly over [-2.5, 2.5], libm is **1.65× slower**. That sweep is the wrong input distribution. Capturing the arguments inference actually passes to `erfc` during a converge (508 848 calls): ``` |x| < 0.5 96.16% 0.5 - 0.84 2.05% 0.84 - 1.25 1.24% 1.25 - 2 0.54% 2 - 6 0.00% > 6 0.00% ``` **98% fall below 0.84375** — exactly the threshold under which FDLIBM's `erfc` is a rational polynomial with no exponential at all, while the Numerical Recipes form always evaluates one. On the captured trace libm is the *faster* of the two: 3.23 ns/call against 3.78. So the hot-loop concern was real but pointed the wrong way, and it could only be settled by measuring the distribution rather than reasoning about instruction counts. One thing worth recording as a caution: an ad-hoc `Instant`-based harness reported a **16% end-to-end speedup**, which was an artifact of that harness allocating and leaking memory per run. Criterion's "no change" is the trustworthy number. The microbenchmark got the sign wrong in one direction and the ad-hoc harness got the magnitude wrong in the other; only the project's own benchmark and the captured trace agreed with each other. ## What it bought The issue predicted this would mostly tighten test tolerances. It did more: - **`compute_margin` against exact standard-normal quantiles: 8.4e-8 → 1.7e-16.** This was the residual left over after the `erfc_inv` sign fix in `683813e`, and it is now at round-off. The draw margin is exact. - `cdf(mu, mu, sigma)` is now **exactly 0.5**. - `sf + cdf` sums to one within a ULP, from 3e-8 out. - **`erfcx`'s two branches now agree to round-off across the x = 2 crossover**, where they were ~1e-7 apart. That matters beyond tolerances: the log-space evidence path (`ln_sf`, `ln_interval`) goes through the continued fraction while the linear path goes through `erfc`, and they were previously inconsistent with each other by more than either's own error. - Ten test tolerances tightened from `1e-6` to `1e-13`..`1e-15` — a 10^7 to 10^9 tightening of the regression net, which is the real durable win. - `MIN_GRID_POINTS` no longer sits on a `cdf` floor. Its documentation said refining past 8 192 nodes "buys nothing"; that was true only because of `erfc`. The error there is now purely the grid's own O(h²), so the constant is an accuracy/cost choice rather than a wall. It is unchanged, but for a different and now-correct reason. ## Cost One dependency: `libm` 0.2.16 — zero transitive dependencies, rust-lang maintained. `statrs` was ruled out as advised (its `erf` is only ~1e-10); `spec_math` was not needed once `libm` proved free. Five numerical goldens moved, by 2.4e-9 to 6e-7 — the magnitude of the removed error. Re-recorded with more digits so future drift stays visible, and verified as movement toward truth rather than re-baselined.
Sign in to join this conversation.
No Branch/Tag Specified
Branches Tags
Labels
Clear labels
api
Shape or ergonomics of the public surface
 blocked
Waiting on another issue or an external decision
 breaking
Resolving it changes the public API; wants a major bump
 bug
Behaviour is wrong, including silently wrong
 decision
Needs a call to be made, not code to be written
 docs
Documentation, including doc comments and the README
 enhancement
New capability the crate does not have
 infra
CI, release process, packaging
 numerics
Floating-point accuracy, cancellation, convergence maths
 performance
Speed or memory
 refactor
Internal cleanup with no intended behaviour change
 testing
Test coverage and test infrastructure
No labels numerics performance
Milestone
No items
No Milestone
Projects
Clear projects
No projects
Assignees
Clear assignees
logaritmisk
No Assignees
1 Participants
Notifications
Due Date
No due date set.
Dependencies

No dependencies set.

Reference: logaritmisk/trueskill-tt#41