Pagina's

2013/10/19

Uniform pseudo random numbers in a range

/*

Random number generators (RNG) are divided in two categories,
hardware RNG, that provide “true” random numbers, and algorithmic RNG, that generate
pseudo random numbers (PRNG). Both types usually generate random numbers "X" as independent
uniform samples in a range 0, . . . 2^b − 1, with b = 8, 16, 32 or b = 64. In applications,
it is instead sometimes desirable to draw random numbers as independent uniform samples "U"
in a range " 1, . . . M ", where moreover "M" may change between drawings.
Transforming "X" to "U" , or many "X" to many "U" is sometimes known as scaling.

Some text is copied from other sources. You find them in "References".
I found it quite hard to explain some of the underlying ideas.
I am grateful to the original authors. Their words saved me a lot of time.

Limits:

    0 <= M <= 2^32 - 1 = uint.MaxValue  
    0 <= U <= M

UPDATE: part 2 

Example: 
Suppose a PRNG produces one random bit each time it is called, and we have: M = 5 ,
binary: M = 101. Call the generator three times to get three successive bits. 
Result: 0 <= U <= 7, binary: 000 <= U <= 111. If "U > M" , than we take the next three bits
from the generator, and repeat to do so until "U <= M".

Another way to find a valid value for "U".
-Take the first, the most significant, bit from the generator.
If it's zero, take two other bits from the generator,
and return a valid value for "U": 000 or 001 or 010 or 011.
If the first bit isn't zero, so it's one, 
take the second bit from the generator.
If the second bit is zero, take the third bit from the generator,
and return a valid value for "U": 100 or 101.
Finally: the first and the second bit are both set,
no need to take a third bit from the generator,
a random bit is saved, we have to try again. 
-Take the first, the most significant, bit from the generator...

              /  U = 0..  return 000 or 001 or 010 or 011
             /
         ---<    U = 10.  return 100 or 101
        |    \ 
        |     \  U = 11   --try again--
        |                              |
        |                              |
         ----<<--------<<--------<<---- 

The PRNG used: .Net 4.0 System.Random
It might have a bug: The period is not guaranteed to be longer than: 2^55-1
                     (I do have some doubts, but no proof)
 
Some important properties of System.Random:

    Random Constructors:

        Random()          Initializes a new instance of the Random class, using a time-dependent value.
                   
        Random(int seed)  Initializes a new instance of the Random class, using the specified seed.

    To improve performance, create one Random object to generate many random numbers over time,
    instead of repeatedly creating new Random objects to generate one random number, so:

        private static Random rand = new Random()

    Random Methods:

        rand.Next()                    Returns a positive random number, greater than or equal to zero,
                                       less than int.MaxValue

        rand.Next(int max)             Returns a positive random number, greater than or equal to zero,
                                       less than max ( max >= 0 ) , ........ rand.Next(0) returns zero 
 
        rand.Next(int min, int max)    Returns a positive random number, greater than or equal to min,
                                       less than max ( max >= min ) , ...... 

    rand.Next() is about twice as fast as rand.Next(int max) , 
    rand.Next(int max) is almost as fast as rand.Next(int min, int max),
    especially when you use something like: return min + rand.Next(max - min)

    It's impossible to create more than 30 random bits at once.
    To create those 30 bits, use: rand.Next(1 << 30)  

Various sources warn about the lower bits of numbers generated by PRNG's, they shouldn't be random.
One way to examine a PRNG is to create a visualisation of the numbers it produces.
Humans are really good at spotting patterns, and visualisation allows you to use your eyes and 
brain directly for this purpose. While you shouldn't consider this type of approach an exhaustive 
or formal analysis, it is a nice and quick way to get a rough impression of a given generator's
performance. The background of the image at the top looks pretty random, no recognizable patterns.
It was created with the code at the very bottom: "random_bitmap()"

Description of the code:

gen_nd_bits(int nd) returns the needed amount of random bits. It uses a 32 bits buffer.
The buffer is replenished when necessary.

gen_uint_fast(uint u) returns a random value between zero (incl.) and u (incl.)
It calls gen_nd_bits(int nd) . When the value is too large another call is done.

gen_uint_eff(uint u) uses random bits more efficiently.
It works bit by bit, from most to least significant bit. When it becomes clear that the (partial)
formed value is smaller than the upper limit ( uint u ) , remaining needed random bits are
added relatively fast. When the (partial) formed value is too large another call is done.

On my trustworthy PC the average (full range) results are:

    |-----------------------|---------------|
    |     function          | * 10^9 bits/s | 
    |-----------------------|---------------|
    |    rand.Next()        |      2.9      | <-- 1 CD in 2 seconds!
    |    rand.Next(1 << 30) |      1.5      |  
    |   gen_nd_bits(int nd) |      1.2      | 
    | gen_uint_fast(uint u) |      0.9      | 
    |  gen_uint_eff(uint u) |      0.7      |
    |-----------------------|---------------|

The efficiency, random bits used per bit: 

      gen_uint_fast(uint u)   1.386 bits/bit, remarkably close to ln(4)
       gen_uint_eff(uint u)   1.044 bits/bit

References:

    José's Daylight Dices 
    RANDOM.ORG 
    Ask Dr. Math, Probability and Random Numbers
    Bit recycling for scaling random number generators
    Math.NET Numerics
    MSDN
    Bug ? System.Random (january 2011)

*/

using System;
using System.Drawing;                          // add a ref!
class scaled_random
{
    private static Random rand = new Random();
    private static uint buf = 0;               // 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;               // copy of nd (needed) = highest set bit of u

    private static uint gen_nd_bits(int nd)    // generate nd bits, 0 <= nd <= 32
    {
        uint x;
        if (nd == 0) return 0;
        if (av <= 2)                           // buffer should contain at least 2 bits, other-
        {                                      // wise there may be too few bits when nd > 30 
            buf <<= 30;
            buf |= (uint)rand.Next(1 << 30);
            av += 30;
        }
        if (av < nd)
        {
            x = buf;
            buf = (uint)rand.Next(1 << 30);
            nd -= av;
            x <<= nd;
            x |= buf << 32 - nd >> 32 - nd;
            buf >>= nd;
            av = 30 - nd;
        }
        else
        {
            x = buf << 32 - nd >> 32 - nd;
            buf >>= nd;
            av -= nd;
        }
        return x;
    }

    private static uint gen_uint_fast(uint u)  // 0 <= returned value <= u  
    {
        uint x;
        if (u == 0) return 0;
        if (u != g_u)                          // if u changes, g_nd might change 
        {
            g_u = u;
            g_nd = nob(u);                     // needed: number of bits of u   
        }
    L0: x = gen_nd_bits(g_nd);
        if (x <= u) return x;
        goto L0;
    }

    private static bool ones = false;          // true if u is (binary) 1, 11, 111, 1111, ...
    private static uint gen_uint_eff(uint u)   // 0 <= returned value <= u  
    {
        uint x; int nd;
        if (u == 0) return 0;
        if (u != g_u)
        {
            g_u = u;
            g_nd = nob(u);
            ones = (u + 1 & u) == 0;
        }
        if (ones) return gen_nd_bits(g_nd);
    L0: nd = g_nd - 1;
        if (gen_nd_bits(1) == 0) return gen_nd_bits(nd);
        x = 1u << nd;
        nd--;
        while (nd > 0)
        {
            x |= gen_nd_bits(1) << nd;
            if (x >> nd < u >> nd) return x | gen_nd_bits(nd);
            if (x > u) goto L0;
            nd--;
        }
        x |= gen_nd_bits(1);
        if (x <= u) return x;
        goto L0;
    }

    private static int nob(uint u) // number of bits, position of the highest set bit, u > 0  
    {                              // afaik fastest way in C#
        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()
    {
        test_gen_uint_fast(17);
        test_gen_uint_eff(54);
        test_4_buckets_fast(123456789);
        test_4_buckets_eff(~0u);
        random_bitmap();
    }

    private static void test_gen_uint_fast(uint u)
    {
        int i; int[] c = new int[u + 1];
        for (i = 0; i < 1000000; i++) c[gen_uint_fast(u)]++;
        for (i = 0; i <= u; i++) Console.WriteLine("{0,3}{1,7}", i, c[i]);
        Console.ReadLine(); Console.Clear();
    }
    private static void test_gen_uint_eff(uint u)
    {
        int i; int[] c = new int[u + 1];
        for (i = 0; i < 1000000; i++) c[gen_uint_eff(u)]++;
        for (i = 0; i <= u; i++) Console.WriteLine("{0,3}{1,7}", i, c[i]);
        Console.ReadLine(); Console.Clear();
    }
    private static void test_4_buckets_fast(uint u) // rounding errors!
    {
        uint x; int i, c0 = 0, c1 = 0, c2 = 0, c3 = 0;
        for (i = 0; i < 1000000; i++)
        {
            x = gen_uint_fast(u);
            if (x < u / 4) c0++;
            else if (x < u / 2) c1++;
            else if (x < u / 2 + u / 4) c2++; else c3++;
        }
        Console.WriteLine(c0); Console.WriteLine(c1);
        Console.WriteLine(c2); Console.WriteLine(c3);
        Console.ReadLine(); Console.Clear();
    }
    private static void test_4_buckets_eff(uint u)
    {
        uint x; int i, c0 = 0, c1 = 0, c2 = 0, c3 = 0;
        for (i = 0; i < 1000000; i++)
        {
            x = gen_uint_eff(u);
            if (x < u / 4) c0++;
            else if (x < u / 2) c1++;
            else if (x < u / 2 + u / 4) c2++; else c3++;
        }
        Console.WriteLine(c0); Console.WriteLine(c1);
        Console.WriteLine(c2); Console.WriteLine(c3);
        Console.ReadLine(); Console.Clear();
    }

    private static void random_bitmap()
    {
        int xmax = 900;
        int ymax = 900;
        Bitmap bmp = new Bitmap(xmax + 1, ymax + 1);
        int x, y, z;
        for (x = 0; x <= xmax; x++)
        {
            for (y = 0; y <= ymax; y++)
            {
                z = 7 & rand.Next(1 << 30);
                switch (z)
                {
                    case 0: bmp.SetPixel(x, y, Color.Yellow); break;
                    case 1: bmp.SetPixel(x, y, Color.Orange); break;
                    case 2: bmp.SetPixel(x, y, Color.Purple); break;
                    case 3: bmp.SetPixel(x, y, Color.White); break;
                    case 4: bmp.SetPixel(x, y, Color.Green); break;
                    case 5: bmp.SetPixel(x, y, Color.Black); break;
                    case 6: bmp.SetPixel(x, y, Color.Blue); break;
                    case 7: bmp.SetPixel(x, y, Color.Red); break;
                }
            }
        }
        bmp.Save(xmax.ToString() + ".bmp");
    }
}