Keep the geometric CDF and quantile precise for a rare success

CLAUDE.md

@@ -45,7 +45,7 @@ This is a .NET library (`ktsu.Essentials`) providing high-performance interfaces
 - `Essentials/IContinuousDistribution.cs` - `IDistribution<double>` plus `Pdf`/`LogPdf`
 - `Essentials/IDiscreteDistribution.cs` - `IDistribution<int>` plus `Pmf`/`LogPmf`, and a default `Quantile` that bisects the CDF. That default is an explicit reimplementation of `IDistribution<int>.Quantile`, which CA1033 misreads as an explicit implementation on a class; the suppression on it says why
 - `Shared/DeflateBound.cs` - zlib's `deflateBound` plus container overhead, linked into the Deflate, Gzip and ZLib providers as their `GetMaxCompressedLength`
-- `Shared/SpecialFunctions.cs` - Incomplete gamma, incomplete beta, error function, normal quantile and cancellation-free `ExpM1`/`Log1P`, linked into the Normal, LogNormal, Exponential, Binomial and Poisson providers rather than placed in the interfaces-only package
+- `Shared/SpecialFunctions.cs` - Incomplete gamma, incomplete beta, error function, normal quantile and cancellation-free `ExpM1`/`Log1P`, linked into the Normal, LogNormal, Exponential, Binomial, Poisson and Geometric providers rather than placed in the interfaces-only package
 - `Shared/DistributionArguments.cs` - Parameter validation shared by every distribution provider, linked in the same way
 - `Shared/PersistenceDefaults.cs` - The `IPersistenceProvider` default `RetrieveOrCreateAsync` body, linked into the FileSystem and Temp providers, which declare that member themselves and so hide the default
 - `Essentials/ISerializationProvider.cs` - Object serialization/deserialization interface

Essentials.DistributionProviders.Geometric/Essentials.DistributionProviders.Geometric.csproj

@@ -13,6 +13,7 @@
   </ItemGroup>
 
   <ItemGroup>
+    <Compile Include="..\Shared\SpecialFunctions.cs" Link="SpecialFunctions.cs" />
     <Compile Include="..\Shared\DistributionArguments.cs" Link="DistributionArguments.cs" />
   </ItemGroup>
 

Essentials.DistributionProviders.Geometric/GeometricDistributionProvider.cs

@@ -1,4 +1,4 @@
-// Copyright (c) 2023-2026 ktsu-dev contributors
+// Copyright (c) 2023-2026 ktsu-dev contributors
 
 namespace ktsu.Essentials.DistributionProviders.Geometric;
 
@@ -23,6 +23,11 @@ namespace ktsu.Essentials.DistributionProviders.Geometric;
 /// </remarks>
 public sealed class GeometricDistributionProvider : IDiscreteDistribution
 {
+	/// <summary>
+	/// How far either way the closed-form quantile may be walked before giving up on it and bisecting.
+	/// </summary>
+	private const int MaxQuantileCorrection = 2;
+
 	private readonly double failureProbability;
 	private readonly double logFailureProbability;
 
@@ -52,7 +57,11 @@ public GeometricDistributionProvider(double probability)
 
 		Probability = probability;
 		failureProbability = 1.0 - probability;
-		logFailureProbability = Math.Log(failureProbability);
+
+		// Taken as log1p(-p) rather than the log of 1 - p, which rounds to exactly one below p of about
+		// 1e-16 and leaves a logarithm of zero: every CDF value would then be zero and every quantile a
+		// walk to int.MaxValue. Above that it still costs digits in proportion to how small p is.
+		logFailureProbability = SpecialFunctions.Log1P(-probability);
 	}
 
 	/// <summary>
@@ -73,10 +82,29 @@ public GeometricDistributionProvider(double probability)
 	public double Variance => failureProbability / (Probability * Probability);
 
 	/// <inheritdoc />
-	public double Pmf(int value) => value < 0 ? 0.0 : Math.Pow(failureProbability, value) * Probability;
+	public double Pmf(int value)
+	{
+		if (value < 0)
+		{
+			return 0.0;
+		}
+
+		// Zero failures is the success probability itself. Taking it through the logarithm would multiply a
+		// zero count by the infinite logarithm of a certain success, which is a NaN rather than one.
+		if (value == 0)
+		{
+			return Probability;
+		}
+
+		return Math.Exp(value * logFailureProbability) * Probability;
+	}
 
 	/// <inheritdoc />
-	public double Cdf(int value) => value < 0 ? 0.0 : 1.0 - Math.Pow(failureProbability, value + 1.0);
+	/// <remarks>
+	/// One minus <c>(1 - p)^(k + 1)</c>, written as <c>-expm1((k + 1) log1p(-p))</c> so that neither the
+	/// failure probability nor the subtraction from one throws away the digits a small CDF lives in.
+	/// </remarks>
+	public double Cdf(int value) => value < 0 ? 0.0 : -SpecialFunctions.ExpM1((value + 1.0) * logFailureProbability);
 
 	/// <summary>
 	/// Evaluates the survival function, the probability of more than <paramref name="value"/> failures.
@@ -87,7 +115,7 @@ public GeometricDistributionProvider(double probability)
 	/// </remarks>
 	/// <param name="value">The failure count to evaluate at.</param>
 	/// <returns>A probability in the range [0, 1].</returns>
-	public double SurvivalFunction(int value) => value < 0 ? 1.0 : Math.Pow(failureProbability, value + 1.0);
+	public double SurvivalFunction(int value) => value < 0 ? 1.0 : Math.Exp((value + 1.0) * logFailureProbability);
 
 	/// <inheritdoc />
 	public int Quantile(double probability)
@@ -106,23 +134,51 @@ public int Quantile(double probability)
 
 		// The smallest count k whose CDF reaches the target satisfies (k + 1) >= log(1 - p) / log(1 - q),
 		// with the inequality flipping because the logarithm of the failure probability is negative.
-		double raw = Math.Ceiling(Math.Log(1.0 - probability) / logFailureProbability) - 1.0;
+		double raw = Math.Ceiling(SpecialFunctions.Log1P(-probability) / logFailureProbability) - 1.0;
 		int count = raw <= 0.0 ? 0 : raw >= int.MaxValue ? int.MaxValue : (int)raw;
 
 		// That form is exact in real arithmetic, but where the target is itself a CDF value the ratio of
 		// the two logarithms lands a rounding step either side of a whole number, and the ceiling turns
-		// a step into a whole outcome. One comparison against the CDF in each direction puts it back, and
-		// neither loop can run more than once because the estimate is never off by more than one.
-		while (count > 0 && Cdf(count - 1) >= probability)
+		// a step into a whole outcome. A comparison against the CDF in each direction puts it back. The
+		// estimate should never be off by more than one, but the walks are bounded all the same: an
+		// unbounded one is what turned a bad estimate into two billion steps, and bisecting the whole
+		// support costs only 31 CDF evaluations.
+		for (int step = 0; step < MaxQuantileCorrection && count > 0 && Cdf(count - 1) >= probability; step++)
 		{
 			count--;
 		}
 
-		while (count < int.MaxValue && Cdf(count) < probability)
+		for (int step = 0; step < MaxQuantileCorrection && count < int.MaxValue && Cdf(count) < probability; step++)
 		{
 			count++;
 		}
 
-		return count;
+		bool settled = (count == int.MaxValue || Cdf(count) >= probability) && (count == 0 || Cdf(count - 1) < probability);
+		return settled ? count : Bisect(probability);
+	}
+
+	/// <summary>
+	/// Finds the smallest count whose CDF reaches <paramref name="probability"/> by bisecting the support.
+	/// </summary>
+	/// <param name="probability">The target probability.</param>
+	/// <returns>The quantile.</returns>
+	internal int Bisect(double probability)
+	{
+		int low = 0;
+		int high = int.MaxValue;
+		while (low < high)
+		{
+			int mid = low + ((high - low) / 2);
+			if (Cdf(mid) >= probability)
+			{
+				high = mid;
+			}
+			else
+			{
+				low = mid + 1;
+			}
+		}
+
+		return low;
 	}
 }

Essentials.Tests/DistributionProviderTests.cs

@@ -1,4 +1,4 @@
-// Copyright (c) 2023-2026 ktsu-dev contributors
+// Copyright (c) 2023-2026 ktsu-dev contributors
 
 namespace ktsu.Essentials.Tests;
 
@@ -672,6 +672,44 @@ public void Geometric_Matches_Its_Closed_Form()
 		Assert.AreEqual(0, certain.Quantile(1.0));
 	}
 
+	[TestMethod]
+	public void Geometric_Keeps_Its_Precision_For_A_Rare_Success()
+	{
+		// Below p of about 1e-16, 1 - p rounds to exactly one, so a CDF built on it is zero everywhere and
+		// the quantile walks one outcome at a time to int.MaxValue. The reference values were computed to
+		// 60 digits outside this codebase.
+		GeometricDistributionProvider rare = new(1e-17);
+		Assert.AreEqual(1.0000000095e-9, rare.Cdf(100_000_000), 1.0000000095e-9 * 1e-12);
+		Assert.AreEqual(1.0 - 1.0000000095e-9, rare.SurvivalFunction(100_000_000), 1e-15);
+		Assert.AreEqual(100_000_000, rare.Quantile(1e-9));
+
+		// Above that threshold 1 - p survives but carries a relative error of about ε/p into the CDF.
+		GeometricDistributionProvider uncommon = new(1e-12);
+		Assert.AreEqual(9.995001676245083e-4, uncommon.Cdf(1_000_000_000), 9.995001676245083e-4 * 1e-12);
+
+		// At the bottom of the double range the CDF is (k + 1) p to every digit that can be stored.
+		GeometricDistributionProvider vanishing = new(1e-300);
+		Assert.AreEqual(1001e-300, vanishing.Cdf(1000), 1001e-300 * 1e-15);
+		Assert.AreEqual(1e-300, vanishing.Pmf(1000), 1e-300 * 1e-15);
+		Assert.AreEqual(0.0, vanishing.Pmf(-1), Tight);
+	}
+
+	[TestMethod]
+	public void Geometric_Bisection_Fallback_Agrees_With_The_Closed_Form_Quantile()
+	{
+		// The quantile falls back to bisecting the support only when its closed-form estimate is off by more
+		// than the bounded correction allows, which no ordinary parameter provokes. Checked directly, so the
+		// fallback is known to land on the same outcome whenever it is needed.
+		foreach (double p in new[] { 0.25, 1e-6, 1e-17 })
+		{
+			GeometricDistributionProvider distribution = new(p);
+			foreach (double q in new[] { 0.0, 1e-9, 0.25, 0.5, 0.999 })
+			{
+				Assert.AreEqual(distribution.Quantile(q), distribution.Bisect(q), $"Geometric({p}) at {q}");
+			}
+		}
+	}
+
 	[TestMethod]
 	public void Categorical_Normalises_Its_Weights()
 	{