Geometric with tiny p: Cdf is 0 everywhere and Quantile/Median/Sample spin ~22 s to return int.MaxValue (Geometric(1e-17).Quantile(1e-9), true ≈ 1e8)

[matt-edmondson]

What's wrong

GeometricDistributionProvider (Essentials.DistributionProviders.Geometric/GeometricDistributionProvider.cs) works with 1 − p directly:

• ~L54–55: failureProbability = 1.0 - probability; logFailureProbability = Math.Log(failureProbability);
• ~L79: Cdf = 1.0 - Math.Pow(failureProbability, value + 1.0)
• ~L109: raw = Math.Ceiling(Math.Log(1.0 - probability) / logFailureProbability) - 1.0
• ~L121: the upward correction loop while (count < int.MaxValue && Cdf(count) < probability) count++;

For p ≲ 1.1e-16, 1.0 - p rounds to exactly 1:

• logFailureProbability becomes 0.
• Cdf is 0 and SurvivalFunction is 1 for every k.
• The quantile estimate is -inf, so count starts at 0.
• The correction loop, whose comment says it "can't run more than once", walks all the way to int.MaxValue, making about 2.1e9 Math.Pow calls.

For larger small p, Cdf still carries relative error of about eps/p. For example, Geometric(1e-12).Cdf(1e9) returns 9.99478e-4 against a true 9.995002e-4.

Repro

IDistribution<int> g = new GeometricDistributionProvider(1e-17);
g.Cdf(100_000_000);  // 0             true ≈ 1e-9
g.Quantile(1e-9);    // 2147483647    after ~21.9 s; true ≈ 1e8
_ = g.Median;        // ~22.8 s

The constructor accepts any p > 0, so nothing stops this.

Why it matters

Median and Sample (which goes through the inverse CDF) effectively hang a thread for tens of seconds and then return a wrong answer. The CDF is silently wrong for small success probabilities, which is a normal use of a geometric model (waiting for a rare event).

Suggested fix

• Store log1p(-p) instead of log(1 - p).
• Cdf = -expm1((k + 1) * log1p(-p)) and SurvivalFunction = exp((k + 1) * log1p(-p)).
• In Quantile, use log1p(-probability) / log1p(-p).
• Add accurate Log1P/ExpM1 helpers to SpecialFunctions, because .NET's double.LogP1/double.ExpM1 are implemented naively as Log(x + 1) / Exp(x) - 1.
• Bound the correction loops, for example to a few steps, falling back to bisection, so an estimate that is off can never turn into a 2e9-step walk.

Acceptance: Geometric(1e-17).Quantile(1e-9) returns ≈ 1e8 in microseconds, and Cdf matches -expm1((k+1)·log1p(-p)) to ~1e-15 relative for p down to 1e-300.

[matt-edmondson]

Triage

• Category: Bug
• Priority: Medium. With p ≲ 1e-16, Quantile, Median and Sample spin for about 22 s and return int.MaxValue, and Cdf is 0 everywhere. The constructor accepts such p, and larger small p still gets a noticeably wrong CDF.
• Area / suggested owner: Essentials distributions, GeometricDistributionProvider, plus SpecialFunctions
• Duplicates / in progress: None, and no open PR. It shares the ExpM1/Log1P helper work with Exponential Cdf/Quantile lose all precision near zero: Exponential(2).Cdf(1e-17) = 0 and Quantile(1e-17) = −0 (expected 2e-17 and 5e-18) #52 (Exponential near zero), so fix the two together.

Next step: Store log1p(-p), rewrite Cdf, SurvivalFunction and Quantile on it, and bound the correction loop (a few steps, then bisection) so a bad estimate can never become a 2e9-step walk.

---

Generated by Claude Code