Keep NextDoubleExclusive strictly below 1 at the top draw [patch]

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<double> normal = new NormalDistributionProvider();
+		IDistribution<double> 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<ulong, IRandomProvider> factory, string name)
 	{
 		ulong[] first = Draw(factory(99UL));
@@ -446,4 +468,17 @@ private static void AssertSeedBehaviour(Func<ulong, IRandomProvider> factory, st
 	}
 
 	private static ulong[] Draw(IRandomProvider random) => [.. Enumerable.Range(0, 32).Select(_ => random.NextUInt64())];
+
+	/// <summary>A provider that fills every buffer with the same 64-bit word, for pinning boundary draws.</summary>
+	private sealed class ConstantRandomProvider(ulong word) : IRandomProvider
+	{
+		public void NextBytes(Span<byte> destination)
+		{
+			byte[] bytes = BitConverter.GetBytes(word);
+			for (int i = 0; i < destination.Length; i++)
+			{
+				destination[i] = bytes[i % bytes.Length];
+			}
+		}
+	}
 }

Essentials/IRandomProvider.cs

@@ -158,7 +158,7 @@
 	/// can carry on [0, 1).
 	/// </remarks>
 	/// <returns>A value in the range [0, 1).</returns>
 	public double NextDouble() => (NextUInt64() >> 11) * RandomHelpers.UnitScale;
 
 	/// <summary>
 	/// Returns a uniformly distributed double in the open unit interval, excluding both endpoints.
@@ -166,11 +166,12 @@
 	/// <remarks>
 	/// 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.
 	/// </remarks>
 	/// <returns>A value in the range (0, 1).</returns>
-	public double NextDoubleExclusive() => ((NextUInt64() >> 11) + 0.5) * RandomHelpers.UnitScale;
+	public double NextDoubleExclusive() => ((NextUInt64() >> 12) + 0.5) * RandomHelpers.OpenUnitScale;
 
 	/// <summary>
 	/// Returns a uniformly distributed double within the specified range.

Essentials/RandomHelpers.cs

@@ -14,6 +14,9 @@ internal static class RandomHelpers
 	/// <summary>The step between adjacent 53-bit fixed-point values on the unit interval.</summary>
 	internal const double UnitScale = 1.0 / 9007199254740992.0;
 
+	/// <summary>The step between adjacent 52-bit fixed-point values, used for the open unit interval.</summary>
+	internal const double OpenUnitScale = 1.0 / 4503599627370496.0;
+
 	/// <summary>The step between adjacent 24-bit fixed-point values on the unit interval.</summary>
 	internal const float UnitScaleSingle = 1.0f / 16777216.0f;