diff --git a/Essentials.Tests/RandomProviderTests.cs b/Essentials.Tests/RandomProviderTests.cs index cd0c873..fa77fe1 100644 --- a/Essentials.Tests/RandomProviderTests.cs +++ b/Essentials.Tests/RandomProviderTests.cs @@ -3,6 +3,8 @@ namespace ktsu.Essentials.Tests; using ktsu.Essentials; +using ktsu.Essentials.DistributionProviders.LogNormal; +using ktsu.Essentials.DistributionProviders.Normal; using ktsu.Essentials.RandomProviders.Crypto; using ktsu.Essentials.RandomProviders.Native; using ktsu.Essentials.RandomProviders.Pcg; @@ -435,6 +437,26 @@ public void Word_Sized_Draws_Agree_With_The_Byte_Buffer() } } + [TestMethod] + [DataRow(ulong.MaxValue)] + [DataRow(0UL)] + public void NextDoubleExclusive_Excludes_Both_Ends_At_The_Extreme_Draws(ulong word) + { + // The top 53-bit draw plus half a step is a tie between 2^53 - 1 and 2^53, which rounds to + // even and so returned exactly 1.0; the quantile of an unbounded distribution is infinite there. + IRandomProvider random = new ConstantRandomProvider(word); + + double u = random.NextDoubleExclusive(); + Assert.IsTrue(u is > 0.0 and < 1.0, $"NextDoubleExclusive returned {u:R} for 0x{word:X16}"); + + IDistribution normal = new NormalDistributionProvider(); + IDistribution logNormal = new LogNormalDistributionProvider(); + double normalSample = normal.Sample(random); + double logNormalSample = logNormal.Sample(random); + Assert.IsTrue(double.IsFinite(normalSample), $"Normal(0, 1) sampled {normalSample} for 0x{word:X16}"); + Assert.IsTrue(double.IsFinite(logNormalSample), $"LogNormal(0, 1) sampled {logNormalSample} for 0x{word:X16}"); + } + private static void AssertSeedBehaviour(Func factory, string name) { ulong[] first = Draw(factory(99UL)); @@ -446,4 +468,17 @@ private static void AssertSeedBehaviour(Func factory, st } private static ulong[] Draw(IRandomProvider random) => [.. Enumerable.Range(0, 32).Select(_ => random.NextUInt64())]; + + /// A provider that fills every buffer with the same 64-bit word, for pinning boundary draws. + private sealed class ConstantRandomProvider(ulong word) : IRandomProvider + { + public void NextBytes(Span destination) + { + byte[] bytes = BitConverter.GetBytes(word); + for (int i = 0; i < destination.Length; i++) + { + destination[i] = bytes[i % bytes.Length]; + } + } + } } diff --git a/Essentials/IRandomProvider.cs b/Essentials/IRandomProvider.cs index 41b44c8..d190c4e 100644 --- a/Essentials/IRandomProvider.cs +++ b/Essentials/IRandomProvider.cs @@ -166,11 +166,12 @@ public long NextInt64(long minInclusive, long maxExclusive) /// /// Inverse-transform sampling needs this: the quantile function of an unbounded distribution is /// infinite at 0 and at 1, so a draw that can return either endpoint would produce an infinity - /// rather than a sample. The 53-bit draw is offset by half a step, which excludes both ends - /// without distorting the spacing of the values in between. + /// rather than a sample. A 52-bit draw is offset by half a step, which excludes both ends + /// without distorting the spacing of the values in between. It is 52 bits rather than 53 because + /// the half-step sum of the largest 53-bit draw is not representable and rounds up to exactly 1. /// /// A value in the range (0, 1). - public double NextDoubleExclusive() => ((NextUInt64() >> 11) + 0.5) * RandomHelpers.UnitScale; + public double NextDoubleExclusive() => ((NextUInt64() >> 12) + 0.5) * RandomHelpers.OpenUnitScale; /// /// Returns a uniformly distributed double within the specified range. diff --git a/Essentials/RandomHelpers.cs b/Essentials/RandomHelpers.cs index 89711f6..a343735 100644 --- a/Essentials/RandomHelpers.cs +++ b/Essentials/RandomHelpers.cs @@ -14,6 +14,9 @@ internal static class RandomHelpers /// The step between adjacent 53-bit fixed-point values on the unit interval. internal const double UnitScale = 1.0 / 9007199254740992.0; + /// The step between adjacent 52-bit fixed-point values, used for the open unit interval. + internal const double OpenUnitScale = 1.0 / 4503599627370496.0; + /// The step between adjacent 24-bit fixed-point values on the unit interval. internal const float UnitScaleSingle = 1.0f / 16777216.0f;