fix: compute Cauchy::{cdf,sf} tails as atan(scale / d) / pi

[FBruzzesi]

Closes #473.

@day01 thanks for suggesting I should give it a try. I asked Claude Opus 5.5 to point me to the scipy and the R implementations:

• SciPy uses arctan2(1, -z) / pi everywhere
• R's pcauchy uses the same split at |z| > 1 and the same atan(1 / |z|) tail, but forms z first, so it shares the tiny-scale zero.

I tried (and measured locally) these alternatives:

• atan2 everywhere: within one scale it is 1 ulp further from the true value on 11 of the existing exact assertions (mean error 0.31 to 0.39 ulp over 20k points with |z| <= 1). The split keeps those assertions as they are.

• atan2(scale, d) / pi in the tail: same accuracy in the far tail, and at most 0.3 ulp better near the switch point (worst case 1.65 vs 1.96 ulps).
• atan(scale / d) is 1.15x to 2.5x faster per call.

Comparison

Relative error against mpmath (50 digits, each input taken as the exact f64):

call	statrs 0.19.1	R pcauchy	SciPy 1.18.1	this PR
cdf(-100)	2.1e-14	4.8e-18	4.8e-18	4.8e-18
cdf(-1e12)	2.0e-4	5.3e-18	5.3e-18	5.3e-18
cdf(-1e16)	1.0	2.7e-17	2.7e-17	2.7e-17
cdf(-1e300)	1.0	1.1e-16	1.1e-16	1.1e-16
sf(1e300)	1.0	1.1e-16	1.1e-16	1.1e-16
cdf(-1e10), scale = 1e-300	1.0	1.0	1.0	4.4e-14

The last row's true value is subnormal (3.18e-311), where f64 spacing limits any result to about 7.8e-14 relative; this PR returns the closest f64 (4.4e-14 off). R and SciPy both divide by the scale first, so both return 0.0 there.

---

Note that if x - location itself overflows (both near ±f64::MAX), the result is still 0.0 or 1.0, the same as R and SciPy.

Summary by CodeRabbit

• Bug Fixes
  • Improved Cauchy distribution cumulative and survival probabilities for values far from the location, including extremely large distances and very small scales. Results are more reliable in the tails, where previous calculations could overflow or lose precision.

[coderabbitai]

No actionable comments were generated in the recent review. 🎉

ℹ️ Recent review info

⚙️ Run configuration

• Configuration used: defaults
• Review profile: CHILL
• Plan: Advanced
• Run ID: 627ac3b7-89f2-422c-8335-da5329b81420

📥 Commits

Reviewing files that changed from the base of the PR and between 55dfea2 and 22fe0b8.

📒 Files selected for processing: 1

• src/distribution/cauchy.rs

Included review availability: This review used your included allowance. Your plan provides up to 4 included reviews per hour; 3 remain after this review.

📝 Walkthrough

Cauchy::cdf and Cauchy::sf now use reciprocal-distance arctangent formulas in their respective far tails. Tests update survival-function expectations and check tail values at extreme distances, including subnormal results.

Changes: Cauchy tail calculations

Layer / File(s)	Summary
Far-tail formulas and validation src/distribution/cauchy.rs	cdf uses a reciprocal-distance arctangent when the lower-tail distance exceeds the scale. sf uses the corresponding formula when the upper-tail distance exceeds the scale. Tests update far-positive survival expectations and check both tails through distances of 1e300, including subnormal results.

Priority: ➖ Normal

Estimated code review effort: 2 (Simple) | ~10 minutes

Change: Bug fix · Severity of issue fixed: Medium

Suggested reviewers: day01

Merge Risk: ⚪ Minimal · up to 22fe0

This change improves the precision of Cauchy tail probabilities and adds tests for extreme and subnormal cases. No actionable merge-blocking risk was identified.

Pre-merge checks: 5

Passed checks: 5 passed

Check name	Status	Explanation
Description Check		Check skipped - CodeRabbit’s high-level summary is enabled.
Title check		The title clearly and concisely describes the main change: computing Cauchy CDF and survival-function tails with a reciprocal-argument arctangent formula.
Linked Issues check		Issue #473 requires accurate, nonzero Cauchy::cdf and Cauchy::sf tail results where f64 can represent them. The reviewed code uses (scale / distance).atan() / PI when the distance exceeds the …
Out of Scope Changes check		The changes stay within issue #473. They modify only Cauchy CDF and survival-function calculations and their tests. No unrelated distribution methods or unrelated behavior are changed.
Docstring Coverage		Docstring coverage is 100.00% which is sufficient. The required threshold is 80.00%. Docstring coverage is scoped to functions touched by this diff. Analyzed 5 functions across 1 files.

✨ Finishing Touches

🧪 Generate unit tests (beta)

• Create a new PR

• Autofix · Keep fixing CodeRabbit findings and required CI, and resolving merge conflicts

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

• X
• Mastodon
• Reddit
• LinkedIn

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

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.87%. Comparing base (55dfea2) to head (22fe0b8).

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #476   +/-   ##
=======================================
  Coverage   95.86%   95.87%           
=======================================
  Files          68       68           
  Lines       16675    16697   +22     
=======================================
+ Hits        15986    16008   +22     
  Misses        689      689           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[day01]

Thanks @FBruzzesi for the contribution :)

IMO it makes sense to address Cauchy as a whole here: #471, #472 and #473, including entropy and the sampling transform. Separate commits would keep this reviewable, with shared accuracy tests and benchmarks. This is a useful starting point, but there are still cases where intermediate calculations lose a representable result.

The goal for me is accuracy and performance. R/SciPy are useful comparisons, but they also have numerical problems on some of these inputs. I would aim for the closest result to a multiprecision reference and measure execution time before/after.

A few concrete things worth covering:

• Cauchy(1e308, 1e308).cdf(-1e308) still returns 0, while the reference is approximately 0.14758361765043326. The subtraction overflows before the new tail formula helps. The distance/scale calculation needs attention too.
• PDF/log-PDF still have bug: Cauchy::pdf returns 0.0 and ln_pdf returns -inf where the density is representable #472. Scaled ratios should help PDF; log-PDF needs its own calculation in log space, so an overflowing denominator or underflowed PDF does not destroy the result.
• Cauchy(0, 1e308).entropy() returns infinity instead of approximately 711.7272328891354. Avoiding 4*pi*scale before the logarithm addresses this case, but subnormal scales and cancellation near zero entropy also need tests.
• Inverse-CDF still has bug: Cauchy::inverse_cdf returns the same value for every quantile at or below 2^-55 #471, including cancellation between location and the scaled quantile. Additional precision needs to survive the final affine step. Sampling should share the improved transformation, with explicit tests for the RNG interval and endpoints.

There is also a small rounding failure in the current CDF:

let d = Cauchy::new(0.0, f64::from_bits(73)).unwrap();
let observed = d.cdf(-2.0);
let expected = f64::from_bits(12);

observed has bits 11, while the correctly rounded result has bits 12. The error against the mathematical reference is about 0.618311 ULP. The ratio loses precision in the subnormal range before division by pi.

Could we add multiprecision-derived fixtures using exact binary64 inputs and report maximum ULP error, correctly rounded counts, incorrect zeros and finite/infinite mismatches? Fixed-parameter monotonicity would also be useful. Please distinguish mathematical error in ULP from integer distance to the rounded reference. The current 1e-322 tolerance allows 20 subnormal ULP, which can hide failures like this. Correct rounding should be the target; a sampled grid alone cannot prove it globally.

For performance, could we commit reproducible Criterion benchmarks comparing the parent and this PR under the same conditions? Center, branch boundary, tails, subnormals and mixed inputs would give a useful picture. Keep setup outside timing, black-box inputs/results, and include the compiler, flags, CPU and repeated measurements with uncertainty. Check accuracy separately for std and libm.

My local release measurements were on an Apple M1 Max. They suggested slower mixed CDF/SF workloads, but repeated runs were variable, so I would not attach a fixed slowdown percentage to this. A before/after comparison on your machine is useful here. Comparing the isolated reciprocal-atan kernel with atan2 does not show the change in full API performance against the parent.

There may still be room to improve both. The ARM64 assembly I inspected repeats the distance subtraction in the central branch and has two divisions in the tail branch. Reusing the distance is worth measuring. Compensated reciprocal-pi multiplication and a small-angle specialization with an error bound are also candidates, but each needs both ULP and runtime results before calling it an improvement.

IMO one coherent Cauchy change with this evidence would be much stronger. I would like us to improve accuracy without a material regression in ordinary and mixed workloads.