Pagina's

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

No comments:

Post a Comment