Keep NextDouble(min, max) below max and finite for wide bounds [patch]

Essentials.DistributionProviders.Uniform/UniformDistributionProvider.cs

@@ -17,6 +17,9 @@ public sealed class UniformDistributionProvider : IContinuousDistribution
 {
 	private readonly double width;
 
+	/// <summary>Half the width, which stays finite when the full width overflows.</summary>
+	private readonly double halfWidth;
+
 	/// <summary>
 	/// Initializes the standard uniform distribution on the unit interval.
 	/// </summary>
@@ -40,6 +43,7 @@ public UniformDistributionProvider(double minimum, double maximum)
 		Minimum = minimum;
 		Maximum = maximum;
 		width = maximum - minimum;
+		halfWidth = (maximum / 2.0) - (minimum / 2.0);
 	}
 
 	/// <inheritdoc />
@@ -49,13 +53,13 @@ public UniformDistributionProvider(double minimum, double maximum)
 	public double Maximum { get; }
 
 	/// <inheritdoc />
-	public double Mean => Minimum + (width / 2.0);
+	public double Mean => (Minimum / 2.0) + (Maximum / 2.0);
 
 	/// <inheritdoc />
 	public double Variance => width * width / 12.0;
 
 	/// <inheritdoc />
-	public double Pdf(double value) => value >= Minimum && value <= Maximum ? 1.0 / width : 0.0;
+	public double Pdf(double value) => value >= Minimum && value <= Maximum ? 0.5 / halfWidth : 0.0;
 
 	/// <inheritdoc />
 	public double Cdf(double value)
@@ -65,15 +69,24 @@ public double Cdf(double value)
 			return 0.0;
 		}
 
-		return value >= Maximum ? 1.0 : (value - Minimum) / width;
+		if (value >= Maximum)
+		{
+			return 1.0;
+		}
+
+		// Bounds of opposite sign can be far enough apart that the width overflows; halving both
+		// sides of the ratio keeps every term finite.
+		return double.IsInfinity(width) ? ((value / 2.0) - (Minimum / 2.0)) / halfWidth : (value - Minimum) / width;
 	}
 
 	/// <inheritdoc />
 	public double Quantile(double probability)
 	{
 		DistributionArguments.QuantileProbability(probability, nameof(probability));
 
-		return Minimum + (probability * width);
+		return double.IsInfinity(width)
+			? (Minimum * (1.0 - probability)) + (Maximum * probability)
+			: Minimum + (probability * width);
 	}
 
 	/// <inheritdoc />

Essentials.Tests/RandomRangeBoundaryTests.cs

@@ -0,0 +1,89 @@
+// Copyright (c) 2023-2026 ktsu-dev contributors
+
+namespace ktsu.Essentials.Tests;
+
+using ktsu.Essentials;
+using ktsu.Essentials.DistributionProviders.Uniform;
+using Microsoft.VisualStudio.TestTools.UnitTesting;
+using System;
+
+/// <summary>
+/// Pins <see cref="IRandomProvider.NextDouble(double, double)"/> and the uniform distribution at the edges
+/// a random sample is unlikely to reach: the top unit draw, and finite bounds whose width overflows.
+/// </summary>
+[TestClass]
+public class RandomRangeBoundaryTests
+{
+	[TestMethod]
+	[DataRow(1.0, 1.0000000000000002)]
+	[DataRow(1.0, 2.0)]
+	[DataRow(0.3, 0.7)]
+	[DataRow(5.0, 5.5)]
+	[DataRow(0.9, 1.0)]
+	public void NextDouble_Never_Returns_The_Exclusive_Bound_At_The_Top_Draw(double minInclusive, double maxExclusive)
+	{
+		IRandomProvider random = new ConstantRandomProvider(ulong.MaxValue);
+
+		double value = random.NextDouble(minInclusive, maxExclusive);
+
+		Assert.IsTrue(value >= minInclusive && value < maxExclusive, $"NextDouble({minInclusive:R}, {maxExclusive:R}) returned {value:R}");
+	}
+
+	[TestMethod]
+	[DataRow(ulong.MaxValue)]
+	[DataRow(0UL)]
+	[DataRow(0x8000000000000000UL)]
+	public void NextDouble_Stays_Finite_When_The_Width_Overflows(ulong word)
+	{
+		IRandomProvider random = new ConstantRandomProvider(word);
+
+		double value = random.NextDouble(double.MinValue, double.MaxValue);
+
+		Assert.IsTrue(double.IsFinite(value), $"NextDouble(MinValue, MaxValue) returned {value} for 0x{word:X16}");
+		Assert.IsTrue(value < double.MaxValue, $"NextDouble(MinValue, MaxValue) returned the exclusive bound for 0x{word:X16}");
+	}
+
+	[TestMethod]
+	public void Uniform_Handles_Bounds_Whose_Width_Overflows()
+	{
+		IContinuousDistribution uniform = new UniformDistributionProvider(-1e308, 1e308);
+
+		Assert.AreEqual(0.0, uniform.Mean);
+		Assert.AreEqual(0.0, uniform.Median);
+		Assert.AreEqual(0.5, uniform.Cdf(0.0));
+		Assert.AreEqual(0.75, uniform.Cdf(1e308 / 2.0), 1e-12);
+		Assert.AreEqual(0.0, uniform.Quantile(0.5));
+		Assert.AreEqual(-1e308, uniform.Quantile(0.0));
+		Assert.AreEqual(1e308, uniform.Quantile(1.0));
+		Assert.IsTrue(uniform.Pdf(0.0) > 0.0, "Pdf underflowed to 0 inside the support");
+		Assert.IsTrue(double.IsFinite(uniform.Sample(new ConstantRandomProvider(ulong.MaxValue))));
+		Assert.IsTrue(double.IsFinite(uniform.Sample(new ConstantRandomProvider(0UL))));
+	}
+
+	[TestMethod]
+	[DataRow(1.0)]
+	[DataRow(-1.0)]
+	[DataRow(0.0)]
+	[DataRow(-0.0)]
+	[DataRow(double.Epsilon)]
+	[DataRow(-double.Epsilon)]
+	[DataRow(double.MaxValue)]
+	[DataRow(double.MinValue + 1e292)]
+	public void NextDown_Matches_Math_BitDecrement(double value)
+	{
+		Assert.AreEqual(Math.BitDecrement(value), RandomHelpers.NextDown(value), $"NextDown({value:R})");
+	}
+
+	/// <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.
@@ -185,7 +185,17 @@
 		RandomHelpers.EnsureFinite(maxExclusive, nameof(maxExclusive));
 		Ensure.NotLessThan(maxExclusive, minInclusive);
 
-		return minInclusive + (NextDouble() * (maxExclusive - minInclusive));
+		double unit = NextDouble();
+		double width = maxExclusive - minInclusive;
+
+		// Bounds of opposite sign can be far enough apart that the width overflows, so interpolate
+		// between them instead; each term is then no larger in magnitude than its bound.
+		double value = double.IsInfinity(width)
+			? (minInclusive * (1.0 - unit)) + (maxExclusive * unit)
+			: minInclusive + (unit * width);
+
+		// Rounding the final addition can land on the excluded bound when the unit draw is close to 1.
+		return value >= maxExclusive && maxExclusive > minInclusive ? RandomHelpers.NextDown(maxExclusive) : value;
 	}
 
 	/// <summary>

Essentials/RandomHelpers.cs

@@ -17,6 +17,30 @@ internal static class RandomHelpers
 	/// <summary>The step between adjacent 24-bit fixed-point values on the unit interval.</summary>
 	internal const float UnitScaleSingle = 1.0f / 16777216.0f;
 
+	/// <summary>
+	/// Returns the largest double below a finite <paramref name="value"/>.
+	/// </summary>
+	/// <remarks>
+	/// <c>Math.BitDecrement</c>, which netstandard2.1 does not have. For a finite double the ordering of
+	/// the values follows the ordering of their bit patterns within each sign, so a step down is one
+	/// less in magnitude above zero and one more below it.
+	/// </remarks>
+	/// <param name="value">A finite value.</param>
+	/// <returns>The adjacent representable value toward negative infinity.</returns>
+	internal static double NextDown(double value)
+	{
+		long bits = BitConverter.DoubleToInt64Bits(value);
+
+		// Both zeros, tested on the bits so that -0.0 is caught too: the step down from either is the
+		// smallest negative subnormal.
+		if ((bits & long.MaxValue) == 0)
+		{
+			return -double.Epsilon;
+		}
+
+		return BitConverter.Int64BitsToDouble(bits > 0 ? bits - 1 : bits + 1);
+	}
+
 	/// <summary>
 	/// Draws a 32-bit value uniformly from [0, <paramref name="range"/>) without modulo bias.
 	/// </summary>