From b209ed9d58f3fbf816ff95aa035445e0c2ef3a86 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Wed, 30 Dec 2015 16:29:55 +0100 Subject: [PATCH] Random: RandomSource.DoSampleInteger for direct integer sampling without floating points #359 MersenneTwister integer sampling now matches standard implementation. --- src/Numerics/Random/CryptoRandomSource.cs | 24 ++++-- src/Numerics/Random/Mcg31m1.cs | 7 +- src/Numerics/Random/Mcg59.cs | 7 +- src/Numerics/Random/MersenneTwister.cs | 44 ++++++---- src/Numerics/Random/Mrg32k3a.cs | 7 +- src/Numerics/Random/Palf.cs | 36 ++++++--- src/Numerics/Random/RandomExtensions.cs | 6 +- src/Numerics/Random/RandomSource.cs | 84 ++++++++++++-------- src/Numerics/Random/SystemRandomSource.cs | 25 ++++-- src/Numerics/Random/WH1982.cs | 9 +-- src/Numerics/Random/WH2006.cs | 7 +- src/Numerics/Random/Xorshift.cs | 29 +++++-- src/UnitTests/Random/MersenneTwisterTests.cs | 12 ++- 13 files changed, 193 insertions(+), 104 deletions(-) diff --git a/src/Numerics/Random/CryptoRandomSource.cs b/src/Numerics/Random/CryptoRandomSource.cs index 700f5645..f17d71b2 100644 --- a/src/Numerics/Random/CryptoRandomSource.cs +++ b/src/Numerics/Random/CryptoRandomSource.cs @@ -86,18 +86,32 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { var bytes = new byte[4]; _crypto.GetBytes(bytes); return BitConverter.ToUInt32(bytes, 0)*Reciprocal; } + /// + /// Returns a random 32-bit signed integer greater than or equal to zero and less than + /// + 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 diff --git a/src/Numerics/Random/Mcg31m1.cs b/src/Numerics/Random/Mcg31m1.cs index 44ac6b0c..37d66564 100644 --- a/src/Numerics/Random/Mcg31m1.cs +++ b/src/Numerics/Random/Mcg31m1.cs @@ -96,12 +96,9 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { double ret = _xn*Reciprocal; _xn = (_xn*Multiplier)%Modulus; diff --git a/src/Numerics/Random/Mcg59.cs b/src/Numerics/Random/Mcg59.cs index 87d7cc61..3b0c5791 100644 --- a/src/Numerics/Random/Mcg59.cs +++ b/src/Numerics/Random/Mcg59.cs @@ -97,12 +97,9 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { double ret = _xn*Reciprocal; _xn = (_xn*Multiplier)%Modulus; diff --git a/src/Numerics/Random/MersenneTwister.cs b/src/Numerics/Random/MersenneTwister.cs index 50ffc5ac..8d38d4f0 100644 --- a/src/Numerics/Random/MersenneTwister.cs +++ b/src/Numerics/Random/MersenneTwister.cs @@ -28,10 +28,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // -/* +/* Original code's copyright and license: Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, - All rights reserved. + All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions @@ -44,8 +44,8 @@ notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. - 3. The names of its contributors may not be used to endorse or promote - products derived from this software without specific prior written + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS @@ -105,7 +105,7 @@ namespace MathNet.Numerics.Random /// /// Mersenne twister constant. /// - const double Reciprocal = 1.0/4294967296.0; + const double Reciprocal = 1.0/4294967296.0; // 1.0/(uint.MaxValue + 1.0) /// /// Mersenne twister constant. @@ -152,7 +152,7 @@ namespace MathNet.Numerics.Random /// /// The seed value. /// Uses the value of to - /// set whether the instance is thread safe. + /// set whether the instance is thread safe. public MersenneTwister(int seed) { init_genrand((uint)seed); @@ -246,8 +246,8 @@ namespace MathNet.Numerics.Random uint k = (_n > key_length ? _n : key_length); for (; k > 0; k--) { - _mt[i] = (_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30))*1664525)) + init_key[j] + j; //non linear - _mt[i] &= 0xffffffff; // for WORDSIZE > 32 machines + _mt[i] = (_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30))*1664525)) + init_key[j] + j; //non linear + _mt[i] &= 0xffffffff; // for WORDSIZE > 32 machines i++; j++; if (i >= _n) @@ -259,8 +259,8 @@ namespace MathNet.Numerics.Random } for (k = _n - 1; k > 0; k--) { - _mt[i] = (_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30))*1566083941)) - i; // non linear - _mt[i] &= 0xffffffff; // for WORDSIZE > 32 machines + _mt[i] = (_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30))*1566083941)) - i; // non linear + _mt[i] &= 0xffffffff; // for WORDSIZE > 32 machines i++; if (i >= _n) { @@ -269,7 +269,7 @@ namespace MathNet.Numerics.Random } } - _mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array + _mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array }*/ /* generates a random number on [0,0xffffffff]-interval */ @@ -320,16 +320,28 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { return genrand_int32()*Reciprocal; } + /// + /// Returns a random 32-bit signed integer greater than or equal to zero and less than + /// + protected sealed override int DoSampleInteger() + { + uint uint32 = genrand_int32(); + int int31 = (int)(uint32 >> 1); + if (int31 == int.MaxValue) + { + return DoSampleInteger(); + } + + return int31; + } + /* /// /// Generates a random number on [0,1) with 53-bit resolution. /// diff --git a/src/Numerics/Random/Mrg32k3a.cs b/src/Numerics/Random/Mrg32k3a.cs index 51c58e6b..424cbac1 100644 --- a/src/Numerics/Random/Mrg32k3a.cs +++ b/src/Numerics/Random/Mrg32k3a.cs @@ -110,12 +110,9 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { double xn = A12*_xn2 - A13*_xn3; double k = (long)(xn/Modulus1); diff --git a/src/Numerics/Random/Palf.cs b/src/Numerics/Random/Palf.cs index a1073fbf..99c34b40 100644 --- a/src/Numerics/Random/Palf.cs +++ b/src/Numerics/Random/Palf.cs @@ -62,7 +62,7 @@ namespace MathNet.Numerics.Random /// /// The multiplier to compute a double-precision floating point number [0, 1) /// - const double IntToDoubleMultiplier = 1.0/(int.MaxValue + 1.0); + const double Reciprocal = 1.0/4294967296.0; // 1.0/(uint.MaxValue + 1.0) /// /// Initializes a new instance of the class using @@ -211,20 +211,36 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { if (_k >= LongLag) { Fill(); } - var x = _x[_k++]; - return (int)(x >> 1)*IntToDoubleMultiplier; + return _x[_k++] * Reciprocal; + } + + /// + /// Returns a random 32-bit signed integer greater than or equal to zero and less than + /// + 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; } /// @@ -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; } } } diff --git a/src/Numerics/Random/RandomExtensions.cs b/src/Numerics/Random/RandomExtensions.cs index 0a165524..57c5dea9 100644 --- a/src/Numerics/Random/RandomExtensions.cs +++ b/src/Numerics/Random/RandomExtensions.cs @@ -188,7 +188,7 @@ namespace MathNet.Numerics.Random /// 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 /// 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 /// 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); } diff --git a/src/Numerics/Random/RandomSource.cs b/src/Numerics/Random/RandomSource.cs index f40ab8ee..bdb87086 100644 --- a/src/Numerics/Random/RandomSource.cs +++ b/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 /// /// A 32-bit signed integer greater than or equal to zero and less than . /// - 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(); } /// @@ -136,36 +136,41 @@ namespace MathNet.Numerics.Random /// The exclusive upper bound of the random number returned. /// A 32-bit signed integer less than . /// is negative. - 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); } /// /// Returns a random number within a specified range. /// - /// The inclusive lower bound of the random number returned. - /// The exclusive upper bound of the random number returned. must be greater than or equal to . + /// The inclusive lower bound of the random number returned. + /// The exclusive upper bound of the random number returned. must be greater than or equal to . /// - /// A 32-bit signed integer greater than or equal to and less than ; that is, the range of return values includes but not . If equals , is returned. + /// A 32-bit signed integer greater than or equal to and less than ; that is, the range of return values includes but not . If equals , is returned. /// - /// is greater than . - public override sealed int Next(int minValue, int maxValue) + /// is greater than . + 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); } /// /// Fills an array with random numbers within a specified range. /// /// The array to fill with random values. - /// The inclusive lower bound of the random number returned. - /// The exclusive upper bound of the random number returned. must be greater than or equal to . - public void NextInt32s(int[] values, int minValue, int maxValue) + /// The inclusive lower bound of the random number returned. + /// The exclusive upper bound of the random number returned. must be greater than or equal to . + 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 /// /// Returns an infinite sequence of random numbers within a specified range. /// - /// The inclusive lower bound of the random number returned. - /// The exclusive upper bound of the random number returned. must be greater than or equal to . - public IEnumerable NextInt32Sequence(int minValue, int maxValue) + /// The inclusive lower bound of the random number returned. + /// The exclusive upper bound of the random number returned. must be greater than or equal to . + public IEnumerable 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. /// /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - protected override sealed double Sample() + protected sealed override double Sample() { if (_threadSafe) { @@ -280,11 +285,26 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// protected abstract double DoSample(); + + /// + /// Returns a random 32-bit signed integer greater than or equal to zero and less than + /// + protected virtual int DoSampleInteger() + { + return (int)(DoSample() * int.MaxValue); + } + + /// + /// Returns a random 32-bit signed integer within the specified range. + /// + /// The inclusive lower bound of the random number returned. + /// The exclusive upper bound of the random number returned. must be greater than or equal to . + protected virtual int DoSampleInteger(int minInclusive, int maxExclusive) + { + return (int)(DoSample()*(maxExclusive - minInclusive)) + minInclusive; + } } } diff --git a/src/Numerics/Random/SystemRandomSource.cs b/src/Numerics/Random/SystemRandomSource.cs index b4967f70..ffa42add 100644 --- a/src/Numerics/Random/SystemRandomSource.cs +++ b/src/Numerics/Random/SystemRandomSource.cs @@ -112,16 +112,31 @@ namespace MathNet.Numerics.Random #endif /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { return _random.NextDouble(); } + /// + /// Returns a random 32-bit signed integer greater than or equal to zero and less than + /// + protected override int DoSampleInteger() + { + return _random.Next(); + } + + /// + /// Returns a random 32-bit signed integer within the specified range. + /// + /// The inclusive lower bound of the random number returned. + /// The exclusive upper bound of the random number returned. must be greater than or equal to . + protected override int DoSampleInteger(int minInclusive, int maxExclusive) + { + return _random.Next(minInclusive, maxExclusive); + } + /// /// 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. diff --git a/src/Numerics/Random/WH1982.cs b/src/Numerics/Random/WH1982.cs index 26045bd2..cafe970b 100644 --- a/src/Numerics/Random/WH1982.cs +++ b/src/Numerics/Random/WH1982.cs @@ -37,7 +37,7 @@ using System.Runtime; namespace MathNet.Numerics.Random { /// - /// Wichmann-Hill’s 1982 combined multiplicative congruential generator. + /// Wichmann-Hill’s 1982 combined multiplicative congruential generator. /// /// See: Wichmann, B. A. & Hill, I. D. (1982), "Algorithm AS 183: /// An efficient and portable pseudo-random number generator". Applied Statistics 31 (1982) 188-190 @@ -106,12 +106,9 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { _xn = (171*_xn)%Modx; _yn = (172*_yn)%Mody; diff --git a/src/Numerics/Random/WH2006.cs b/src/Numerics/Random/WH2006.cs index e73f9f55..27cb3b5e 100644 --- a/src/Numerics/Random/WH2006.cs +++ b/src/Numerics/Random/WH2006.cs @@ -108,12 +108,9 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { _xn = 11600*_xn%Modx; _yn = 47003*_yn%Mody; diff --git a/src/Numerics/Random/Xorshift.cs b/src/Numerics/Random/Xorshift.cs index d5d1d7f9..18c4ca11 100644 --- a/src/Numerics/Random/Xorshift.cs +++ b/src/Numerics/Random/Xorshift.cs @@ -68,7 +68,7 @@ namespace MathNet.Numerics.Random /// /// The multiplier to compute a double-precision floating point number [0, 1) /// - const double UlongToDoubleMultiplier = 1.0/(uint.MaxValue + 1.0); + const double UlongToDoubleMultiplier = 1.0/4294967296.0; // 1.0/(uint.MaxValue + 1.0) /// /// Seed or last but three unsigned random number. @@ -252,12 +252,9 @@ namespace MathNet.Numerics.Random } /// - /// 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. /// - /// - /// A double-precision floating point number greater than or equal to 0.0, and less than 1.0. - /// - protected override sealed double DoSample() + protected sealed override double DoSample() { var t = (_a*_x) + _c; _x = _y; @@ -267,6 +264,26 @@ namespace MathNet.Numerics.Random return _z*UlongToDoubleMultiplier; } + /// + /// Returns a random 32-bit signed integer greater than or equal to zero and less than + /// + 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; + } + /// /// Fills an array with random numbers greater than or equal to 0.0 and less than 1.0. /// diff --git a/src/UnitTests/Random/MersenneTwisterTests.cs b/src/UnitTests/Random/MersenneTwisterTests.cs index 1a41bd07..d3c669a7 100644 --- a/src/UnitTests/Random/MersenneTwisterTests.cs +++ b/src/UnitTests/Random/MersenneTwisterTests.cs @@ -50,7 +50,17 @@ namespace MathNet.Numerics.UnitTests.Random /// Sample known values. /// [Test] - public void SampleKnownValues() + public void SampleKnownIntegerValuesSeed42() + { + var mt = new MersenneTwister(42); + Assert.AreEqual(mt.Next(), 804318771); + } + + /// + /// Sample known values. + /// + [Test] + public void SampleKnownFloatingPointValues() { var mt = new MersenneTwister(0); Assert.AreEqual(mt.NextDouble(), 0.5488135023042560);