diff --git a/Essentials.DistributionProviders.Uniform/UniformDistributionProvider.cs b/Essentials.DistributionProviders.Uniform/UniformDistributionProvider.cs index 32ff525..48079c2 100644 --- a/Essentials.DistributionProviders.Uniform/UniformDistributionProvider.cs +++ b/Essentials.DistributionProviders.Uniform/UniformDistributionProvider.cs @@ -17,6 +17,9 @@ public sealed class UniformDistributionProvider : IContinuousDistribution { private readonly double width; + /// Half the width, which stays finite when the full width overflows. + private readonly double halfWidth; + /// /// Initializes the standard uniform distribution on the unit interval. /// @@ -40,6 +43,7 @@ public UniformDistributionProvider(double minimum, double maximum) Minimum = minimum; Maximum = maximum; width = maximum - minimum; + halfWidth = (maximum / 2.0) - (minimum / 2.0); } /// @@ -49,13 +53,13 @@ public UniformDistributionProvider(double minimum, double maximum) public double Maximum { get; } /// - public double Mean => Minimum + (width / 2.0); + public double Mean => (Minimum / 2.0) + (Maximum / 2.0); /// public double Variance => width * width / 12.0; /// - 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; /// public double Cdf(double value) @@ -65,7 +69,14 @@ 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; } /// @@ -73,7 +84,9 @@ 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); } /// diff --git a/Essentials.Tests/RandomRangeBoundaryTests.cs b/Essentials.Tests/RandomRangeBoundaryTests.cs new file mode 100644 index 0000000..039d1d1 --- /dev/null +++ b/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; + +/// +/// Pins and the uniform distribution at the edges +/// a random sample is unlikely to reach: the top unit draw, and finite bounds whose width overflows. +/// +[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})"); + } + + /// 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..bc815e4 100644 --- a/Essentials/IRandomProvider.cs +++ b/Essentials/IRandomProvider.cs @@ -185,7 +185,17 @@ public double NextDouble(double minInclusive, double maxExclusive) 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; } /// diff --git a/Essentials/RandomHelpers.cs b/Essentials/RandomHelpers.cs index 89711f6..3e2b048 100644 --- a/Essentials/RandomHelpers.cs +++ b/Essentials/RandomHelpers.cs @@ -17,6 +17,30 @@ internal static class RandomHelpers /// The step between adjacent 24-bit fixed-point values on the unit interval. internal const float UnitScaleSingle = 1.0f / 16777216.0f; + /// + /// Returns the largest double below a finite . + /// + /// + /// Math.BitDecrement, 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. + /// + /// A finite value. + /// The adjacent representable value toward negative infinity. + 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); + } + /// /// Draws a 32-bit value uniformly from [0, ) without modulo bias. ///