Random: RandomSource.DoSampleInteger for direct integer sampling without floating points #359

MersenneTwister integer sampling now matches standard implementation.

src/Numerics/Random/CryptoRandomSource.cs

@@ -86,18 +86,32 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             var bytes = new byte[4];
             _crypto.GetBytes(bytes);
             return BitConverter.ToUInt32(bytes, 0)*Reciprocal;
         }
 
+        /// <summary>
+        /// Returns a random 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>
+        /// </summary>
+        protected sealed override int DoSampleInteger()
+        {
+            var bytes = new byte[4];
+            _crypto.GetBytes(bytes);
+            uint uint32 = BitConverter.ToUInt32(bytes, 0);
+            int int31 = (int)uint32 >> 1;
+            if (int31 == int.MaxValue)
+            {
+                return DoSampleInteger();
+            }
+
+            return int31;
+        }
+
         public void Dispose()
         {
 #if !NET35

src/Numerics/Random/Mcg31m1.cs

@@ -96,12 +96,9 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             double ret = _xn*Reciprocal;
             _xn = (_xn*Multiplier)%Modulus;

src/Numerics/Random/Mcg59.cs

@@ -97,12 +97,9 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             double ret = _xn*Reciprocal;
             _xn = (_xn*Multiplier)%Modulus;

src/Numerics/Random/MersenneTwister.cs

@@ -105,7 +105,7 @@ namespace MathNet.Numerics.Random
         /// <summary>
         /// Mersenne twister constant.
         /// </summary>
-        const double Reciprocal = 1.0/4294967296.0;
+        const double Reciprocal = 1.0/4294967296.0; // 1.0/(uint.MaxValue + 1.0)
 
         /// <summary>
         /// Mersenne twister constant.
@@ -320,16 +320,28 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             return genrand_int32()*Reciprocal;
         }
 
+        /// <summary>
+        /// Returns a random 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>
+        /// </summary>
+        protected sealed override int DoSampleInteger()
+        {
+            uint uint32 = genrand_int32();
+            int int31 = (int)(uint32 >> 1);
+            if (int31 == int.MaxValue)
+            {
+                return DoSampleInteger();
+            }
+
+            return int31;
+        }
+
         /* /// <summary>
         /// Generates a random number on [0,1) with 53-bit resolution.
         /// </summary>

src/Numerics/Random/Mrg32k3a.cs

@@ -110,12 +110,9 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             double xn = A12*_xn2 - A13*_xn3;
             double k = (long)(xn/Modulus1);

src/Numerics/Random/Palf.cs

@@ -62,7 +62,7 @@ namespace MathNet.Numerics.Random
         /// <summary>
         /// The multiplier to compute a double-precision floating point number [0, 1)
         /// </summary>
-        const double IntToDoubleMultiplier = 1.0/(int.MaxValue + 1.0);
+        const double Reciprocal = 1.0/4294967296.0; // 1.0/(uint.MaxValue + 1.0)
 
         /// <summary>
         /// Initializes a new instance of the <see cref="Palf"/> class using
@@ -211,20 +211,36 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             if (_k >= LongLag)
             {
                 Fill();
             }
 
-            var x = _x[_k++];
+            return _x[_k++] * Reciprocal;
-            return (int)(x >> 1)*IntToDoubleMultiplier;
+        }
+
+        /// <summary>
+        /// Returns a random 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>
+        /// </summary>
+        protected override int DoSampleInteger()
+        {
+            if (_k >= LongLag)
+            {
+                Fill();
+            }
+
+            uint uint32 = _x[_k++];
+            int int31 = (int)(uint32 >> 1);
+            if (int31 == int.MaxValue)
+            {
+                return DoSampleInteger();
+            }
+
+            return int31;
         }
 
         /// <summary>
@@ -272,7 +288,7 @@ namespace MathNet.Numerics.Random
                     k = 0;
                 }
 
-                values[i] = (int)(x[k++] >> 1)*IntToDoubleMultiplier;
+                values[i] = x[k++]*Reciprocal;
             }
         }
 
@@ -333,7 +349,7 @@ namespace MathNet.Numerics.Random
                     k = 0;
                 }
 
-                yield return (int)(x[k++] >> 1)*IntToDoubleMultiplier;
+                yield return x[k++]*Reciprocal;
             }
         }
     }

src/Numerics/Random/RandomExtensions.cs

@@ -188,7 +188,7 @@ namespace MathNet.Numerics.Random
         /// </remarks>
         public static long NextInt64(this System.Random rnd)
         {
-            var buffer = new byte[sizeof (long)];
+            var buffer = new byte[8];
 
             rnd.NextBytes(buffer);
             var candidate = BitConverter.ToInt64(buffer, 0);
@@ -216,7 +216,7 @@ namespace MathNet.Numerics.Random
         /// </remarks>
         public static int NextFullRangeInt32(this System.Random rnd)
         {
-            var buffer = new byte[sizeof (int)];
+            var buffer = new byte[4];
             rnd.NextBytes(buffer);
             return BitConverter.ToInt32(buffer, 0);
         }
@@ -236,7 +236,7 @@ namespace MathNet.Numerics.Random
         /// </remarks>
         public static long NextFullRangeInt64(this System.Random rnd)
         {
-            var buffer = new byte[sizeof (long)];
+            var buffer = new byte[8];
             rnd.NextBytes(buffer);
             return BitConverter.ToInt64(buffer, 0);
         }

src/Numerics/Random/RandomSource.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2015 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -117,17 +117,17 @@ namespace MathNet.Numerics.Random
         /// <returns>
         /// A 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>.
         /// </returns>
-        public override sealed int Next()
+        public sealed override int Next()
         {
             if (_threadSafe)
             {
                 lock (_lock)
                 {
-                    return (int)(DoSample()*int.MaxValue);
+                    return DoSampleInteger();
                 }
             }
 
-            return (int)(DoSample()*int.MaxValue);
+            return DoSampleInteger();
         }
 
         /// <summary>
@@ -136,36 +136,41 @@ namespace MathNet.Numerics.Random
         /// <param name="maxValue">The exclusive upper bound of the random number returned.</param>
         /// <returns>A 32-bit signed integer less than <paramref name="maxValue"/>.</returns>
         /// <exception cref="T:System.ArgumentOutOfRangeException"><paramref name="maxValue"/> is negative. </exception>
-        public override sealed int Next(int maxValue)
+        public sealed override int Next(int maxValue)
         {
             if (maxValue <= 0)
             {
                 throw new ArgumentException(Resources.ArgumentMustBePositive);
             }
 
+            if (maxValue == int.MaxValue)
+            {
+                return Next();
+            }
+
             if (_threadSafe)
             {
                 lock (_lock)
                 {
-                    return (int)(DoSample()*maxValue);
+                    return DoSampleInteger(0, maxValue);
                 }
             }
 
-            return (int)(DoSample()*maxValue);
+            return DoSampleInteger(0, maxValue);
         }
 
         /// <summary>
         /// Returns a random number within a specified range.
         /// </summary>
-        /// <param name="minValue">The inclusive lower bound of the random number returned.</param>
+        /// <param name="minInclusive">The inclusive lower bound of the random number returned.</param>
-        /// <param name="maxValue">The exclusive upper bound of the random number returned. <paramref name="maxValue"/> must be greater than or equal to <paramref name="minValue"/>.</param>
+        /// <param name="maxExclusive">The exclusive upper bound of the random number returned. <paramref name="maxExclusive"/> must be greater than or equal to <paramref name="minInclusive"/>.</param>
         /// <returns>
-        /// A 32-bit signed integer greater than or equal to <paramref name="minValue"/> and less than <paramref name="maxValue"/>; that is, the range of return values includes <paramref name="minValue"/> but not <paramref name="maxValue"/>. If <paramref name="minValue"/> equals <paramref name="maxValue"/>, <paramref name="minValue"/> is returned.
+        /// A 32-bit signed integer greater than or equal to <paramref name="minInclusive"/> and less than <paramref name="maxExclusive"/>; that is, the range of return values includes <paramref name="minInclusive"/> but not <paramref name="maxExclusive"/>. If <paramref name="minInclusive"/> equals <paramref name="maxExclusive"/>, <paramref name="minInclusive"/> is returned.
         /// </returns>
-        /// <exception cref="T:System.ArgumentOutOfRangeException"><paramref name="minValue"/> is greater than <paramref name="maxValue"/>. </exception>
+        /// <exception cref="T:System.ArgumentOutOfRangeException"><paramref name="minInclusive"/> is greater than <paramref name="maxExclusive"/>. </exception>
-        public override sealed int Next(int minValue, int maxValue)
+        public sealed override int Next(int minInclusive, int maxExclusive)
         {
-            if (minValue > maxValue)
+            if (minInclusive > maxExclusive)
             {
                 throw new ArgumentException(Resources.ArgumentMinValueGreaterThanMaxValue);
             }
@@ -174,20 +179,20 @@ namespace MathNet.Numerics.Random
             {
                 lock (_lock)
                 {
-                    return (int)(DoSample()*(maxValue - minValue)) + minValue;
+                    return DoSampleInteger(minInclusive, maxExclusive);
                 }
             }
 
-            return (int)(DoSample()*(maxValue - minValue)) + minValue;
+            return DoSampleInteger(minInclusive, maxExclusive);
         }
 
         /// <summary>
         /// Fills an array with random numbers within a specified range.
         /// </summary>
         /// <param name="values">The array to fill with random values.</param>
-        /// <param name="minValue">The inclusive lower bound of the random number returned.</param>
+        /// <param name="minInclusive">The inclusive lower bound of the random number returned.</param>
-        /// <param name="maxValue">The exclusive upper bound of the random number returned. <paramref name="maxValue"/> must be greater than or equal to <paramref name="minValue"/>.</param>
+        /// <param name="maxExclusive">The exclusive upper bound of the random number returned. <paramref name="maxExclusive"/> must be greater than or equal to <paramref name="minInclusive"/>.</param>
-        public void NextInt32s(int[] values, int minValue, int maxValue)
+        public void NextInt32s(int[] values, int minInclusive, int maxExclusive)
         {
             if (_threadSafe)
             {
@@ -195,7 +200,7 @@ namespace MathNet.Numerics.Random
                 {
                     for (var i = 0; i < values.Length; i++)
                     {
-                        values[i] = (int)(DoSample()*(maxValue - minValue)) + minValue;
+                        values[i] = DoSampleInteger(minInclusive, maxExclusive);
                     }
                 }
             }
@@ -203,7 +208,7 @@ namespace MathNet.Numerics.Random
             {
                 for (var i = 0; i < values.Length; i++)
                 {
-                    values[i] = (int)(DoSample()*(maxValue - minValue)) + minValue;
+                    values[i] = DoSampleInteger(minInclusive, maxExclusive);
                 }
             }
         }
@@ -211,19 +216,19 @@ namespace MathNet.Numerics.Random
         /// <summary>
         /// Returns an infinite sequence of random numbers within a specified range.
         /// </summary>
-        /// <param name="minValue">The inclusive lower bound of the random number returned.</param>
+        /// <param name="minInclusive">The inclusive lower bound of the random number returned.</param>
-        /// <param name="maxValue">The exclusive upper bound of the random number returned. <paramref name="maxValue"/> must be greater than or equal to <paramref name="minValue"/>.</param>
+        /// <param name="maxExclusive">The exclusive upper bound of the random number returned. <paramref name="maxExclusive"/> must be greater than or equal to <paramref name="minInclusive"/>.</param>
-        public IEnumerable<int> NextInt32Sequence(int minValue, int maxValue)
+        public IEnumerable<int> NextInt32Sequence(int minInclusive, int maxExclusive)
         {
             for (int i = 0; i < 64; i++)
             {
-                yield return Next(minValue, maxValue);
+                yield return Next(minInclusive, maxExclusive);
             }
 
             var buffer = new int[64];
             while (true)
             {
-                NextInt32s(buffer, minValue, maxValue);
+                NextInt32s(buffer, minInclusive, maxExclusive);
                 for (int i = 0; i < buffer.Length; i++)
                 {
                     yield return buffer[i];
@@ -249,7 +254,7 @@ namespace MathNet.Numerics.Random
                 {
                     for (var i = 0; i < buffer.Length; i++)
                     {
-                        buffer[i] = (byte)(((int)(DoSample()*int.MaxValue))%256);
+                        buffer[i] = (byte)(DoSampleInteger()%256);
                     }
                 }
 
@@ -258,7 +263,7 @@ namespace MathNet.Numerics.Random
 
             for (var i = 0; i < buffer.Length; i++)
             {
-                buffer[i] = (byte)(((int)(DoSample()*int.MaxValue))%256);
+                buffer[i] = (byte)(DoSampleInteger()%256);
             }
         }
 
@@ -266,7 +271,7 @@ namespace MathNet.Numerics.Random
         /// Returns a random number between 0.0 and 1.0.
         /// </summary>
         /// <returns>A double-precision floating point number greater than or equal to 0.0, and less than 1.0.</returns>
-        protected override sealed double Sample()
+        protected sealed override double Sample()
         {
             if (_threadSafe)
             {
@@ -280,11 +285,26 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
         protected abstract double DoSample();
+
+        /// <summary>
+        /// Returns a random 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>
+        /// </summary>
+        protected virtual int DoSampleInteger()
+        {
+            return (int)(DoSample() * int.MaxValue);
+        }
+
+        /// <summary>
+        /// Returns a random 32-bit signed integer within the specified range.
+        /// </summary>
+        /// <param name="minInclusive">The inclusive lower bound of the random number returned.</param>
+        /// <param name="maxExclusive">The exclusive upper bound of the random number returned. <paramref name="maxExclusive"/> must be greater than or equal to <paramref name="minInclusive"/>.</param>
+        protected virtual int DoSampleInteger(int minInclusive, int maxExclusive)
+        {
+            return (int)(DoSample()*(maxExclusive - minInclusive)) + minInclusive;
+        }
     }
 }

src/Numerics/Random/SystemRandomSource.cs

@@ -112,16 +112,31 @@ namespace MathNet.Numerics.Random
 #endif
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             return _random.NextDouble();
         }
 
+        /// <summary>
+        /// Returns a random 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>
+        /// </summary>
+        protected override int DoSampleInteger()
+        {
+            return _random.Next();
+        }
+
+        /// <summary>
+        /// Returns a random 32-bit signed integer within the specified range.
+        /// </summary>
+        /// <param name="minInclusive">The inclusive lower bound of the random number returned.</param>
+        /// <param name="maxExclusive">The exclusive upper bound of the random number returned. <paramref name="maxExclusive"/> must be greater than or equal to <paramref name="minInclusive"/>.</param>
+        protected override int DoSampleInteger(int minInclusive, int maxExclusive)
+        {
+            return _random.Next(minInclusive, maxExclusive);
+        }
+
         /// <summary>
         /// Fill an array with uniform random numbers greater than or equal to 0.0 and less than 1.0.
         /// WARNING: potentially very short random sequence length, can generate repeated partial sequences.

src/Numerics/Random/WH1982.cs

@@ -106,12 +106,9 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             _xn = (171*_xn)%Modx;
             _yn = (172*_yn)%Mody;

src/Numerics/Random/WH2006.cs

@@ -108,12 +108,9 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             _xn = 11600*_xn%Modx;
             _yn = 47003*_yn%Mody;

src/Numerics/Random/Xorshift.cs

@@ -68,7 +68,7 @@ namespace MathNet.Numerics.Random
         /// <summary>
         /// The multiplier to compute a double-precision floating point number [0, 1)
         /// </summary>
-        const double UlongToDoubleMultiplier = 1.0/(uint.MaxValue + 1.0);
+        const double UlongToDoubleMultiplier = 1.0/4294967296.0; // 1.0/(uint.MaxValue + 1.0)
 
         /// <summary>
         /// Seed or last but three unsigned random number.
@@ -252,12 +252,9 @@ namespace MathNet.Numerics.Random
         }
 
         /// <summary>
-        /// Returns a random number between 0.0 and 1.0.
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
         /// </summary>
-        /// <returns>
+        protected sealed override double DoSample()
-        /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0.
-        /// </returns>
-        protected override sealed double DoSample()
         {
             var t = (_a*_x) + _c;
             _x = _y;
@@ -267,6 +264,26 @@ namespace MathNet.Numerics.Random
             return _z*UlongToDoubleMultiplier;
         }
 
+        /// <summary>
+        /// Returns a random 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>
+        /// </summary>
+        protected override int DoSampleInteger()
+        {
+            var t = (_a * _x) + _c;
+            _x = _y;
+            _y = _z;
+            _c = t >> 32;
+            _z = t & 0xffffffff;
+            uint uint32 = (uint)_z;
+            int int31 = (int)(uint32 >> 1);
+            if (int31 == int.MaxValue)
+            {
+                return DoSampleInteger();
+            }
+
+            return int31;
+        }
+
         /// <summary>
         /// Fills an array with random numbers greater than or equal to 0.0 and less than 1.0.
         /// </summary>

src/UnitTests/Random/MersenneTwisterTests.cs

@@ -50,7 +50,17 @@ namespace MathNet.Numerics.UnitTests.Random
         /// Sample known values.
         /// </summary>
         [Test]
-        public void SampleKnownValues()
+        public void SampleKnownIntegerValuesSeed42()
+        {
+            var mt = new MersenneTwister(42);
+            Assert.AreEqual(mt.Next(), 804318771);
+        }
+
+        /// <summary>
+        /// Sample known values.
+        /// </summary>
+        [Test]
+        public void SampleKnownFloatingPointValues()
         {
             var mt = new MersenneTwister(0);
             Assert.AreEqual(mt.NextDouble(), 0.5488135023042560);