Improved Random Number Generators

archived comments edit

Via Dare's blog, I found this interesting post on Random Number Generation on Michael Brundage's website. My undergrad thesis was on the topic of pseudorandom number generation so I thought I'd take the two classes he provided for a quick spin.

Unfortunately, the C# samples did not compile as is. In his post he discusses how the C++ samples are optimized. I figured I might be able to use them to guide changes to the C# port and could post the results here. Please note that I have not tested them yet and need to verify that my changes were correct. Enjoy and let me know if I got anything wrong.

[CLSCompliant(false)]

public class MersenneTwister

{

    private ulong _index;

    private ulong[] _buffer = new ulong[624];

 

    /// <summary>

    /// Creates a new <see cref="MersenneTwister"/> instance.

    /// </summary>

    public MersenneTwister()

    {

        Random r = new Random();

        for (int i = 0; i < 624; i++)

            _buffer[i] = (ulong)r.Next();

        _index = 0;

    }

 

    /// <summary>

    /// Returns a random long integer.

    /// </summary>

    /// <returns></returns>

    public ulong Random()

    {

        if (_index == 624)

        {

            _index = 0;

            long i = 0;

            ulong s;

            for (; i < 624 - 397; i++)

            {

                s = (_buffer[i] & 0x80000000) | (_buffer[i+1] & 0x7FFFFFFF);

                _buffer[i] = _buffer[i + 397] ^ (s >> 1) ^ ((s & 1) * 0x9908B0DF);

            }

            for (; i < 623; i++)

            {

                s = (_buffer[i] & 0x80000000) | (_buffer[i+1] & 0x7FFFFFFF);

                _buffer[i] = _buffer[i - (624 - 397)] ^ (s >> 1) ^ ((s & 1) * 0x9908B0DF);

            }

 

            s = (_buffer[623] & 0x80000000) | (_buffer[0] & 0x7FFFFFFF);

            _buffer[623] = _buffer[396] ^ (s >> 1) ^ ((s & 1) * 0x9908B0DF);

        }

        return _buffer[_index++];

    }

}

[CLSCompliant(false)]

public sealed class R250Combined521

{

    private ulong r250_index;

    private ulong r521_index;

    private ulong[] r250_buffer = new ulong[250];

    private ulong[] r521_buffer = new ulong[521];

 

    /// <summary>

    /// Creates a new <see cref="R250Combined521"/> instance.

    /// </summary>

    public R250Combined521()

    {

        Random r = new Random();

        ulong i = 521;

        ulong mask1 = 1;

        ulong mask2 = 0xFFFFFFFF;

 

        while (i-- > 250)

        {

            r521_buffer[i] = (ulong)r.Next();

        }

        while (i-- > 31)

        {

            r250_buffer[i] = (ulong)r.Next();

            r521_buffer[i] = (ulong)r.Next();

        }

 

        /*

        Establish linear independence of the bit columns

        by setting the diagonal bits and clearing all bits above

        */

        while (i-- > 0)

        {

            r250_buffer[i] = (((uint)r.Next()) | mask1) & mask2;

            r521_buffer[i] = (((uint)r.Next()) | mask1) & mask2;

            mask2 = mask2 ^ mask1;

            //mask2 ^= mask1;

            mask1 >>= 1;

        }

        r250_buffer[0] = mask1;

        r521_buffer[0] = mask2;

        r250_index = 0;

        r521_index = 0;

    }

 

    /// <summary>

    /// Returns a random long integer.

    /// </summary>

    /// <returns></returns>

    public ulong random()

    {

        ulong i1 = r250_index;

        ulong i2 = r521_index;

 

        ulong j1 = i1 - (250-103);

        if (j1 < 0)

            j1 = i1 + 103;

        ulong j2 = i2 - (521-168);

        if (j2 < 0)

            j2 = i2 + 168;

 

        ulong r = (r250_buffer[j1] ^ r250_buffer[i1]);

        r250_buffer[i1] = r;

        ulong s = (r521_buffer[j2] ^ r521_buffer[i2]);

        r521_buffer[i2] = s;

 

        i1 = (i1 != 249) ? (i1 + 1) : 0;

        r250_index = i1;

        i2 = (i2 != 521) ? (i2 + 1) : 0;

        r521_index = i2;

 

        return r ^ s;

    }

}

Comments