Keep the exponential CDF and quantile precise near zero

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 and normal quantile, linked into the Normal, LogNormal, 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 and Poisson providers rather than placed in the interfaces-only package
 - `Shared/DistributionArguments.cs` - Parameter validation shared by every distribution provider, linked in the same way
 - `Essentials/ISerializationProvider.cs` - Object serialization/deserialization interface
 - `Essentials/ISerializationOptions.cs` - Configurable serialization options (naming, inclusion, boxing policies)

Essentials.DistributionProviders.Exponential/Essentials.DistributionProviders.Exponential.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.Exponential/ExponentialDistributionProvider.cs

@@ -1,4 +1,4 @@
-// Copyright (c) 2023-2026 ktsu-dev contributors
+// Copyright (c) 2023-2026 ktsu-dev contributors
 
 namespace ktsu.Essentials.DistributionProviders.Exponential;
 
@@ -50,7 +50,12 @@ public ExponentialDistributionProvider()
 	public double LogPdf(double value) => value < 0.0 ? double.NegativeInfinity : Math.Log(Rate) - (Rate * value);
 
 	/// <inheritdoc />
-	public double Cdf(double value) => value <= 0.0 ? 0.0 : 1.0 - Math.Exp(-Rate * value);
+	/// <remarks>
+	/// Evaluated as <c>-expm1(-λx)</c> rather than <c>1 - exp(-λx)</c>: the subtraction keeps only the
+	/// digits of a number near one, so short waits lost their precision and anything below about 1e-16
+	/// came out as zero.
+	/// </remarks>
+	public double Cdf(double value) => value <= 0.0 ? 0.0 : -SpecialFunctions.ExpM1(-Rate * value);
 
 	/// <summary>
 	/// Evaluates the survival function, the probability that the wait exceeds <paramref name="value"/>.
@@ -69,7 +74,9 @@ public double Quantile(double probability)
 	{
 		DistributionArguments.QuantileProbability(probability, nameof(probability));
 
-		return probability >= 1.0 ? double.PositiveInfinity : -Math.Log(1.0 - probability) / Rate;
+		// log1p rather than log(1 - p), for the same reason Cdf uses expm1: a small probability would
+		// otherwise round to a quantile of zero, and the negation of log(1) to one of negative zero.
+		return probability >= 1.0 ? double.PositiveInfinity : -SpecialFunctions.Log1P(-probability) / Rate;
 	}
 
 	/// <inheritdoc />

Essentials.Tests/DistributionProviderTests.cs

@@ -477,6 +477,29 @@ public void Exponential_Matches_Its_Closed_Form()
 		Assert.AreEqual(distribution.SurvivalFunction(1.0), distribution.SurvivalFunction(3.0) / distribution.SurvivalFunction(2.0), 1e-14);
 	}
 
+	[TestMethod]
+	public void Exponential_Keeps_Its_Precision_Near_Zero()
+	{
+		// 1 - exp(-λx) and -log(1 - p) keep only the digits of a number near one: below about 1e-16 both
+		// rounded to zero, and the quantile to negative zero, outside the support.
+		ExponentialDistributionProvider distribution = new(2.0);
+		Assert.AreEqual(2e-17, distribution.Cdf(1e-17), 2e-17 * 1e-15);
+		Assert.AreEqual(5e-18, distribution.Quantile(1e-17), 5e-18 * 1e-15);
+
+		ExponentialDistributionProvider unit = new();
+		Assert.AreEqual(9.9999999995e-11, unit.Cdf(1e-10), 9.9999999995e-11 * 1e-15);
+		Assert.AreEqual(1.00000000005e-10, unit.Quantile(1e-10), 1.00000000005e-10 * 1e-15);
+
+		Assert.IsFalse(double.IsNegative(distribution.Quantile(0.0)), "Quantile(0) is negative zero");
+		Assert.IsFalse(double.IsNegative(distribution.Quantile(1e-17)), "Quantile(1e-17) is negative");
+
+		// The rewrite must not disturb the ordinary range or the far tail.
+		Assert.AreEqual(0.6321205588285577, distribution.Cdf(0.5), 1e-15);
+		Assert.AreEqual(1.0, distribution.Cdf(1000.0), Tight);
+		Assert.AreEqual(1.0, distribution.Cdf(double.PositiveInfinity), Tight);
+		Assert.AreEqual(0.5, distribution.Quantile(distribution.Cdf(0.5)), 1e-14);
+	}
+
 	[TestMethod]
 	public void LogNormal_Matches_Its_Closed_Form()
 	{

Shared/SpecialFunctions.cs

@@ -1,4 +1,4 @@
-// Copyright (c) 2023-2026 ktsu-dev contributors
+// Copyright (c) 2023-2026 ktsu-dev contributors
 
 namespace ktsu.Essentials;
 
@@ -133,6 +133,63 @@ internal static double LogGamma(double value)
 	/// <returns>The natural logarithm of <paramref name="value"/> factorial.</returns>
 	internal static double LogFactorial(int value) => LogGamma(value + 1.0);
 
+	/// <summary>
+	/// Evaluates <c>exp(x) - 1</c> without losing precision for arguments near zero.
+	/// </summary>
+	/// <remarks>
+	/// Subtracting one from <c>exp(x)</c> keeps only the digits that survived the rounding of a number
+	/// near one, so the relative error grows like ε/|x| and the result is zero below about 1e-16. Kahan's
+	/// correction divides the rounded difference by the logarithm of the same rounded exponential, which
+	/// carries the identical rounding error and cancels it, leaving an answer good to a few ulps. Written
+	/// out because <c>double.ExpM1</c> is not available on every target and on .NET is the naive
+	/// subtraction anyway.
+	/// </remarks>
+	/// <param name="value">The argument.</param>
+	/// <returns><c>exp(value) - 1</c>.</returns>
+	internal static double ExpM1(double value)
+	{
+		double exponential = Math.Exp(value);
+		if (exponential == 1.0)
+		{
+			return value;
+		}
+
+		double difference = exponential - 1.0;
+		if (difference == -1.0 || double.IsPositiveInfinity(exponential))
+		{
+			return difference;
+		}
+
+		return difference * value / Math.Log(exponential);
+	}
+
+	/// <summary>
+	/// Evaluates <c>log(1 + x)</c> without losing precision for arguments near zero.
+	/// </summary>
+	/// <remarks>
+	/// The mirror of <see cref="ExpM1"/>: forming <c>1 + x</c> rounds away the low digits of a small
+	/// <c>x</c> before the logarithm sees them. Scaling by <c>x / ((1 + x) - 1)</c> — the ratio of the
+	/// argument to what actually survived the addition — restores them. Returns <paramref name="value"/>
+	/// itself, sign of zero included, when the addition rounds to exactly one.
+	/// </remarks>
+	/// <param name="value">The argument. Must be greater than or equal to -1.</param>
+	/// <returns><c>log(1 + value)</c>.</returns>
+	internal static double Log1P(double value)
+	{
+		double sum = 1.0 + value;
+		if (sum == 1.0)
+		{
+			return value;
+		}
+
+		if (double.IsPositiveInfinity(sum) || sum == 0.0)
+		{
+			return Math.Log(sum);
+		}
+
+		return Math.Log(sum) * value / (sum - 1.0);
+	}
+
 	/// <summary>
 	/// Evaluates the natural logarithm of a binomial coefficient.
 	/// </summary>