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

Fast Cube Root ( < 2^32 )

/*                                                                                                            **
 The time to take a cube root from a 32 bits uint was: < 26 ns (mean: < 22 ns)
                                             new time: < 14 ns (mean: < 13 ns)
 
 The trick used for small square roots is used again: "unroll the loop".
 First step: A binary search, the "if" statements, large values ( >= 1u << 24 )
 are found relatively faster, small values ( < 64 ) are handled during the search.
 Second step: A jump into the unrolled loop, the "labels".
 
      |--------------------------------| 
      |         cro32(x)               |                         
      |-----|-----||------------|------|
      |   x |  ns ||          x |   ns |
      |-----|-----||------------|------|
      |   0 | 2,4 ||        511 |  4,6 |
      |   1 | 2,3 ||       1023 |  5,3 |
      |   3 | 2,3 ||       2047 |  5,3 |
      |   7 | 2,3 ||       4095 |  5,0 |
      |  15 | 3,6 ||      32768 |  7,2 |
      |  31 | 3,0 ||    1048576 |  8,5 |
      |  63 | 3,0 ||  268435456 | 12,0 |
      |  64 | 5,0 ||  536870912 | 12,3 |     heaviest case,
      |  65 | 5,0 || 1070599167 | 13,7 | <== 1070599167 is the
      | 100 | 5,0 || 1070599168 | 13,7 |     3th power of 1023, 
      | 127 | 5,0 || 1073741824 | 12,3 |     binary: 1111111111
      | 128 | 5,0 || 1073774592 | 12,0 |
      | 255 | 4,6 || 2147483648 | 12,3 |
      | 256 | 4,6 || 4294967295 | 12,0 |
      |-----|-----||------------|------|

       4294967295 roots in 53884 ms, mean time 12,55 ns, faster than square root!
         UPDATE 2013/11/04 50686 ms, mean time 11,80 ns, see: part 2
 
 ** Original version: Figure 2: Architecture of the iteration (alternative A)
    A digit-by-digit algorithm for radix-2 cube root and its implementation (2004)
    J.-A.Pineiro, J.D.Bruguera, L.Ciminiera, P.Montuschi
  
*/
using System;
using System.Diagnostics;

class cro_32
{
    private static uint cro32(uint x)
    {
        uint y = 4u, z = 16u, b = 0u;
        if (x < 1u << 24)
            if (x < 1u << 12)
                if (x < 1u << 06)
                    if (x < 1u << 03)
                        return x == 0u ? 0u : 1u;
                    else
                        return x < 27u ? 2u : 3u;
                else
                    if (x < 1u << 09) goto L8; else goto L7;
            else
                if (x < 1u << 18)
                    if (x < 1u << 15) goto L6; else goto L5;
                else
                    if (x < 1u << 21) goto L4; else goto L3;
        else
            if (x < 1u << 30)
                if (x < 1u << 27) goto L2; else goto L1;

        if (x >= 27u << 27) { x -= 27u << 27; z = 36u; y = 6u; } else { x -= 1u << 30; } goto M0;
    L1: if (x >= 27u << 24) { x -= 27u << 24; z = 36u; y = 6u; } else { x -= 1u << 27; } goto M1;
    L2: if (x >= 27u << 21) { x -= 27u << 21; z = 36u; y = 6u; } else { x -= 1u << 24; } goto M2;
    L3: if (x >= 27u << 18) { x -= 27u << 18; z = 36u; y = 6u; } else { x -= 1u << 21; } goto M3;
    L4: if (x >= 27u << 15) { x -= 27u << 15; z = 36u; y = 6u; } else { x -= 1u << 18; } goto M4;
    L5: if (x >= 27u << 12) { x -= 27u << 12; z = 36u; y = 6u; } else { x -= 1u << 15; } goto M5;
    L6: if (x >= 27u << 09) { x -= 27u << 09; z = 36u; y = 6u; } else { x -= 1u << 12; } goto M6;
    L7: if (x >= 27u << 06) { x -= 27u << 06; z = 36u; y = 6u; } else { x -= 1u << 09; } goto M7;
    L8: if (x >= 27u << 03) { x -= 27u << 03; z = 36u; y = 6u; } else { x -= 1u << 06; } goto M8;

    M0: if (x >= 61u << 24) { x -= 61u << 24; z += y * 2 + 1u; y += 1u; } y *= 2; z *= 4;
    M1: b = (y + z) * 3 + 1u << 21; if (x >= b) { x -= b; z += y * 2 + 1u; y += 1u; } y *= 2; z *= 4;
    M2: b = (y + z) * 3 + 1u << 18; if (x >= b) { x -= b; z += y * 2 + 1u; y += 1u; } y *= 2; z *= 4;
    M3: b = (y + z) * 3 + 1u << 15; if (x >= b) { x -= b; z += y * 2 + 1u; y += 1u; } y *= 2; z *= 4;
    M4: b = (y + z) * 3 + 1u << 12; if (x >= b) { x -= b; z += y * 2 + 1u; y += 1u; } y *= 2; z *= 4;
    M5: b = (y + z) * 3 + 1u << 09; if (x >= b) { x -= b; z += y * 2 + 1u; y += 1u; } y *= 2; z *= 4;
    M6: b = (y + z) * 3 + 1u << 06; if (x >= b) { x -= b; z += y * 2 + 1u; y += 1u; } y *= 2; z *= 4;
    M7: b = (y + z) * 3 + 1u << 03; if (x >= b) { x -= b; z += y * 2 + 1u; y += 1u; } y *= 2; z *= 4;
    M8: return x <= (y + z) * 3 ? y : y + 1u;
    }

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        cro32(225);
        cro32(1070599167);
        check_perfect_cubes32();
        time_cro32();
        cro32_all();
        Console.ReadLine();
    }

    private static void check_perfect_cubes32()
    {
        uint n = 0, a = 0, b = 1, c = 6;
        while (n <= 1625)
        {
            if (n != cro32(a)) Console.WriteLine("WRONG");
            a += b; b += c; c += 6; n++;
        }
        Console.WriteLine("CHECKED PERFECT CUBES");
    }

    private static void time_cro32()
    {
        int i, j; uint x; double t;
        uint[] a = { 0, 1, 3, 7, 15, 31, 63, 64, 65,
                     100, 127, 128, 255, 256, 511, 1023, 2047, 4095,
                     1u << 15, 1u << 20, 1u << 28 ,1u << 29 , 1070599167,
                     1070599168,1u << 30,(1u << 30) + (1u << 15), 1u << 31, ~0u };
        Console.WriteLine("|------------|------|");
        Console.WriteLine("|          x |  ns  |");
        Console.WriteLine("|------------|------|");
        for (j = 0; j < a.Length; j++)
        {
            x = a[j];
            cro32(x);
            sw.Restart();
            for (i = 0; i < 10000000; i++) cro32(x);
            sw.Stop();
            t = sw.ElapsedMilliseconds;
            sw.Restart();
            for (i = 0; i < 10000000; i++) ; // nada
            sw.Stop();
            t -= sw.ElapsedMilliseconds;
            Console.WriteLine("| {0,10} | {1:00.0} |", x, t / 10);
        }
        Console.WriteLine("|------------|------|");
        Console.WriteLine();
    }

    private static void cro32_all()
    {
        uint x; double t;
        cro32(~0u);
        sw.Restart();
        for (x = 0; x < ~0u; x++) cro32(x);
        sw.Stop();
        t = sw.ElapsedMilliseconds;
        sw.Restart();
        for (x = 0; x < ~0u; x++) ; // nada
        sw.Stop();
        t -= sw.ElapsedMilliseconds;
        Console.Write(x + " roots in " + t + " ms, ");
        Console.WriteLine("mean time {0:.00} ns", t * 1000000 / x);
        // 4294967295 roots in 53884 ms, mean time 12,55 ns
    }
}

Fast Square Root ( < 2^32 )

/*
The time to take a square root from a 32 bits uint was: ~35 ns, new time: <16 ns (mean: <14 ns).
 
It's based on: A square root hardware algorithm from Guy L. Steele,Jr. 

The slightly optimized C# version below is nearly twice as fast. A few "goto Label" statements
are used, they are fairly easy to remove, but I doubt it will be any faster (at all).
  
             |--------------------------------|
             |        Times "cro32(x)"        |
             |      Athlon X4, XP, 2 GB       |        
             |-----|-----||------------|------|
             |   x |  ns ||          x |   ns |
             |-----|-----||------------|------|
             |   0 | 1,6 ||        511 |  4,3 |
             |   1 | 1,6 ||       1023 |  4,6 |
             |   3 | 1,6 ||       2047 |  5,0 |
             |   7 | 2,3 ||       4095 |  5,6 |           4294967295 roots in 56962 ms, 
             |  15 | 2,3 ||      32768 |  6,6 |                     mean time 13,26 ns.
             |  31 | 3,6 ||    1048576 | 10,0 |
             |  63 | 3,0 ||  268435456 | 13,7 |
             |  64 | 4,3 ||  536870912 | 11,4 |
             |  65 | 4,3 || 1070599167 | 13,0 |
             | 100 | 5,0 || 1073741824 | 15,1 |
             | 127 | 4,3 || 1073741824 | 15,1 |
             | 128 | 4,3 || 1073774592 | 15,1 |
             | 255 | 4,0 || 2147483648 | 13,1 |
             | 256 | 4,7 || 4294967295 | 15,7 |
             |-----|-----||------------|------|

*/
using System;                                                     
using System.Diagnostics;                                         
                                                                  
class Square_Root                                                 
{                                                                 
    private static uint sro32(uint x)
    {
        uint y, z;
        if (x < 1u << 16)
            if (x < 1u << 08)
                if (x < 1u << 04) return x < 1u << 02 ? x + 3u >> 2 : x + 15u >> 3;
                else
                    if (x < 1u << 06)
                    { y = 1u << 03; x -= 1u << 04; if (x >= 5u << 02) { x -= 5u << 02; y |= 1u << 02; } goto L0; }
                    else
                    { y = 1u << 05; x -= 1u << 06; if (x >= 5u << 04) { x -= 5u << 04; y |= 1u << 04; } goto L1; }
            else                                             // slower (on my pc): .... y = 3u << 04; } goto L1; }
                if (x < 1u << 12)
                    if (x < 1u << 10)
                    { y = 1u << 07; x -= 1u << 08; if (x >= 5u << 06) { x -= 5u << 06; y |= 1u << 06; } goto L2; }
                    else
                    { y = 1u << 09; x -= 1u << 10; if (x >= 5u << 08) { x -= 5u << 08; y |= 1u << 08; } goto L3; }
                else
                    if (x < 1u << 14)
                    { y = 1u << 11; x -= 1u << 12; if (x >= 5u << 10) { x -= 5u << 10; y |= 1u << 10; } goto L4; }
                    else
                    { y = 1u << 13; x -= 1u << 14; if (x >= 5u << 12) { x -= 5u << 12; y |= 1u << 12; } goto L5; }
        else
            if (x < 1u << 24)
                if (x < 1u << 20)
                    if (x < 1u << 18)
                    { y = 1u << 15; x -= 1u << 16; if (x >= 5u << 14) { x -= 5u << 14; y |= 1u << 14; } goto L6; }
                    else
                    { y = 1u << 17; x -= 1u << 18; if (x >= 5u << 16) { x -= 5u << 16; y |= 1u << 16; } goto L7; }
                else
                    if (x < 1u << 22)
                    { y = 1u << 19; x -= 1u << 20; if (x >= 5u << 18) { x -= 5u << 18; y |= 1u << 18; } goto L8; }
                    else
                    { y = 1u << 21; x -= 1u << 22; if (x >= 5u << 20) { x -= 5u << 20; y |= 1u << 20; } goto L9; }
            else
                if (x < 1u << 28)
                    if (x < 1u << 26)
                    { y = 1u << 23; x -= 1u << 24; if (x >= 5u << 22) { x -= 5u << 22; y |= 1u << 22; } goto La; }
                    else
                    { y = 1u << 25; x -= 1u << 26; if (x >= 5u << 24) { x -= 5u << 24; y |= 1u << 24; } goto Lb; }
                else
                    if (x < 1u << 30)
                    { y = 1u << 27; x -= 1u << 28; if (x >= 5u << 26) { x -= 5u << 26; y |= 1u << 26; } goto Lc; }
                    else
                    { y = 1u << 29; x -= 1u << 30; if (x >= 5u << 28) { x -= 5u << 28; y |= 1u << 28; } }

        z = y | 1u << 26; y /= 2; if (x >= z) { x -= z; y |= 1u << 26; }
    Lc: z = y | 1u << 24; y /= 2; if (x >= z) { x -= z; y |= 1u << 24; }
    Lb: z = y | 1u << 22; y /= 2; if (x >= z) { x -= z; y |= 1u << 22; }
    La: z = y | 1u << 20; y /= 2; if (x >= z) { x -= z; y |= 1u << 20; }
    L9: z = y | 1u << 18; y /= 2; if (x >= z) { x -= z; y |= 1u << 18; }
    L8: z = y | 1u << 16; y /= 2; if (x >= z) { x -= z; y |= 1u << 16; }
    L7: z = y | 1u << 14; y /= 2; if (x >= z) { x -= z; y |= 1u << 14; }
    L6: z = y | 1u << 12; y /= 2; if (x >= z) { x -= z; y |= 1u << 12; }
    L5: z = y | 1u << 10; y /= 2; if (x >= z) { x -= z; y |= 1u << 10; }
    L4: z = y | 1u << 08; y /= 2; if (x >= z) { x -= z; y |= 1u << 08; }
    L3: z = y | 1u << 06; y /= 2; if (x >= z) { x -= z; y |= 1u << 06; }
    L2: z = y | 1u << 04; y /= 2; if (x >= z) { x -= z; y |= 1u << 04; }
    L1: z = y | 1u << 02; y /= 2; if (x >= z) { x -= z; y |= 1u << 02; }
    L0: return x > y ? y / 2 | 1u : y / 2;
    }

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        Console.WriteLine(sro32(4));
        Console.WriteLine(sro32(40));
        Console.WriteLine(sro32(400));
        Console.WriteLine(sro32(4000));
        Console.WriteLine(sro32(40000));
        Console.WriteLine(sro32(400000));
        Console.WriteLine(sro32(4000000));
        Console.WriteLine(sro32(40000000));
        Console.WriteLine(sro32(400000000));
        Console.WriteLine(sro32(4000000000));
        test_perfect_squares();
        time_sro32();
        sro32_all(); // ~1 minute                                 
        Console.ReadLine();
    }
    private static void test_perfect_squares()
    {
        uint n = 0, a = 0, b = 1;
        do
        {
            if (sro32(a) != n) Console.WriteLine("WRONG");
            a += b;
            b += 2;
        }
        while (n++ < 65535);
        Console.WriteLine("Perfect squares tested");
    }
    private static void time_sro32()
    {
        int i, j; uint x; double t;
        //uint[] a = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 };
        uint[] a = { 0, 1, 3, 7, 15, 31, 63, 64, 65,
                     100, 127, 128, 255, 256, 511, 1023, 2047, 4095,
                     1u << 15, 1u << 20, 1u << 28 ,1u << 29 , 1070599167,
                     1u << 30,1u << 30,(1u << 30) + (1u << 15), 1u << 31, 0xfffe0001u, ~0u };
        Console.WriteLine("|------------|------|");
        Console.WriteLine("|          x |  ns  |");
        Console.WriteLine("|------------|------|");
        for (j = 0; j < a.Length; j++)
        {
            x = a[j];
            sro32(x);
            sw.Restart();
            for (i = 0; i < 10000000; i++) sro32(x);
            sw.Stop();
            t = sw.ElapsedMilliseconds;
            sw.Restart();
            for (i = 0; i < 10000000; i++) ; // nada
            sw.Stop();
            t -= sw.ElapsedMilliseconds;
            Console.WriteLine("| {0,10} | {1:00.0} |", x, t / 10);
        }
        Console.WriteLine("|------------|------|");
        Console.WriteLine();
    }
    private static void sro32_all()
    {
        uint x; double t;
        sro32(~0u);
        sw.Restart();
        for (x = 0; x < ~0u; x++) sro32(x);
        sw.Stop();
        t = sw.ElapsedMilliseconds;
        sw.Restart();
        for (x = 0; x < ~0u; x++) ; // nada
        sw.Stop();
        t -= sw.ElapsedMilliseconds;
        Console.Write(x + " roots in " + t + " ms, ");
        Console.WriteLine("mean time {0:.00} ns",
                          (double)(t) * 1000000 / x);
        // 4294967295 roots in 56962 ms, mean time 13,26 ns
    }
}

//------------------------------------------------------------
// Another one, slower, the line "t = ..." is done many times.
//
// private static uint sro32s(uint x)
//     uint y1 = 0, y2 = 0, t; int k = 15;
//     for (; k >= 0; k--)
//     {
//         t = y2 + (y1 << k + 1) | (1u << k * 2);
//         if (x >= t) { y2 = t; y1 |= 1u << k; }
//     }
//     return y1;
//-----------------------------------------------------------

Small Cube Roots ( < 2^64 )

/*
The time to take a cube root from a 64 bits ulong was: ~330 ns, new time: < 75 ns
The time to take a cube root from a 32 bits  uint was:  ~25 ns, new time: < 22 ns (mean)
* UPDATE (2013-09-22) new time for 64 bits: < 70 ns, casting float to int is faster than 
                                                     casting float to uint.

The 64 bits version: 
 
    www.hackersdelight.org/hdcodetxt/acbrt.c.txt
 
    Why? 0x2a5137a0

A union data type is used, not available in C#, first work around: "BitConverter", at the
bottom: "cro64old(x)". It takes ~110 ns, "BitConverter" takes ~10 ns, it's used four times,
so 40 ns are gone. A better solution:
 
    Unions (or an equivalent) in C# - Sairama's Tip of the Day

A 63 bits version might be slightly faster, no overflow checks, less correction loops?
 
 
The 32 bits version:
 
A while loop is replaced by a "do while", leading (triple) zero bits are counted to avoid
useless iterations, all (4294967295) roots in 92,7 s, mean time 21,6 ns.
** UPDATE (2013-09-26) new time for 32 bits: < 15 ns, see: Fast Cube Root 
 
  
    |-------------||--------------------------||-------------------------------------------| 
    |  cro12(x)   ||       cr032(x)           ||               cro64(random x)             |
    |------|------||------------|-------------||------|----------------------|-------------|
    |    x |  ns  ||          x |     ns   ** || bits |                    x |     ns   *  |  
    |------|------||------------|------|------||------|----------------------|------|------|  
    |    0 |  3,3 ||          0 |  5,8 |  2,4 ||   64 | 13397031575470575331 | 74,7 | 69,8 |  
    |    1 |  3,7 ||          1 |  6,5 |  2,3 ||   63 |  6512002244030281881 | 70,4 | 66,8 |  
    |    3 |  5,2 ||          3 |  6,5 |  2,3 ||   62 |  4508122487018397160 | 73,6 | 68,5 |  
    |    7 |  5,3 ||          7 |  6,5 |  2,3 ||   61 |  1458795770683441400 | 74,0 | 66,8 |  
    |   15 |  6,3 ||         15 |  7,7 |  3,6 ||   60 |   994257496981222664 | 74,4 | 66,8 |  
    |   31 |  7,3 ||         31 |  7,5 |  3,0 ||   59 |   470915095588726850 | 74,0 | 66,8 |  
    |   63 |  7,3 ||         63 |  7,5 |  3,0 ||   58 |   201745106158251712 | 74,0 | 66,8 |  
    |   64 |  7,3 ||         64 |  9,4 |  4,8 ||   57 |    78804375734046440 | 74,0 | 66,8 |  
    |   65 |  8,3 ||         65 |  9,4 |  4,6 ||   56 |    57442872892622787 | 74,0 | 66,8 |  
    |  100 |  8,3 ||        100 |  9,4 |  4,6 ||   55 |    19109171792964730 | 74,0 | 66,8 |  
    |  127 |  9,3 ||        127 |  9,4 |  4,6 ||   54 |     9099755032911121 | 74,0 | 66,8 |  
    |  128 |  9,3 ||        128 |  9,4 |  4,6 ||   53 |     8170442697918990 | 74,0 | 66,8 |  
    |  255 | 10,3 ||        255 |  9,4 |  5,0 ||   52 |     2913033827815339 | 74,0 | 66,8 |  
    |  256 | 10,3 ||        256 |  9,4 |  5,0 ||   51 |     1395743698554671 | 74,0 | 66,8 |  
    |  511 | 11,3 ||        511 |  9,6 |  5,0 ||   50 |      929340477894354 | 74,0 | 66,8 |  
    | 1023 | 14,3 ||       1023 | 10,7 |  5,0 ||   49 |      351768287627885 | 74,1 | 66,8 |  
    | 2047 | 16,3 ||       2047 | 10,5 |  5,7 ||   48 |      165110734465934 | 70,9 | 66,5 |  
    | 4095 | 23,3 ||       4095 | 11,5 |  5,3 ||   47 |      116364019662847 | 73,6 | 66,8 |  
    |      |      ||      32768 | 13,7 |  6,4 ||   46 |       45894645838798 | 70,9 | 66,5 |  
    |      |      ||    1048576 | 15,7 |  8,7 ||   45 |       25652614759428 | 70,9 | 65,5 |  
    |      |      ||  268435456 | 19,7 | 11,0 ||   44 |       12123210057895 | 71,7 | 65,5 |  
    |      |      ||  536870912 | 20,1 | 12,0 ||   43 |        5902555891371 | 71,7 | 65,5 |  
    |      |      || 1073741824 | 25,8 | 11,3 ||   42 |        3490450430912 | 71,7 | 65,5 |  
    |      |      || 1073741824 | 25,8 | 11,3 ||   41 |        1149630654634 | 71,7 | 65,5 |  
    |      |      || 1073774592 | 25,8 | 11,3 ||   40 |         902265632303 | 71,7 | 65,5 |  
    |      |      || 2147483648 | 21,1 | 12,3 ||   39 |         343393482192 | 71,7 | 65,5 |  
    |      |      || 4294967295 | 21,8 | 12,3 ||   38 |         198281228228 | 71,7 | 65,5 |  
    |      |      ||            |      |      ||   37 |          93252157892 | 71,7 | 65,5 |  
    |      |      ||            |      |      ||   36 |          41910484903 | 71,7 | 65,5 |
    |      |      ||            |      |      ||   35 |          23282355985 | 71,8 | 65,5 |
    |      |      ||            |      |      ||   34 |           8677323163 | 71,7 | 65,5 |
    |      |      ||            |      |      ||   33 |           6714386675 | 71,7 | 65,5 |
    |------|------||------------|------|------||------|----------------------|------|------|
 
*/
using System;
using System.Diagnostics;
using System.Runtime.InteropServices;

[StructLayout(LayoutKind.Explicit)]

public struct fu_32   // float <==> uint
{
    [FieldOffset(0)]
    public float f;
    [FieldOffset(0)]
    public uint u;
}
class cro_64_32
{
    private static uint cro64(ulong x)
    {
        if (x >= 18446724184312856125) return 2642245;
        float fx = (float)x;
        fu_32 fu32 = new fu_32();
        fu32.f = fx;
        uint uy = fu32.u / 4;
        uy += uy / 4;
        uy += uy / 16;
        uy += uy / 256;
        uy += 0x2a5137a0;
        fu32.u = uy;
        float fy = fu32.f;
        fy = 0.33333333f * (fx / (fy * fy) + 2.0f * fy);

        // uint y1 = (uint)
        //    (0.33333333f * (fx / (fy * fy) + 2.0f * fy));

        int y0 = (int)                                      // 2013-09-22
            (0.33333333f * (fx / (fy * fy) + 2.0f * fy));   // 5 ns 
        uint y1 = (uint)y0;                                 // faster

        ulong y2, y3;
        if (y1 >= 2642245)
        {
            y1 = 2642245;
            y2 = 6981458640025;
            y3 = 18446724184312856125;
        }
        else
        {
            y2 = (ulong)y1 * y1;
            y3 = y2 * y1;
        }
        if (y3 > x)
        {
            y1 -= 1;
            y2 -= 2 * y1 + 1;
            y3 -= 3 * y2 + 3 * y1 + 1;
            while (y3 > x)
            {
                y1 -= 1;
                y2 -= 2 * y1 + 1;
                y3 -= 3 * y2 + 3 * y1 + 1;
            }
            return y1;
        }
        do
        {
            y3 += 3 * y2 + 3 * y1 + 1;
            y2 += 2 * y1 + 1;
            y1 += 1;
        }
        while (y3 <= x);
        return y1 - 1;
    }

    private static uint cro32(uint x)
    {
        uint y = 0, z = 0, b = 0;
        int s = x < 1u << 24 ? x < 1u << 12 ? x < 1u << 06 ? x < 1u << 03 ? 00 : 03 :
                                                             x < 1u << 09 ? 06 : 09 :
                                              x < 1u << 18 ? x < 1u << 15 ? 12 : 15 :
                                                             x < 1u << 21 ? 18 : 21 :
                               x >= 1u << 30 ? 30 : x < 1u << 27 ? 24 : 27;
        do
        {
            y *= 2;
            z *= 4;
            b = 3 * y + 3 * z + 1 << s;
            if (x >= b)
            {
                x -= b;
                z += 2 * y + 1;
                y += 1;
            }
            s -= 3;
        }
        while (s >= 0);
        return y;
    }

    private static uint cro12(uint x)
    {
        uint y = 0, a = 0, b = 1, c = 0;
        while (a < x)
        {
            y++;
            b += c;
            a += b;
            c += 6;
        }
        if (a != x) y--;
        return y;
    }

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        cro64(long.MaxValue);
        check_perfect_cubes64();
        check_perfect_cubes32();
        time_cro64();
        time_cro32();
        time_cro12();
        cro32_all();              // ~90 s, sorry, no cro64_all() ;) it would take ~44 millennia.  
        Console.ReadLine();
    }
    private static void check_perfect_cubes64()
    {
        ulong n = 0, a = 0, b = 1, c = 6;
        while (n <= 2642245)
        {
            if (n != cro64(a)) Console.WriteLine("WRONG");
            a += b; b += c; c += 6; n++;
        }
        Console.WriteLine("CHECKED64");
    }
    private static void check_perfect_cubes32()
    {
        uint n = 0, a = 0, b = 1, c = 6;
        while (n <= 1625)
        {
            if (n != cro32(a)) Console.WriteLine("WRONG");
            a += b; b += c; c += 6; n++;
        }
        Console.WriteLine("CHECKED32");
    }

    private static Random random = new Random(0);       // Random(0)/Random(), same/random sequence 
    private static uint rnd32()                         // random uint
    {
        return (uint)(random.Next(1 << 16)) << 16 |
               (uint)(random.Next(1 << 16));
    }
    private static uint rnd32(uint n)                   // random uint of n bits
    {
        if (n < 2) return n;
        if (n > 32) n = 32;
        return (((uint)(random.Next(1 << 16)) << 16 |
                 (uint)(random.Next(1 << 16)))
                >> 32 - (int)n) | (1u << (int)n - 1);
    }
    private static ulong rnd64(uint n)                 // random ulong of n bits
    {
        if (n < 2) return n;
        if (n > 64) n = 64;
        ulong u = rnd32();
        u = u << 32 | rnd32();
        u >>= 64 - (int)n;
        u |= 1uL << (int)n - 1;
        return u;
    }

    private static void time_cro64()
    {
        int i; uint j; ulong x; double t;
        Console.WriteLine("|------|----------------------|------|");
        Console.WriteLine("| bits |                    x |  ns  |");
        Console.WriteLine("|------|----------------------|------|");
        Console.WriteLine("|             PRESS A KEY            |");
        Console.ReadLine();
        for (j = 64; j > 32; j -= 1)
        {
            x = rnd64(j);
            cro64(x);
            sw.Restart();
            for (i = 0; i < 10000000; i++) cro64(x);
            sw.Stop();
            t = sw.ElapsedMilliseconds;
            sw.Restart();
            for (i = 0; i < 10000000; i++) ; // nada
            sw.Stop();
            t -= sw.ElapsedMilliseconds;
            Console.WriteLine("|   {0,2} | {1,20} | {2:00.0} |", j, x, t / 10);
        }
        Console.WriteLine("|------|----------------------|------|");
        Console.WriteLine();
    }
    private static void time_cro32()
    {
        int i, j; uint x; double t;
        uint[] a = { 0, 1, 3, 7, 15, 31, 63, 64, 65,
                     100, 127, 128, 255, 256, 511, 1023, 2047, 4095,
                     1u << 15, 1u << 20, 1u << 28 ,1u << 29 , 1070599167,
                     1u << 30,1u << 30,(1u << 30) + (1u << 15), 1u << 31, ~0u };
        Console.WriteLine("|------------|------|");
        Console.WriteLine("|          x |  ns  |");
        Console.WriteLine("|------------|------|");
        for (j = 0; j < a.Length; j++)
        {
            x = a[j];
            cro32(x);
            sw.Restart();
            for (i = 0; i < 10000000; i++) cro32(x);
            sw.Stop();
            t = sw.ElapsedMilliseconds;
            sw.Restart();
            for (i = 0; i < 10000000; i++) ; // nada
            sw.Stop();
            t -= sw.ElapsedMilliseconds;
            Console.WriteLine("| {0,10} | {1:00.0} |", x, t / 10);
        }
        Console.WriteLine("|------------|------|");
        Console.WriteLine();
    }
    private static void time_cro12()
    {
        int i, j; uint x; double t;
        uint[] a = { 0, 1, 3, 7, 15, 31, 63, 64, 65,
                     100, 127, 128, 255, 256, 511, 1023, 2047, 4095 };
        Console.WriteLine("|------------|------|");
        Console.WriteLine("|          x |  ns  |");
        Console.WriteLine("|------------|------|");
        for (j = 0; j < a.Length; j++)
        {
            x = a[j];
            cro12(x);
            sw.Restart();
            for (i = 0; i < 10000000; i++) cro12(x);
            sw.Stop();
            t = sw.ElapsedMilliseconds;
            sw.Restart();
            for (i = 0; i < 10000000; i++) ; // nada
            sw.Stop();
            t -= sw.ElapsedMilliseconds;
            Console.WriteLine("| {0,10} | {1:00.0} |", x, t / 10);
        }
        Console.WriteLine("|------------|------|");
        Console.WriteLine();
    }
    private static void cro32_all()
    {
        uint x; double t;
        cro32(~0u);
        sw.Restart();
        for (x = 0; x < ~0u; x++) cro32(x);
        sw.Stop();
        t = sw.ElapsedMilliseconds;
        sw.Restart();
        for (x = 0; x < ~0u; x++) ; // nada
        sw.Stop();
        t -= sw.ElapsedMilliseconds;
        Console.Write(x + " roots in " + t + " ms, ");
        Console.WriteLine("mean time {0:.00} ns",
                          (double)(t) * 1000000 / x);
        // 4294967295 roots in 92732 ms, mean time 21,59 ns
    }
    private static uint cro64old(ulong x)
    {
        if (x >= 18446724184312856125) return 2642245;                       // avoid overflow
        float fx = (float)x;           //------------------------------------------------------\
        uint uy = BitConverter.ToUInt32(BitConverter.GetBytes(fx), 0) / 4;   // ~20 ns          \
        uy += uy / 4;                                                        //                  \
        uy += uy / 16;                                                       //                   \
        uy += uy / 256;                                                      //                    \
        uy += 0x2a5137a0;                                                    //                     \
        float fy = BitConverter.ToSingle(BitConverter.GetBytes(uy), 0);      // ~20 ns              32 bits
        fy = 0.33333333f * (fx / (fy * fy) + 2.0f * fy);                     //                     /
        uint y1 = (uint)(0.33333333f * (fx / (fy * fy) + 2.0f * fy));        //                    /
        ulong y2, y3;                                                        //                   /
        if (y1 >= 2642245)                                                   // avoid overflow   /
        {                                                                    //                 /
            y1 = 2642245;              //------------------------------------------------------/                                                  
            y2 = 6981458640025;
            y3 = 18446724184312856125;
        }
        else
        {
            y2 = (ulong)y1 * y1;
            y3 = y2 * y1;
        }
        if (y3 > x)
        {
            y1 -= 1;
            y2 -= 2 * y1 + 1;
            y3 -= 3 * y2 + 3 * y1 + 1;
            while (y3 > x)
            {
                y1 -= 1;
                y2 -= 2 * y1 + 1;
                y3 -= 3 * y2 + 3 * y1 + 1;
            }
            return y1;
        }
        do
        {
            y3 += 3 * y2 + 3 * y1 + 1;
            y2 += 2 * y1 + 1;
            y1 += 1;
        }
        while (y3 <= x);
        return y1 - 1;
    }
}