Pagina's

2014/01/23

The kth Permutation

/*

It is based on:

    "Using Permutations in .Net for Improved System Security"

Factoradics are used to find the kth permutation of order n.

*/

class Permutation
{
    private static uint[] kth_perm(uint k, uint n)  //  k=0,n=3 ==> {0,1,2}
    {                                               //  k=1,n=3 ==> {0,2,1}
        if (n == 0) return new uint[] { };
        uint[] perm = new uint[n++];
        uint[] fdic = new uint[n];                  //  in factoradic form, no factorials!
        uint i = 1, d, m = n - 1;
        do
        {
            d = k / i;
            fdic[m--] = k - d * i++;
            k = d;
        }
        while (m != ~0u);
        for (k = n-- - 2; k != ~0u; k--)
        {
            m = perm[k++] = fdic[k];
            for (i = k--; i < n; i++) if (perm[i] >= m) perm[i]++;
        }
        return perm;
    }

    private static uint[] next_perm(uint[] perm)    //  {0,1,2} ==> {0,2,1}
    {
        int len = perm.Length - 1;
        if (len <= 0) return perm;
        int i = len;
        while (i > 0 && perm[i--] <= perm[i]) ;
        uint swap = perm[i];
        int j = len;
        while (j > 0 && perm[j] <= swap) j--;
        if (j == 0)
            do
            {
                swap = perm[len];
                perm[len--] = perm[j];
                perm[j++] = swap;
            }
            while (j < len);
        else
        {
            perm[i++] = perm[j];
            perm[j] = swap;
            while (len > i)
            {
                swap = perm[len];
                perm[len--] = perm[i];
                perm[i++] = swap;
            }
        }
        return perm;
    }

    static void Main()
    {
        uint[] perm = kth_perm(999999, 10);  //            {2,7,8,3,9,1,5,4,6,0}
        perm = next_perm(perm);              //            {2,7,8,3,9,1,5,6,0,4}
        perm = new uint[] { 3, 1, 3, 1 };    //  {3,1,3,1}
        perm = next_perm(perm);              //  {3,3,1,1}
        perm = next_perm(perm);              //  {1,1,3,3}
        perm = next_perm(perm);              //  {1,3,1,3}
        perm = next_perm(perm);              //  {1,3,3,1}
        perm = next_perm(perm);              //  {3,1,1,3}
        perm = next_perm(perm);              //  {3,1,3,1}
    }
}

2014/01/20

The mth Combination

/*
In "www.MWC58 Fast RNG" appeared an off topic item: 
The mth 2-element subset of the set { 0, 1 , 2 , 3 , ... , n } , or

    "Generating the mth Lexicographical Element of a Mathematical Combination"

The author, Dr. James McCaffrey, remarks: "Huh?" Also by McCaffrey:

    "Using Permutations in .Net for Improved Systems Security"

Factorials (factoradics) are used to find the mth permutation of order n.
It might be a good idea to read that first. A few years ago I published code 
based on it, sorry, the site disappeared, but: The kth Permutation
 
Binomial Coefficients (combinadics) are used to find the mth combination of order n.
McCaffrey observed: "when you map from an index value to a vector value, alternate number
representations (factoradics, combinadics, ...) are often the key to an elegant solution".
 
A few examples:

A set contains 1 item (n=1): "0" (it will become clear that it is easier to work zero
based, like indices in an array), so we can choose 1 subset (k=1) , the set itself: "0"
With n=2 , the set is: {0,1} and k=1 , the first (m=0) subset is: "0" , the second one
(m=1) : "1" . With k=2 , 1 subset, the set itself.
    
     n=1                n=2          n=2          n=3           n=3          n=3 
     k=1                k=1          k=2          k=1           k=2          k=3
    s=set={0}          s={0,1}      s={0,1}      s={0,1,2}     s={0,1,2}    s={0,1,2}
    m  ss(=subset)     m  ss        m  ss        m  ss         m  ss        m  ss
    0  0               0  0         0  0,1       0  0          0  0,1       0  0,1,2
    1  0               1  1         1  0,1       1  1          1  0,2       1  0,1,2
    ....               2  0         ......       2  2          2  1,2       ........
                       ....                      3  0          3  0,1       
                                                 ....          ......        

The number of possible combinations: " n! / k! / (n-k)! " = Choose(n,k) = bnm(n,k)

         
Another combination algorithm:

    Gosper's hack (item 175) 

A description: A Novel Application (snoob, page 14)    

A remark from the author, Henry S. Warren: ... the problem of finding the next higher number
after a given number that has the same number of 1-bits. You are forgiven if you are asking,
"Why on earth would anyone want to compute that?" The "\\\\////" image above is made with it,
turn it 90 degrees, spot the pattern.
*/

class Combination
{
    private static byte snoob(byte b)
    {
        byte smallest, ripple, ones;              //   x = 00001010  {1,3}
        smallest = (byte)(b & (~b + 1));          //       00000010
        ripple = (byte)(b + smallest);            //       00001100
        ones = (byte)(b ^ ripple);                //       00000110
        ones = (byte)((ones >> 2) / smallest);    //       00000000
        return (byte)(ripple | ones);             //       00001100  {2,3}
    }

    private static uint[] comb = new uint[] { };
    private static uint g_n = 0, g_k = 0;

    private static void next_comb(uint n, uint k)
    {
        if (k == 0 || k > n) comb = new uint[] { };
        else
        {
            if (n != g_n || k != g_k)
            {
                g_n = n; g_k = k; comb = new uint[k];
                for (n = 1; n < k; ) comb[n] = n++;
            }
            else
            {
                if (comb[0] == n - k) for (n = 0; n < k; ) comb[n] = n++;
                else
                {
                    k--; n--; while (comb[k--] == n--) ;
                    k++; n = ++comb[k]; while (++k < g_k) comb[k] = ++n;
                    n = g_n;
                }
            }
        }
    }

    private static void next_comb() { next_comb(g_n, g_k); }

    private static uint[] mth_comb(uint m, uint n, uint k)
    {
        if (k == 0 || k > n) return new uint[] { };
        uint i = 0, nk = n - k;
        g_n = n - 1; g_k = k - 1;
        comb = new uint[k];
        uint c = bnm(n, k);
        m = c - m % c - 1;
        c = c * nk-- / n--;
        while (i < g_k)
        {
            while (c > m) c = c * nk-- / n--;  //  while (c > m) c = c * nk-- / n--;
            m -= c;                            //  McCaffrey: A binary search approach for
            comb[i++] = g_n - n;               //  large values of n could improve performance.
            c = c * k-- / n--;                 //      www.bitlength factorial  
        }                                      //  ==> approximation bitlength binomial,
        while (c > m) c = c * nk-- / n--;      //  ==> a way to narrow the search ?
        comb[i] = g_n - n;
        g_n++; g_k++;
        return comb;
    }

    private static uint bnm(uint n, uint k)
    {
        if (k > n) return 0;
        if (k > n / 2) k = n - k;
        uint nk = n - k, c = k = 1;
        while (n > nk) c = c * n-- / k++;
        return c;
    }

    static void Main()
    {
        mth_comb(5, 5, 3);  //  { 0, 3 , 4 }
        next_comb();        //  { 1, 2 , 3 }
        next_comb(4, 2);    //  { 0, 1 }
        next_comb();        //  { 0, 2 }
    }

    private static ulong bnm64(uint n, uint k)  //  Binomial Coefficient needs an update
    {
        checked
        {
            if (k > n) return 0;
            if (k > n / 2) k = n - k;
            uint nk = n - k;
            uint c = k = 1;
            ulong cc;
            try
            {
                while (n > nk) c = c * n-- / k++;
                return c;
            }
            catch
            {
                try
                {
                    n++;
                    cc = c;
                    while (n > nk) cc = cc * n-- / k++;
                    return cc;
                }
                catch
                {
                    unchecked
                    {
                        return 0;
                        // using System.Numerics.BigInteger;
                        // n++; 
                        // BigInteger C = cc; 
                        // while (n > nk) C = C * n-- / k++;
                        // return C;
                    }
                }
            }
        }
    }
}

2014/01/06

CMR63 Fast Scaled Random Number Generator part 5

/*
"CMR" stands for "Constant Multiply Rotate", 
"63" stands for the period: > 2^63 (almost 2^64).
It's  based on a publication by Mark A. Overton(2011):

    http://www. Fast, High-Quality, Parallel Random Number Generators
 
CMR63:

    private static uint cmr()
    {
        z0 *= m0;
        z1 *= m1;
        z0 = z0 << s0 | z0 >> r0;
        z1 = z1 << s1 | z1 >> r1;
        return z0 ^ z1;
    }

The constants can be found in: http://www.OvertonTables.zip (cmr-32.txt) 

A small part of the cmr-32 table:

    cmr-32:  x = Constant*x; x = rotl(x,RL); 
    Seed value can be 1, or SeedStart..(SeedStart+SeedCnt-1). 
    SB = number of bits in seed (2^SB <= SeedCnt).
    SeedCnt  SeedStart  SB   Period    Constant  RL Prime factors of period 
    -------- ---------- -- ---------- ---------- -- ----------------------- 
    62973467 1377002680 25 4294966876 3563976171 16 2^2 1073741719 
    33525602 3202323436 24 4294965919 1422968075 16 307 757 18481 
    32022541 1180658772 24 4294966152 1977089609 19 2^3 3 41 1051 4153 
    31280659  554048116 24 4294966449  433149435 17 3 1259 1137137
     2749407 2613580375 21 4294950337  272690735 19 4294950337 (prime) 
     1095230 3119024045 20 4294928147   64333559 18 4294928147 (prime) 
      998895 3255944955 19 4294915769 3152644205 13 4294915769 (prime) 
      664885 3993266363 19 4294881427 4031235431 15 4294881427 (prime) 

The full table contains 5586 data rows, 272 periods are prime. 
The smallest period is 4284340603. Thus the smallest period of two
combined generators (with periods having no common factors) is at least:
    4284340603^2  
    4284340603^2 > 2^63.99  ==>  2^63.99 < period < 2^64
I use seed values to get unique pairs of constants (multipliers),
assuming random sequences generated by those unique pairs are not
correlated. Constants belonging to the 272 prime periods can be combined
with other constants (belonging to 272 of the 5586 - 272 non prime periods).
A more refined scheme might get the number of combinations closer to 5586/2. 
On my pc (Athlon II X4 640, 3GHz, 2GB, XP) 10^9 iterations of CMR63 take
6.27 seconds, so the full period takes about 3600 years! 
    2^63.99 / 10^9 * 6.27 / (60 * 60 * 24 * 365) ~ 3600
 
The scaled version:

"cmr_uint(u)" returns a value x ( 0 <= x <= u ).
The time to generate one scaled random number is 21 ns, mean 14 ns,
more accurate, mean maximum time: 21 ns, mean mean is mean: 14 ns.
Lots of random bits are spoiled, it does not matter: 3600 years! 
 
CMR16:  
 
Uses rotation (shift) value 16. Periods of the sub generators
( z2 = m2 * ... , z3 = m3 * ... ) are prime.
Result: 22 Divergent High Quality Random Number Sequences.
*/

class cmr63
{
    private static int s0 = 16, r0 = 16, s1 = 15, r1 = 17; // init_cmr(0);
    private static uint m0 = 3563976171, z0 = 4125873261,
                        m1 = 4031235431, z1 = 3803445283;

    private static void init_cmr(uint seed)
    {
        byte[] s = { 16, 16, 19, 17, 19, 18, 13, 15 };
        uint[] m = { 3563976171, 1422968075, 1977089609,  433149435,
                      272690735,   64333559, 3152644205, 4031235431};
        seed &= 3; m0 = m[seed]; s0 = s[seed]; r0 = 32 - s0; z0 = 1;
        seed ^= 7; m1 = m[seed]; s1 = s[seed]; r1 = 32 - s1; z1 = 1;
        cmr();
    }

    private static uint cmr()
    {
        z0 *= m0;
        z1 *= m1;
        z0 = z0 << s0 | z0 >> r0;
        z1 = z1 << s1 | z1 >> r1;
        return z0 ^ z1;
    }

    private static uint g_u = 0;
    private static int g_s = 0;

    private static uint cmr_uint(uint u)
    {
        uint x;
        if (u == 0) return 0;
        if (u != g_u)
        {
            g_u = u;
            g_s = 32 - nob(u);
        }
        do
        {
            z0 *= m0;
            z1 *= m1;
            z0 = z0 << s0 | z0 >> r0;
            z1 = z1 << s1 | z1 >> r1;
            x = (z0 ^ z1) >> g_s;
        }
        while (x > u);
        return x;
    }

    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()
    {
        // http://www.TestU01.BigCrush ?
    }

    private static uint m2 = 3745979853, z2 = 500031303,
                        m3 = 0623716905, z3 = 707339565;

    private static void init_cmr16(uint seed)
    {
        uint[] m = { 3745979853, 4055716687, 3693386591, 3542220329,
                     1775851103, 1866916287, 4188393139, 4141129223,
                     1173908643, 3198474053,   11119693, 1282266473,
                     4076777453, 3908725387, 3293562383, 2492630213,
                     1818407027,  608828557,  872259061, 2075607481,
                     1573125557, 2615661665, 1402711077, 3212405133,
                      680154359, 2023590663, 3458456891, 4184846215,
                     2408125305, 2558924297, 3008413683,  466035855,
                     1647905439, 2930730743,  733571709, 3997625831,
                     1919196763, 3392242035,  100431167,  579587817,
                     3074845609, 1931914705, 3131462569,  623716905};
        seed %= 22;
        m2 = m[seed];
        m3 = m[43 - seed];
        z2 = z3 = 1;
        cmr16();
    }

    private static uint cmr16()  //  2^63.99 < period < 2^64
    {
        z2 *= m2;
        z3 *= m3;
        z2 = z2 << 16 | z2 >> 16;
        z3 = z3 << 16 | z3 >> 16;
        return z2 ^ z3;
    }
}