Ported dnA random number generators.

Signed-off-by: jvangael <jurgen.vangael@gmail.com> Signed-off-by: jvangael <jurgen.vangael@gmail.com>
17 years ago · be85f223b4
16 changed files with 976 additions and 7 deletions
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@ -103,8 +103,14 @@
      <DependentUpon>Resources.resx</DependentUpon>
    </Compile>
    <Compile Include="Random\AbstractRandomNumberGenerator.cs" />
+    <Compile Include="Random\Mcg31m1.cs" />
+    <Compile Include="Random\Mcg59.cs" />
    <Compile Include="Random\MersenneTwister.cs" />
+    <Compile Include="Random\Mrg32k3a.cs" />
+    <Compile Include="Random\SystemCrypto.cs" />
    <Compile Include="Random\SystemRandomExtensions.cs" />
+    <Compile Include="Random\WH1982.cs" />
+    <Compile Include="Random\WH2006.cs" />
    <Compile Include="Sampling\Sample.Random.cs" />
    <Compile Include="Sampling\Sample.Chebyshev.cs" />
    <Compile Include="Sampling\Sample.Equidistant.cs" />
--- a/src/Numerics/Random/Mcg31m1.cs
+++ b/src/Numerics/Random/Mcg31m1.cs
@ -0,0 +1,99 @@
+// <copyright file="Mcg31m1.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Random
+{
+    using System;
+
+    /// <summary>
+    /// Multiplicative congruential generator using a modulus of 2^31-1 and a multiplier of 1132489760.
+    /// </summary>
+    public partial class Mcg31m1 : AbstractRandomNumberGenerator
+    {
+        private const ulong _modulus = 2147483647;
+        private const ulong _multiplier = 1132489760;
+        private const double _reciprocal = 1.0 / _modulus;
+        private ulong _xn;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg31m1"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        public Mcg31m1() : this((int) DateTime.Now.Ticks)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg31m1"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public Mcg31m1(bool threadSafe) : this((int)DateTime.Now.Ticks, threadSafe)
+        {
+            
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg31m1"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>If the seed value is zero, it is set to one. Uses the
+        /// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
+        /// set whether the instance is thread safe.</remarks>
+        public Mcg31m1(int seed) : this(seed, Control.ThreadSafeRandomNumberGenerators)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg31m1"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <param name="threadSafe">if set to <c>true</c>, the class is thread safe.</param>
+        public Mcg31m1(int seed, bool threadSafe) : base(threadSafe)
+        {
+            if (seed == 0)
+            {
+                seed = 1;
+            }
+            _xn = (uint) seed%_modulus;
+        }
+
+        /// <summary>
+        /// 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 double DoSample()
+        {
+            double ret = _xn*_reciprocal;
+            _xn = (_xn*_multiplier)%_modulus;
+            return ret;
+        }
+    }
+}
--- a/src/Numerics/Random/Mcg59.cs
+++ b/src/Numerics/Random/Mcg59.cs
@ -0,0 +1,100 @@
+// <copyright file="Mcg59.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Random
+{
+    using System;
+
+    /// <summary>
+    /// Multiplicative congruential generator using a modulus of 2^59 and a multiplier of 13^13.
+    /// </summary>
+    public partial class Mcg59 : AbstractRandomNumberGenerator
+    {
+        private const double _reciprocal = 1.0 / _modulus;
+        private const ulong _modulus = 576460752303423488;
+        private const ulong _multiplier = 302875106592253;
+        private ulong _xn;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg59"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        public Mcg59() : this((int) DateTime.Now.Ticks)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg59"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public Mcg59(bool threadSafe) : this((int)DateTime.Now.Ticks, threadSafe)
+        {
+            
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg59"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>If the seed value is zero, it is set to one. Uses the
+        /// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
+        /// set whether the instance is thread safe.</remarks>
+        public Mcg59(int seed) : this(seed, Control.ThreadSafeRandomNumberGenerators)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg59"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>The seed is set to 1, if the zero is used as the seed.</remarks>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public Mcg59(int seed, bool threadSafe) : base(threadSafe)
+        {
+            if( seed == 0)
+            {
+                seed = 1;
+            }
+            _xn = (uint) seed % _modulus;
+        }
+
+        /// <summary>
+        /// 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 double DoSample()
+        {
+            double ret = _xn * _reciprocal;
+            _xn = (_xn * _multiplier) % _modulus;
+            return ret;
+        }
+    }
+}
--- a/src/Numerics/Random/Mrg32k3a.cs
+++ b/src/Numerics/Random/Mrg32k3a.cs
@ -0,0 +1,136 @@
+// <copyright file="Mrg32k3a.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Random
+{
+    using System;
+
+    /// <summary>
+    /// A 32-bit combined multiple recursive generator with 2 components of order 3.
+    /// </summary>
+    ///<remarks>Based off of P. L'Ecuyer, "Combined Multiple Recursive Random Number Generators," Operations Research, 44, 5 (1996), 816--822. </remarks>
+    public partial class Mrg32k3a : AbstractRandomNumberGenerator
+    {
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg31m1"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        /// <remarks>If the seed value is zero, it is set to one. Uses the
+        /// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
+        /// set whether the instance is thread safe.</remarks>
+        public Mrg32k3a() : this((int)DateTime.Now.Ticks)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mcg31m1"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public Mrg32k3a(bool threadSafe) : this((int)DateTime.Now.Ticks, threadSafe)
+        {
+        }
+        private const double _a12 = 1403580;
+        private const double _a13 = 810728;
+        private const double _a21 = 527612;
+        private const double _a23 = 1370589;
+        private const double _modulus1 = 4294967087;
+        private const double _modulus2 = 4294944443;
+
+        private const double _reciprocal = 1.0/_modulus1;
+        private double _xn1 = 1;
+        private double _xn2 = 1;
+        private double _xn3;
+        private double _yn1 = 1;
+        private double _yn2 = 1;
+        private double _yn3 = 1;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mrg32k3a"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>If the seed value is zero, it is set to one. Uses the
+        /// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
+        /// set whether the instance is thread safe.</remarks>
+        public Mrg32k3a(int seed) : this(seed, Control.ThreadSafeRandomNumberGenerators)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Mrg32k3a"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <param name="threadSafe">if set to <c>true</c>, the class is thread safe.</param>
+        public Mrg32k3a(int seed, bool threadSafe) : base(threadSafe)
+        {
+            if (seed == 0)
+            {
+                seed = 1;
+            }
+            _xn3 = (uint)seed;
+        }
+
+
+        /// <summary>
+        /// 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 double DoSample()
+        {
+            double xn = _a12*_xn2 - _a13*_xn3;
+            double k = (long) (xn/_modulus1);
+            xn -= k*_modulus1;
+            if (xn < 0)
+            {
+                xn += _modulus1;
+            }
+
+            double yn = _a21*_yn1 - _a23*_yn3;
+            k = (long) (yn/_modulus2);
+            yn -= k*_modulus2;
+            if (yn < 0)
+            {
+                yn += _modulus2;
+            }
+            _xn3 = _xn2;
+            _xn2 = _xn1;
+            _xn1 = xn;
+            _yn3 = _yn2;
+            _yn2 = _yn1;
+            _yn1 = yn;
+
+            if (xn <= yn)
+            {
+                return (xn - yn + _modulus1)*_reciprocal;
+            }
+            return (xn - yn)*_reciprocal;
+        }
+    }
+}
--- a/src/Numerics/Random/SystemCrypto.cs
+++ b/src/Numerics/Random/SystemCrypto.cs
@ -0,0 +1,98 @@
+// <copyright file="SystemCrypto.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Random
+{
+    using System;
+    using System.Security.Cryptography;
+
+    /// <summary>
+    /// A random number generator based on the <see cref="System.Security.Cryptography.RandomNumberGenerator"/> class in the .NET library.
+    /// </summary>
+    public class SystemCryptoRandomNumberGenerator : AbstractRandomNumberGenerator
+    {
+        private const double mReciprocal = 1.0 / uint.MaxValue;
+        private readonly RandomNumberGenerator mRandom;
+
+        /// <summary>
+        /// Construct a new random number generator with a random seed.
+        /// </summary>
+        /// <remarks>Uses <see cref="System.Security.Cryptography.RNGCryptoServiceProvider"/> and uses the value of 
+        /// <see cref="Control.ThreadSafeRandomNumberGenerators"/> to set whether the instance is thread safe.</remarks>
+        public SystemCryptoRandomNumberGenerator(): this(new RNGCryptoServiceProvider(), Control.ThreadSafeRandomNumberGenerators)
+        {
+            mRandom = new RNGCryptoServiceProvider();
+        }
+
+        /// <summary>
+        /// Construct a new random number generator with random seed.
+        /// </summary>
+        /// <param name="rng">The <see cref="RandomNumberGenerator"/> to use.</param>
+        /// <remarks>Uses the value of  <see cref="Control.ThreadSafeRandomNumberGenerators"/> to set whether the instance is thread safe.</remarks>
+        public SystemCryptoRandomNumberGenerator(RandomNumberGenerator rng) : this(rng, Control.ThreadSafeRandomNumberGenerators)
+        {
+        }
+
+        /// <summary>
+        /// Construct a new random number generator with random seed.
+        /// </summary>
+        /// <remarks>Uses <see cref="System.Security.Cryptography.RNGCryptoServiceProvider"/></remarks>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public SystemCryptoRandomNumberGenerator(bool threadSafe): this(new RNGCryptoServiceProvider(), threadSafe)
+        {
+        }
+
+        /// <summary>
+        /// Construct a new random number generator with random seed.
+        /// </summary>
+        /// <param name="rng">The <see cref="RandomNumberGenerator"/> to use.</param>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public SystemCryptoRandomNumberGenerator(RandomNumberGenerator rng, bool threadSafe) : base(threadSafe)
+        {
+            if (rng == null)
+            {
+                throw new ArgumentNullException("rng");
+            }
+            mRandom = rng;
+        }
+
+
+        /// <summary>
+        /// 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 double DoSample()
+        {
+            byte[] bytes = new byte[4];
+            mRandom.GetBytes(bytes);
+            return BitConverter.ToUInt32(bytes, 0) * mReciprocal;
+        }
+    }
+}
--- a/src/Numerics/Random/WH1982.cs
+++ b/src/Numerics/Random/WH1982.cs
@ -0,0 +1,113 @@
+// <copyright file="WH1982.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Random
+{
+    using System;
+
+    /// <summary>
+    /// Wichmann-Hill’s 1982 combined multiplicative congruential generator. 
+    /// </summary>
+    /// <remarks>See: Wichmann, B. A. &amp; Hill, I. D. (1982), "Algorithm AS 183:
+    /// An efficient and portable pseudo-random number generator". Applied Statistics 31 (1982) 188-190
+    ///</remarks>
+    public class WH1982 : AbstractRandomNumberGenerator
+    {
+        private const uint _modx = 30269;
+        private const double _modx_recip = 1.0/_modx;
+        private const uint _mody = 30307;
+        private const double _mody_recip = 1.0/_mody;
+        private const uint _modz = 30323;
+        private const double _modz_recip = 1.0/_modz;
+        private uint _xn;
+        private uint _yn = 1;
+        private uint _zn = 1;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="WH1982"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        public WH1982() : this((int) DateTime.Now.Ticks)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="WH1982"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public WH1982(bool threadSafe)
+            : this((int) DateTime.Now.Ticks, threadSafe)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="WH1982"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>If the seed value is zero, it is set to one. Uses the
+        /// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
+        /// set whether the instance is thread safe.</remarks>
+        public WH1982(int seed) : this(seed, Control.ThreadSafeRandomNumberGenerators)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="WH1982"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>The seed is set to 1, if the zero is used as the seed.</remarks>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public WH1982(int seed, bool threadSafe)
+            : base(threadSafe)
+        {
+            if (seed == 0)
+            {
+                seed = 1;
+            }
+            _xn = (uint) seed%_modx;
+        }
+
+        /// <summary>
+        /// 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 double DoSample()
+        {
+            _xn = (171*_xn)%_modx;
+            _yn = (172*_yn)%_mody;
+            _zn = (170*_zn)%_modz;
+
+            double w = _xn*_modx_recip + _yn*_mody_recip + _zn*_modz_recip;
+            w -= (int) w;
+            return w;
+        }
+    }
+}
--- a/src/Numerics/Random/WH2006.cs
+++ b/src/Numerics/Random/WH2006.cs
@ -0,0 +1,117 @@
+// <copyright file="WH2006.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.Random
+{
+    using System;
+
+    /// <summary>
+    /// 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
+    /// </remarks>
+    public class WH2006 : AbstractRandomNumberGenerator
+    {
+        private const uint _modw = 2147483123;
+        private const double _modw_recip = 1.0/_modw;
+        private const uint _modx = 2147483579;
+        private const double _modx_recip = 1.0/_modx;
+        private const uint _mody = 2147483543;
+        private const double _mody_recip = 1.0/_mody;
+        private const uint _modz = 2147483423;
+        private const double _modz_recip = 1.0/_modz;
+        private ulong _wn = 1;
+        private ulong _xn;
+        private ulong _yn = 1;
+        private ulong _zn = 1;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="WH2006"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        public WH2006() : this((int) DateTime.Now.Ticks)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="WH2006"/> class using
+        /// the current time as the seed.
+        /// </summary>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public WH2006(bool threadSafe)
+            : this((int) DateTime.Now.Ticks, threadSafe)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="WH2006"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>If the seed value is zero, it is set to one. Uses the
+        /// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
+        /// set whether the instance is thread safe.</remarks>
+        public WH2006(int seed) : this(seed, Control.ThreadSafeRandomNumberGenerators)
+        {
+        }
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="WH2006"/> class.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <remarks>The seed is set to 1, if the zero is used as the seed.</remarks>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public WH2006(int seed, bool threadSafe)
+            : base(threadSafe)
+        {
+            if (seed == 0)
+            {
+                seed = 1;
+            }
+            _xn = (uint) seed%_modx;
+        }
+
+        /// <summary>
+        /// 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 double DoSample()
+        {
+            _xn = 11600*_xn%_modx;
+            _yn = 47003*_yn%_mody;
+            _zn = 23000*_zn%_modz;
+            _wn = 33000*_wn%_modw;
+
+            double u = _xn*_modx_recip + _yn*_mody_recip + _zn*_modz_recip + _wn*_modw_recip;
+            u -= (int) u;
+            return u;
+        }
+    }
+}
--- a/src/Numerics/SpecialFunctions.cs
+++ b/src/Numerics/SpecialFunctions.cs
@ -29,6 +29,7 @@
 namespace MathNet.Numerics
 {
    using System;
+    using Properties;

    /// <summary>
    /// This class implements a collection of special function evaluations for double precision. This class 
@ -66,13 +67,57 @@ namespace MathNet.Numerics
            };

        /// <summary>
-        /// Static Initializer for all special function routines
+        /// Initializes static members of the SpecialFunctions class.
        /// </summary>
        static SpecialFunctions()
        {
            InitializeFactorial();
        }

+        /// <summary>
+        /// Computes the <paramref name="t"/>'th Harmonic number.
+        /// </summary>
+        /// <param name="t">The Harmonic number which needs to be computed.</param>
+        /// <returns>The t'th Harmonic number.</returns>
+        public static double Harmonic(int t)
+        {
+            return Constants.EulerMascheroni + DiGamma(t + 1.0);
+        }
+
+        /// <summary>
+        /// Computes the logarithm of the Euler Beta function.
+        /// </summary>
+        /// <param name="z">A positive real number.</param>
+        /// <param name="w">A positive real number.</param>
+        /// <returns>The logarithm of the Euler Beta function evaluated at z,w.</returns>
+        /// <exception cref="ArgumentException">If <paramref name="z"/> or <paramref name="w"/> are not positive.</exception>
+        public static double BetaLn(double z, double w)
+        {
+            if (z <= 0.0)
+            {
+                throw new ArgumentException(Resources.ArgumentMustBePositive, "z");
+            }
+
+            if (w <= 0.0)
+            {
+                throw new ArgumentException(Resources.ArgumentMustBePositive, "w");
+            }
+
+            return GammaLn(z) + GammaLn(w) - GammaLn(z + w);
+        }
+
+        /// <summary>
+        /// Computes the Euler Beta function.
+        /// </summary>
+        /// <param name="z">The first Beta parameter, a positive real number.</param>
+        /// <param name="w">The second Beta parameter, a positive real number.</param>
+        /// <returns>The Euler Beta function evaluated at z,w.</returns>
+        /// <exception cref="ArgumentException">If <paramref name="z"/> or <paramref name="w"/> are not positive.</exception>
+        public static double Beta(double z, double w)
+        {
+            return System.Math.Exp(BetaLn(z, w));
+        }
+
        /// <summary>
        /// Computes the logarithm of the Gamma function. 
        /// </summary>
@ -235,6 +280,16 @@ namespace MathNet.Numerics
                return Double.NaN;
            }

+            if (Double.IsNegativeInfinity(p))
+            {
+                return 0.0;
+            }
+
+            if (Double.IsPositiveInfinity(p))
+            {
+                return Double.PositiveInfinity;
+            }
+
            double x = Math.Exp(p);
            for (double d = 1.0; d > 1.0e-15; d /= 2.0)
            {
@ -249,11 +304,6 @@ namespace MathNet.Numerics
            throw new NotImplementedException();
        }

-        public static double BetaLn(double a, double b)
-        {
-            throw new NotImplementedException();
-        }
-
        public static double BetaRegularized(double a, double b, double x)
        {
            throw new NotImplementedException();
--- a/src/UnitTests/Random/Mcg31m1Tests.cs
+++ b/src/UnitTests/Random/Mcg31m1Tests.cs
@ -0,0 +1,39 @@
+// <copyright file="Mcg31m1Tests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.RandomTests
+{
+    using MbUnit.Framework;
+    using MathNet.Numerics.Random;
+
+    [TestFixture]
+    public class Mcg31m1Tests : RandomTests
+    {
+        public Mcg31m1Tests() : base(typeof(Mcg31m1)){}
+    }
+}
--- a/src/UnitTests/Random/Mcg59Tests.cs
+++ b/src/UnitTests/Random/Mcg59Tests.cs
@ -0,0 +1,39 @@
+// <copyright file="Mcg59Tests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.RandomTests
+{
+    using MbUnit.Framework;
+    using MathNet.Numerics.Random;
+
+    [TestFixture]
+    public class Mcg59Tests : RandomTests
+    {
+        public Mcg59Tests() : base(typeof(Mcg59)) { }
+    }
+}
--- a/src/UnitTests/Random/MersenneTwisterTests.cs
+++ b/src/UnitTests/Random/MersenneTwisterTests.cs
@ -1,4 +1,4 @@
-// <copyright file="SystemRandomExtensionTests.cs" company="Math.NET">
+// <copyright file="MersenneTwisterTests.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
 //
--- a/src/UnitTests/Random/Mrg32k3aTests.cs
+++ b/src/UnitTests/Random/Mrg32k3aTests.cs
@ -0,0 +1,41 @@
+// <copyright file="Mrg32k3aTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.RandomTests
+{
+    using MbUnit.Framework;
+    using MathNet.Numerics.Random;
+
+    [TestFixture]
+    public class Mrg32k3aTests : RandomTests
+    {
+        public Mrg32k3aTests() : base(typeof (Mrg32k3a))
+        {
+        }
+    }
+}
--- a/src/UnitTests/Random/SystemCryptoTests.cs
+++ b/src/UnitTests/Random/SystemCryptoTests.cs
@ -0,0 +1,41 @@
+// <copyright file="SystemCryptoTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.RandomTests
+{
+    using MbUnit.Framework;
+    using MathNet.Numerics.Random;
+
+    [TestFixture]
+    public class SystemCryptoRandomNumberGeneratorTests : RandomTests
+    {
+        public SystemCryptoRandomNumberGeneratorTests() : base(typeof (SystemCryptoRandomNumberGenerator))
+        {
+        }
+    }
+}
--- a/src/UnitTests/Random/WH1982Tests.cs
+++ b/src/UnitTests/Random/WH1982Tests.cs
@ -0,0 +1,42 @@
+// <copyright file="WH1982Tests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.RandomTests
+{
+    using MbUnit.Framework;
+    using MathNet.Numerics.Random;
+
+    [TestFixture]
+    public class WH1982Tests : RandomTests
+    {
+        public WH1982Tests()
+            : base(typeof (WH1982))
+        {
+        }
+    }
+}
--- a/src/UnitTests/Random/WH2006Tests.cs
+++ b/src/UnitTests/Random/WH2006Tests.cs
@ -0,0 +1,42 @@
+// <copyright file="WH2006Tests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.RandomTests
+{
+    using MbUnit.Framework;
+    using MathNet.Numerics.Random;
+
+    [TestFixture]
+    public class WH2006Tests : RandomTests
+    {
+        public WH2006Tests()
+            : base(typeof(WH2006))
+        {
+        }
+    }
+}
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@ -94,9 +94,15 @@
    <Compile Include="NumberTheoryTests\IntegerTheoryTest.cs" />
    <Compile Include="PrecisionTest.cs" />
    <Compile Include="Properties\AssemblyInfo.cs" />
+    <Compile Include="Random\Mcg31m1Tests.cs" />
+    <Compile Include="Random\Mcg59Tests.cs" />
    <Compile Include="Random\MersenneTwisterTests.cs" />
+    <Compile Include="Random\Mrg32k3aTests.cs" />
    <Compile Include="Random\RandomTests.cs" />
+    <Compile Include="Random\SystemCryptoTests.cs" />
    <Compile Include="Random\SystemRandomExtensionTests.cs" />
+    <Compile Include="Random\WH1982Tests.cs" />
+    <Compile Include="Random\WH2006Tests.cs" />
    <Compile Include="SortingTests.cs" />
    <Compile Include="SpecialFunctionsTests\ErfTests.cs" />
    <Compile Include="SpecialFunctionsTests\FactorialTest.cs" />