Random: static SampleSequence routines

src/Numerics/Random/CryptoRandomSource.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System.Collections.Generic;
+
 #if !PORTABLE
 
 using System;
@@ -103,6 +105,40 @@ namespace MathNet.Numerics.Random
             _crypto.Dispose();
 #endif
         }
+
+        /// <summary>
+        /// Returns an array of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        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;
+        }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> 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;
+                }
+            }
+        }
     }
 }
 

src/Numerics/Random/Mcg31m1.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System.Collections.Generic;
+
 namespace MathNet.Numerics.Random
 {
     /// <summary>
@@ -119,5 +121,23 @@ namespace MathNet.Numerics.Random
             }
             return data;
         }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> SampleSequence(int seed)
+        {
+            if (seed == 0)
+            {
+                seed = 1;
+            }
+            ulong xn = (uint)seed%Modulus;
+
+            while (true)
+            {
+                yield return xn*Reciprocal;
+                xn = (xn*Multiplier)%Modulus;
+            }
+        }
     }
 }

src/Numerics/Random/Mcg59.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System.Collections.Generic;
+
 namespace MathNet.Numerics.Random
 {
     /// <summary>
@@ -121,5 +123,23 @@ namespace MathNet.Numerics.Random
             }
             return data;
         }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> SampleSequence(int seed)
+        {
+            if (seed == 0)
+            {
+                seed = 1;
+            }
+            ulong xn = (uint)seed%Modulus;
+
+            while (true)
+            {
+                yield return xn*Reciprocal;
+                xn = (xn*Multiplier)%Modulus;
+            }
+        }
     }
 }

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;
         }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> 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;
+            }
+        }
     }
 }

src/Numerics/Random/Mrg32k3a.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System.Collections.Generic;
+
 namespace MathNet.Numerics.Random
 {
     /// <summary>
@@ -180,5 +182,45 @@ namespace MathNet.Numerics.Random
             }
             return data;
         }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> 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;
+            }
+        }
     }
 }

src/Numerics/Random/Palf.cs

@@ -29,6 +29,7 @@
 // </copyright>
 
 using System;
+using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 
 namespace MathNet.Numerics.Random
@@ -270,5 +271,52 @@ namespace MathNet.Numerics.Random
             }
             return data;
         }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> 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;
+            }
+        }
     }
 }

src/Numerics/Random/RandomSource.cs

@@ -29,6 +29,7 @@
 // </copyright>
 
 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();
         }
+
+        /// <summary>
+        /// Returns an array of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        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;
+        }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> SampleSequence(int systemRandomSeed)
+        {
+            var rnd = new System.Random(systemRandomSeed);
+
+            while (true)
+            {
+                yield return rnd.NextDouble();
+            }
+        }
     }
 }

src/Numerics/Random/WH1982.cs

@@ -28,6 +28,8 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System.Collections.Generic;
+
 namespace MathNet.Numerics.Random
 {
     /// <summary>
@@ -139,5 +141,29 @@ namespace MathNet.Numerics.Random
             }
             return data;
         }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> 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;
+            }
+        }
     }
 }

src/Numerics/Random/WH2006.cs

@@ -28,10 +28,12 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+using System.Collections.Generic;
+
 namespace MathNet.Numerics.Random
 {
     /// <summary>
-    /// Wichmann-Hill’s 2006 combined multiplicative congruential generator. 
+    /// Wichmann-Hill’s 2006 combined multiplicative congruential generator.
     /// </summary>
     /// <remarks>See: Wichmann, B. A. &amp; Hill, I. D. (2006), "Generating good pseudo-random numbers".
     /// Computational Statistics &amp; Data Analysis 51:3 (2006) 1614-1622
@@ -144,5 +146,31 @@ namespace MathNet.Numerics.Random
             }
             return data;
         }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> 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;
+            }
+        }
     }
 }

src/Numerics/Random/Xorshift.cs

@@ -29,6 +29,7 @@
 // </copyright>
 
 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);
 
         /// <summary>
-        /// Seed or last but three unsigned random number. 
+        /// Seed or last but three unsigned random number.
         /// </summary>
         ulong _x;
 
         /// <summary>
-        /// Last but two unsigned random number. 
+        /// Last but two unsigned random number.
         /// </summary>
         ulong _y;
 
         /// <summary>
-        /// Last but one unsigned random number. 
+        /// Last but one unsigned random number.
         /// </summary>
         ulong _z;
 
         /// <summary>
-        /// The value of the carry over. 
+        /// The value of the carry over.
         /// </summary>
         ulong _c;
 
@@ -290,5 +291,32 @@ namespace MathNet.Numerics.Random
             }
             return data;
         }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        public static IEnumerable<double> 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;
+            }
+        }
     }
 }

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();
             }
         }