diff --git a/src/Numerics/Random/CryptoRandomSource.cs b/src/Numerics/Random/CryptoRandomSource.cs index b1de66b9..8bc864b9 100644 --- a/src/Numerics/Random/CryptoRandomSource.cs +++ b/src/Numerics/Random/CryptoRandomSource.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System.Collections.Generic; + #if !PORTABLE using System; @@ -103,6 +105,40 @@ namespace MathNet.Numerics.Random _crypto.Dispose(); #endif } + + /// + /// Returns an array of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static double[] Samples(int length) + { + var rnd = new RNGCryptoServiceProvider(); + var bytes = new byte[length*4]; + rnd.GetBytes(bytes); + var data = new double[length]; + for (int i = 0; i < data.Length; i++) + { + data[i] = BitConverter.ToUInt32(bytes, i*4)*Reciprocal; + } + return data; + } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence() + { + var rnd = new RNGCryptoServiceProvider(); + var buffer = new byte[1024*4]; + + while (true) + { + rnd.GetBytes(buffer); + for (int i = 0; i < buffer.Length; i += 4) + { + yield return BitConverter.ToUInt32(buffer, i)*Reciprocal; + } + } + } } } diff --git a/src/Numerics/Random/Mcg31m1.cs b/src/Numerics/Random/Mcg31m1.cs index e9d9e402..6770f87c 100644 --- a/src/Numerics/Random/Mcg31m1.cs +++ b/src/Numerics/Random/Mcg31m1.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System.Collections.Generic; + namespace MathNet.Numerics.Random { /// @@ -119,5 +121,23 @@ namespace MathNet.Numerics.Random } return data; } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int seed) + { + if (seed == 0) + { + seed = 1; + } + ulong xn = (uint)seed%Modulus; + + while (true) + { + yield return xn*Reciprocal; + xn = (xn*Multiplier)%Modulus; + } + } } } diff --git a/src/Numerics/Random/Mcg59.cs b/src/Numerics/Random/Mcg59.cs index 0e5d2f01..9b20083c 100644 --- a/src/Numerics/Random/Mcg59.cs +++ b/src/Numerics/Random/Mcg59.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System.Collections.Generic; + namespace MathNet.Numerics.Random { /// @@ -121,5 +123,23 @@ namespace MathNet.Numerics.Random } return data; } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int seed) + { + if (seed == 0) + { + seed = 1; + } + ulong xn = (uint)seed%Modulus; + + while (true) + { + yield return xn*Reciprocal; + xn = (xn*Multiplier)%Modulus; + } + } } } diff --git a/src/Numerics/Random/MersenneTwister.cs b/src/Numerics/Random/MersenneTwister.cs index 2aeeab7f..1621163b 100644 --- a/src/Numerics/Random/MersenneTwister.cs +++ b/src/Numerics/Random/MersenneTwister.cs @@ -67,6 +67,7 @@ */ using System; +using System.Collections.Generic; using System.Threading; namespace MathNet.Numerics.Random @@ -383,5 +384,56 @@ namespace MathNet.Numerics.Random } return data; } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int seed) + { + uint[] t = new uint[624]; + int k; + uint s = (uint)seed; + + t[0] = s & 0xffffffff; + for (k = 1; k < N; k++) + { + t[k] = (1812433253*(t[k - 1] ^ (t[k - 1] >> 30)) + (uint)k); + t[k] &= 0xffffffff; + } + + while (true) + { + uint y; + + if (k >= N) + { + int kk; + for (kk = 0; kk < N - M; kk++) + { + y = (t[kk] & UpperMask) | (t[kk + 1] & LowerMask); + t[kk] = t[kk + M] ^ (y >> 1) ^ Mag01[y & 0x1]; + } + for (; kk < N - 1; kk++) + { + y = (t[kk] & UpperMask) | (t[kk + 1] & LowerMask); + t[kk] = t[kk + (M - N)] ^ (y >> 1) ^ Mag01[y & 0x1]; + } + y = (t[N - 1] & UpperMask) | (t[0] & LowerMask); + t[N - 1] = t[M - 1] ^ (y >> 1) ^ Mag01[y & 0x1]; + + k = 0; + } + + y = t[k++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680; + y ^= (y << 15) & 0xefc60000; + y ^= (y >> 18); + + yield return y*Reciprocal; + } + } } } diff --git a/src/Numerics/Random/Mrg32k3a.cs b/src/Numerics/Random/Mrg32k3a.cs index cba2ad6f..a52b08be 100644 --- a/src/Numerics/Random/Mrg32k3a.cs +++ b/src/Numerics/Random/Mrg32k3a.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System.Collections.Generic; + namespace MathNet.Numerics.Random { /// @@ -180,5 +182,45 @@ namespace MathNet.Numerics.Random } return data; } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int seed) + { + double x1 = 1; + double x2 = 1; + double x3 = (uint)seed; + double y1 = 1; + double y2 = 1; + double y3 = 1; + + while (true) + { + double xn = A12*x2 - A13*x3; + double k = (long)(xn/Modulus1); + xn -= k*Modulus1; + if (xn < 0) + { + xn += Modulus1; + } + + double yn = A21*y1 - A23*y3; + k = (long)(yn/Modulus2); + yn -= k*Modulus2; + if (yn < 0) + { + yn += Modulus2; + } + x3 = x2; + x2 = x1; + x1 = xn; + y3 = y2; + y2 = y1; + y1 = yn; + + yield return xn <= yn ? (xn - yn + Modulus1)*Reciprocal : (xn - yn)*Reciprocal; + } + } } } diff --git a/src/Numerics/Random/Palf.cs b/src/Numerics/Random/Palf.cs index 7a0e8921..40fa6f65 100644 --- a/src/Numerics/Random/Palf.cs +++ b/src/Numerics/Random/Palf.cs @@ -29,6 +29,7 @@ // using System; +using System.Collections.Generic; using MathNet.Numerics.Properties; namespace MathNet.Numerics.Random @@ -270,5 +271,52 @@ namespace MathNet.Numerics.Random } return data; } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int seed) + { + if (seed == 0) + { + seed = 1; + } + + int threads = Control.NumberOfParallelWorkerThreads; + const int shortLag = DefaultShortLag; + var longLag = DefaultLongLag; + + // Align LongLag to number of worker threads. + if (longLag%threads != 0) + { + longLag = ((longLag/threads) + 1)*threads; + } + + var x = Generate.Map(MersenneTwister.Samples(longLag, seed), uniform => (uint)(uniform*uint.MaxValue)); + var k = longLag; + + while (true) + { + if (k >= longLag) + { + for (int index = 0; index < threads; index++) + { + // Two loops to avoid costly modulo operations + for (var j = index; j < shortLag; j = j + threads) + { + x[j] += x[j + (longLag - shortLag)]; + } + + for (var j = shortLag + index; j < longLag; j = j + threads) + { + x[j] += x[j - shortLag - index]; + } + } + k = 0; + } + + yield return (int)(x[k++] >> 1)*IntToDoubleMultiplier; + } + } } } diff --git a/src/Numerics/Random/RandomSource.cs b/src/Numerics/Random/RandomSource.cs index 2534ada7..569fc310 100644 --- a/src/Numerics/Random/RandomSource.cs +++ b/src/Numerics/Random/RandomSource.cs @@ -29,6 +29,7 @@ // using System; +using System.Collections.Generic; using MathNet.Numerics.Properties; namespace MathNet.Numerics.Random @@ -262,5 +263,33 @@ namespace MathNet.Numerics.Random { return base.Sample(); } + + /// + /// Returns an array of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static double[] Samples(int length, int systemRandomSeed) + { + var rnd = new System.Random(systemRandomSeed); + + var data = new double[length]; + for (int i = 0; i < data.Length; i++) + { + data[i] = rnd.NextDouble(); + } + return data; + } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int systemRandomSeed) + { + var rnd = new System.Random(systemRandomSeed); + + while (true) + { + yield return rnd.NextDouble(); + } + } } } diff --git a/src/Numerics/Random/WH1982.cs b/src/Numerics/Random/WH1982.cs index 476f3827..15083501 100644 --- a/src/Numerics/Random/WH1982.cs +++ b/src/Numerics/Random/WH1982.cs @@ -28,6 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System.Collections.Generic; + namespace MathNet.Numerics.Random { /// @@ -139,5 +141,29 @@ namespace MathNet.Numerics.Random } return data; } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int seed) + { + if (seed == 0) + { + seed = 1; + } + uint xn = (uint)seed%Modx; + uint yn = 1; + uint zn = 1; + + while (true) + { + xn = (171*xn)%Modx; + yn = (172*yn)%Mody; + zn = (170*zn)%Modz; + + double w = xn*ModxRecip + yn*ModyRecip + zn*ModzRecip; + yield return w - (int)w; + } + } } } diff --git a/src/Numerics/Random/WH2006.cs b/src/Numerics/Random/WH2006.cs index e1bfff60..37afd8a5 100644 --- a/src/Numerics/Random/WH2006.cs +++ b/src/Numerics/Random/WH2006.cs @@ -28,10 +28,12 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System.Collections.Generic; + namespace MathNet.Numerics.Random { /// - /// Wichmann-Hill’s 2006 combined multiplicative congruential generator. + /// Wichmann-Hill’s 2006 combined multiplicative congruential generator. /// /// See: Wichmann, B. A. & Hill, I. D. (2006), "Generating good pseudo-random numbers". /// Computational Statistics & Data Analysis 51:3 (2006) 1614-1622 @@ -144,5 +146,31 @@ namespace MathNet.Numerics.Random } return data; } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int seed) + { + if (seed == 0) + { + seed = 1; + } + ulong wn = 1; + ulong xn = (uint)seed%Modx; + ulong yn = 1; + ulong zn = 1; + + while (true) + { + xn = 11600*xn%Modx; + yn = 47003*yn%Mody; + zn = 23000*zn%Modz; + wn = 33000*wn%Modw; + + double u = xn*ModxRecip + yn*ModyRecip + zn*ModzRecip + wn*ModwRecip; + yield return u - (int)u; + } + } } } diff --git a/src/Numerics/Random/Xorshift.cs b/src/Numerics/Random/Xorshift.cs index f6df39b4..c42c95cd 100644 --- a/src/Numerics/Random/Xorshift.cs +++ b/src/Numerics/Random/Xorshift.cs @@ -29,6 +29,7 @@ // using System; +using System.Collections.Generic; using MathNet.Numerics.Properties; namespace MathNet.Numerics.Random @@ -66,22 +67,22 @@ namespace MathNet.Numerics.Random const double UlongToDoubleMultiplier = 1.0/(uint.MaxValue + 1.0); /// - /// Seed or last but three unsigned random number. + /// Seed or last but three unsigned random number. /// ulong _x; /// - /// Last but two unsigned random number. + /// Last but two unsigned random number. /// ulong _y; /// - /// Last but one unsigned random number. + /// Last but one unsigned random number. /// ulong _z; /// - /// The value of the carry over. + /// The value of the carry over. /// ulong _c; @@ -290,5 +291,32 @@ namespace MathNet.Numerics.Random } return data; } + + /// + /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0. + /// + public static IEnumerable SampleSequence(int seed, ulong a = ASeed, ulong c = CSeed, ulong x1 = YSeed, ulong x2 = ZSeed) + { + if (a <= c) + { + throw new ArgumentException(string.Format(Resources.ArgumentOutOfRangeGreater, "a", "c"), "a"); + } + + if (seed == 0) + { + seed = 1; + } + ulong x = (uint)seed; + + while (true) + { + var t = (a*x) + c; + x = x1; + x1 = x2; + c = t >> 32; + x2 = t & 0xffffffff; + yield return x2*UlongToDoubleMultiplier; + } + } } } diff --git a/src/UnitTests/Random/RandomTests.cs b/src/UnitTests/Random/RandomTests.cs index 4eeb98b5..7030cd56 100644 --- a/src/UnitTests/Random/RandomTests.cs +++ b/src/UnitTests/Random/RandomTests.cs @@ -78,9 +78,11 @@ namespace MathNet.Numerics.UnitTests.Random // make sure are within 10% of the expected sum. Assert.IsTrue(sum >= (N/2.0) - (.05*N)); Assert.IsTrue(sum <= (N/2.0) + (.05*N)); - if (random is IDisposable) + + var disposable = random as IDisposable; + if (disposable != null) { - ((IDisposable)random).Dispose(); + disposable.Dispose(); } }