Pagina's

2013/12/09

WELL512 scaled random numbers part 3


/*
"WELL" stands for "Well Equidistributed Long-period Linear",
"512" stands for the period: 2^512
Chris Lomont's version is used, visit:

    http://www.lomont.org/Math/Papers/2008/Lomont_PRNG_2008.pdf
  
and read about the advantages of the "WELL512" PRNG.
A "state" array of 16 uints is used, it's seeded by another PRNG: "System.Random".
"System.Random" uses a seed, the absolute value of an int.
Sequences of random numbers generated by System.Random, using different seed values,
are correlated, the left side of the image above clearly shows patterns.
The right side looks pretty random, sequences appear to be less correlated.
In other words: It is possible to regenerate a high quality random sequence,
using the same seed value. Different sequences, using different seed value's, 
seem to be less correlated.


    Results:
    
    10^7 times rand.Next() takes: 107 ms (~11 ns per random number)
    10^7 times rand.Next(int someValue) takes: 236 ms (~24 ns ...)
    
    |--------|------------------------------------------|       rand.Next()    11 ns      
    |        |                 ms 10^7 *                |      rand.Next(i)    24 ns              
    |        |--------------|-------------|-------------|      well_uint(u)    44 ns, mean  ~32 ns
    |      u | well_uint(u) | gen_uint(u) | g_uint_c(u) |       gen_uint(u)    50 ns, mean  ~35 ns    
    |--------|-------|------|-------|-----|------|------|       g_uint_c(u)   280 ns, mean ~250 ns
    |      0 |    24 |      |    24 |     |   24 |      |         
    |      1 |       |   74 |       |  73 |      |  139 |          
    |      2 |   121 |      |   124 |     |  296 |      |            
    |      3 |       |   80 |       |  78 |      |  211 |       
    |      4 |   222 |      |   242 |     |  518 |      |                 
    |      7 |       |   82 |       |  83 |      |  282 |               
    |      8 |   272 |      |   290 |     |  680 |      |     
    |  2^7-1 |       |   98 |       | 102 |      |  562 |     
    |  2^7-0 |   344 |      |   380 |     | 1079 |      |           
    | 2^15-1 |       |  131 |       | 130 |      | 1123 |               
    | 2^15-0 |   379 |      |   414 |     | 1634 |      |               
    | 2^30-1 |       |  193 |       | 193 |      | 2164 |  
    | 2^30-0 |   431 |      |   494 |     | 2672 |      |  
    | 2^31-1 |       |  195 |       | 210 |      | 2251 |  
    | 2^31-0 |   427 |      |   461 |     | 2747 |      |  
    | 2^32-1 |       |  196 |       | 211 |      | 2307 |  
    |--------|-------|------|-------|-----|------|------|


http://www.cacert.at/cgi-bin/rngresults  find and hover mouse over "well512" for the code.

|---------|---------|--------|----------|----------|--------|---------|----------|----------|----------|----------|
|  Gene-  |  Speed  | Upload | Entropy  | Birthday | Matrix | 6*8 Ma- | Minimum  | Random   | Squeeze  | Overlap- |
|  rator  | (MiB/s) |  Size  | (->8)    | Spacing  | Ranks  | trix R  | Distance | Spheres  |          | ping Sums|
|-------- |---------|--------|----------|----------|--------|---------|----------|----------|----------|----------|
| WELL512 |   400   | 19 MiB | 7.999992 | 0.446363 |  0.636 |  0.879  | 0.967495 | 0.690789 | 0.942252 | 0.005987 |
|  Random |   370   | 17 MiB | 7.999989 | 0.000000 |  0.392 |  0.136  | 0.063920 | 0.905097 | 0.098440 | 0.005467 |
|     CSP |    34   | 15 MiB | 7.999987 | 0.032102 |  0.084 |  0.942  | 0.894277 | 0.036221 | 0.706538 | 0.023523 |
|-------- |---------|--------|----------|----------|--------|---------|----------|----------|----------|----------|  
  
*/

using System;
using System.Drawing; // add a ref

class scaled_well512
{
    private static uint[] state = new uint[16];
    private static int index = -1;

    private static void init_well(int seed)
    {
        seed = Math.Abs(seed);
        index = seed & 15;
        var rand = new Random(seed);
        uint u = (uint)rand.Next(1 << 30);
        for (int i = 0; i < 15; i++)
        {
            state[i] = (uint)rand.Next(1 << 30) << 2 | u & 3u;
            state[i] ^= (uint)rand.Next(1 << 30);
            state[i] ^= (uint)rand.Next(1 << 30);
            u >>= 2;
        }
        state[15] = (uint)rand.Next(1 << 30) << 2 | (uint)rand.Next(4);
        state[15] ^= (uint)rand.Next(1 << 30);
        state[15] ^= (uint)rand.Next(1 << 30);
    }

    private static uint well512()
    {
        uint a, b, c, d;
        a = state[index];
        c = state[index + 13 & 15];
        b = a ^ c ^ (a << 16) ^ (c << 15);
        c = state[index + 9 & 15];
        c ^= c >> 11;
        a = state[index] = b ^ c;
        d = a ^ ((a << 5) & 0xda442d24u);
        index = index + 15 & 15;
        a = state[index];
        state[index] = a ^ b ^ d ^ (a << 2) ^ (b << 18) ^ (c << 28);
        return state[index];
    }

    private static ulong buf = 0;              // 64 random bits buffer
    private static int av = 0;                 // available bits in buf
    private static uint g_u = 0;               // copy of u
    private static int g_nd = 0;               // nd (needed) bits = highest set bit of u
    private static int g_s = 0;                // shift
    private static bool ones = false;          // true if u = 10, 110, 1110, ...

    private static uint well_uint(uint u)
    {
        uint x; int s;
        if (u == 0) return 0;
        if (u != g_u)
        {
            g_u = u;
            g_nd = nob(u);
            g_s = 64 - g_nd;
            ones = (u + 1 & u + 2) == 0;
        }
    L0: if (av < g_nd)
        {
            av += 32;
            buf <<= 32;
            buf |= well512();
        }
        x = (uint)buf << g_s >> g_s;
        buf >>= g_nd;
        av -= g_nd;
        if (x <= u) return x;
        if (ones) goto L0;                     // faster when u = 2, 6, 14, ...
        s = nob(x ^ u) - 1;
        av += s;
        buf <<= s;
        buf ^= x;
        goto L0;
    }

    private static int nob(uint u) // number of bits, position of the highest set bit, u > 0  
    {
        return
        u < 1u << 16 ? u < 1u << 08 ? u < 1u << 04 ? u < 1u << 02 ? u < 1u << 01 ? 01 : 02 :
                                                                    u < 1u << 03 ? 03 : 04 :
                                                     u < 1u << 06 ? u < 1u << 05 ? 05 : 06 :
                                                                    u < 1u << 07 ? 07 : 08 :
                                      u < 1u << 12 ? u < 1u << 10 ? u < 1u << 09 ? 09 : 10 :
                                                                    u < 1u << 11 ? 11 : 12 :
                                                     u < 1u << 14 ? u < 1u << 13 ? 13 : 14 :
                                                                    u < 1u << 15 ? 15 : 16 :
                       u < 1u << 24 ? u < 1u << 20 ? u < 1u << 18 ? u < 1u << 17 ? 17 : 18 :
                                                                    u < 1u << 19 ? 19 : 20 :
                                                     u < 1u << 22 ? u < 1u << 21 ? 21 : 22 :
                                                                    u < 1u << 23 ? 23 : 24 :
                                      u < 1u << 28 ? u < 1u << 26 ? u < 1u << 25 ? 25 : 26 :
                                                                    u < 1u << 27 ? 27 : 28 :
                                                     u < 1u << 30 ? u < 1u << 29 ? 29 : 30 :
                                                                    u < 1u << 31 ? 31 : 32;
    }

    static void Main()
    {
        random_bitmap();
        well_bitmap();
    }

    private static void random_bitmap()
    {
        int xmax = 450;
        int ymax = 450;
        var bmp = new Bitmap(xmax + 1, ymax + 1);
        int x, y, i = 0; uint z;
        for (x = 0; x <= xmax; x++)
        {
            for (y = 0; y < ymax; y++)
            {
                var rand = new Random(i++);
                z = (uint)rand.Next(1 << 30);
                if ((z & 1) == 0) bmp.SetPixel(x, y, Color.White);
                else bmp.SetPixel(x, y, Color.Black);
            }
        }
        bmp.Save(xmax.ToString() + "R.bmp");
    }

    private static void well_bitmap()
    {
        int xmax = 450;
        int ymax = 450;
        var bmp = new Bitmap(xmax + 1, ymax + 1);
        int x, y, i = 0; uint z;
        for (x = 0; x <= xmax; x++)
        {
            for (y = 0; y < ymax; y++)
            {
                init_well(i);
                i += 1;
                //z = well512();
                z = well_uint(1 << 30);
                z >>= 0;
                if ((z & 1) == 0) bmp.SetPixel(x, y, Color.White);
                else bmp.SetPixel(x, y, Color.Black);
            }
        }
        bmp.Save(xmax.ToString() + "W.bmp");
    }
}

No comments:

Post a Comment