diff --git a/CLAUDE.md b/CLAUDE.md index 38377ed..2382503 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -45,7 +45,7 @@ This is a .NET library (`ktsu.Essentials`) providing high-performance interfaces - `Essentials/IContinuousDistribution.cs` - `IDistribution` plus `Pdf`/`LogPdf` - `Essentials/IDiscreteDistribution.cs` - `IDistribution` plus `Pmf`/`LogPmf`, and a default `Quantile` that bisects the CDF. That default is an explicit reimplementation of `IDistribution.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) diff --git a/Essentials.DistributionProviders.Exponential/Essentials.DistributionProviders.Exponential.csproj b/Essentials.DistributionProviders.Exponential/Essentials.DistributionProviders.Exponential.csproj index 1be4173..cdeaed3 100644 --- a/Essentials.DistributionProviders.Exponential/Essentials.DistributionProviders.Exponential.csproj +++ b/Essentials.DistributionProviders.Exponential/Essentials.DistributionProviders.Exponential.csproj @@ -13,6 +13,7 @@ + diff --git a/Essentials.DistributionProviders.Exponential/ExponentialDistributionProvider.cs b/Essentials.DistributionProviders.Exponential/ExponentialDistributionProvider.cs index 763b47b..7dba9cf 100644 --- a/Essentials.DistributionProviders.Exponential/ExponentialDistributionProvider.cs +++ b/Essentials.DistributionProviders.Exponential/ExponentialDistributionProvider.cs @@ -50,7 +50,12 @@ public ExponentialDistributionProvider() public double LogPdf(double value) => value < 0.0 ? double.NegativeInfinity : Math.Log(Rate) - (Rate * value); /// - public double Cdf(double value) => value <= 0.0 ? 0.0 : 1.0 - Math.Exp(-Rate * value); + /// + /// Evaluated as -expm1(-λx) rather than 1 - exp(-λx): 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. + /// + public double Cdf(double value) => value <= 0.0 ? 0.0 : -SpecialFunctions.ExpM1(-Rate * value); /// /// Evaluates the survival function, the probability that the wait exceeds . @@ -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; } /// diff --git a/Essentials.Tests/DistributionProviderTests.cs b/Essentials.Tests/DistributionProviderTests.cs index ec6445a..f365e22 100644 --- a/Essentials.Tests/DistributionProviderTests.cs +++ b/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() { diff --git a/Shared/SpecialFunctions.cs b/Shared/SpecialFunctions.cs index 1425b41..6b500ac 100644 --- a/Shared/SpecialFunctions.cs +++ b/Shared/SpecialFunctions.cs @@ -133,6 +133,70 @@ internal static double LogGamma(double value) /// The natural logarithm of factorial. internal static double LogFactorial(int value) => LogGamma(value + 1.0); + /// + /// Evaluates exp(x) - 1 without losing precision for arguments near zero. + /// + /// + /// Subtracting one from exp(x) 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 double.ExpM1 is not available on every target and on .NET is the naive + /// subtraction anyway. + /// + /// The argument. + /// exp(value) - 1. + internal static double ExpM1(double value) + { + // Below ε the answer is the argument itself to within an ulp, and exp rounds to exactly one, which + // would leave the correction below dividing zero by zero. Above ε it never rounds to one. + if (Math.Abs(value) <= Epsilon) + { + return value; + } + + double exponential = Math.Exp(value); + if (exponential <= 0.0) + { + return -1.0; + } + + if (double.IsPositiveInfinity(exponential)) + { + return exponential; + } + + return (exponential - 1.0) * value / Math.Log(exponential); + } + + /// + /// Evaluates log(1 + x) without losing precision for arguments near zero. + /// + /// + /// The mirror of : forming 1 + x rounds away the low digits of a small + /// x before the logarithm sees them. Scaling by x / ((1 + x) - 1) — the ratio of the + /// argument to what actually survived the addition — restores them. Returns + /// itself, sign of zero included, when it is too small for the addition to register. + /// + /// The argument. Must be greater than or equal to -1. + /// log(1 + value). + internal static double Log1P(double value) + { + // As in ExpM1: below ε the answer is the argument, and 1 + x may round to exactly one. + if (Math.Abs(value) <= Epsilon) + { + return value; + } + + double sum = 1.0 + value; + if (double.IsPositiveInfinity(sum) || sum <= 0.0) + { + return Math.Log(sum); + } + + return Math.Log(sum) * value / (sum - 1.0); + } + /// /// Evaluates the natural logarithm of a binomial coefficient. ///