How can i calculate the Nth combo based only on it's index. There should be (n+k-1)!/(k!(n-1)!) combinations with repetitions.
with n=2, k=5 you get:
0|{0,0,0,0,0}
1|{0,0,0,0,1}
2|{0,0,0,1,1}
3|{0,0,1,1,1}
4|{0,1,1,1,1}
5|{1,1,1,1,1}
So black_magic_function(3) should produce {0,0,1,1,1}.
This will be going into a GPU shader, so i want each work-group/thread to be able to figure out their subset of permutations without having to store the sequence globally.
with n=3, k=5 you get:
i=0, {0,0,0,0,0}
i=1, {0,0,0,0,1}
i=2, {0,0,0,0,2}
i=3, {0,0,0,1,1}
i=4, {0,0,0,1,2}
i=5, {0,0,0,2,2}
i=6, {0,0,1,1,1}
i=7, {0,0,1,1,2}
i=8, {0,0,1,2,2}
i=9, {0,0,2,2,2}
i=10, {0,1,1,1,1}
i=11, {0,1,1,1,2}
i=12, {0,1,1,2,2}
i=13, {0,1,2,2,2}
i=14, {0,2,2,2,2}
i=15, {1,1,1,1,1}
i=16, {1,1,1,1,2}
i=17, {1,1,1,2,2}
i=18, {1,1,2,2,2}
i=19, {1,2,2,2,2}
i=20, {2,2,2,2,2}
The algorithm for generating it can be seen as MBnext_multicombination at http://www.martinbroadhurst.com/combinatorial-algorithms.html
Update:
So i thought i'd replace the binomial coefficient in pascals triangle with (n+k-1)!/(k!(n-1)!) to see how it looks.
(* Mathematica code to display pascal and other triangle *)
t1 = Table[Binomial[n, k], {n, 0, 8}, {k, 0, n}];
t2 = Table[(n + k - 1)!/(k! (n - 1)!), {n, 0, 8}, {k, 0, n}];
(*display*)
{Row[#, "\t"]} & /@ t1 // Grid
{Row[#, "\t"]} & /@ t2 // Grid
T1:
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
1 7 21 35 35 21 7 1
1 8 28 56 70 56 28 8 1
T2:
Indeterminate
1 1
1 2 3
1 3 6 10
1 4 10 20 35
1 5 15 35 70 126
1 6 21 56 126 252 462
1 7 28 84 210 462 924 1716
1 8 36 120 330 792 1716 3432 6435
Comparing with the n=3,k=5 console output at the start of this post: the third diagonal {3,6,10,15,21,28,36} gives the index of each roll-over point {0,0,0,1,1} -> {0,0,1,1,1} -> {0,1,1,1,1}, etc. And the diagonal to the left of it seems to show how many values are contained in the previous block (diagonal[2][i] == diagonal[3][i] - diagonal[3][i-1])). And if you read the 5th row of the pyramid horizontally you get the max amount of combinations for increasing values of N in (n+k-1)!/(k!(n-1)!) where K=5.
There is probably a way to use this information to determine the exact combo for an arbitrary index, without enumerating the whole set, but i'm not sure if i need to go that far. The original problem was just to decompose the full combo space into equal subsets, that can be generated locally, and worked on in parallel by the GPU. So the triangle above gives us the starting index of every block, of which the combo can be trivially derived, and all its successive elements incrementally enumerated. It also gives us the block size, and how many total combinations we have. So now it becomes a packing problem of how to fit unevenly sized blocks into groups of equal work load across X amount of threads.
See the example at: https://en.wikipedia.org/wiki/Combinatorial_number_system#Finding_the_k-combination_for_a_given_number
Just replace the binomial Coefficient with (n+k-1)!/(k!(n-1)!).
Assuming n=3,k=5, let's say we want to calculate the 19th combination (id=19).
id=0, {0,0,0,0,0}
id=1, {0,0,0,0,1}
id=2, {0,0,0,0,2}
...
id=16, {1,1,1,1,2}
id=17, {1,1,1,2,2}
id=18, {1,1,2,2,2}
id=19, {1,2,2,2,2}
id=20, {2,2,2,2,2}
The result we're looking for is {1,2,2,2,2}.
Examining our 'T2' triangle: n=3,k=5 points to 21, being the 5th number (top to bottom) of the third diagonal (left to right).
Indeterminate
1 1
1 2 3
1 3 6 10
1 4 10 20 35
1 5 15 35 70 126
1 6 21 56 126 252 462
1 7 28 84 210 462 924 1716
1 8 36 120 330 792 1716 3432 6435
We need to find the largest number in this row (horizontally, not diagonally) that does not exceed our id=19 value. So moving left from 21 we arrive at 6 (this operation is performed by the largest function below). Since 6 is the 2nd number in this row it corresponds to n==2 (or g[2,5] == 6 from the code below).
Now that we've found the 5th number in the combination, we move up a floor in the pyramid, so k-1=4. We also subtract the 6 we encountered below from id, so id=19-6=13. Repeating the entire process we find 5 (n==2 again) to be the largest number less than 13 in this row.
Next: 13-5=8, Largest is 4 in this row (n==2 yet again).
Next: 8-4=4, Largest is 3 in this row (n==2 one more time).
Next: 4-3=1, Largest is 1 in this row (n==1)
So collecting the indices at each stage we get {1,2,2,2,2}
The following Mathematica code does the job:
g[n_, k_] := (n + k - 1)!/(k! (n - 1)!)
largest[i_, nn_, kk_] := With[
{x = g[nn, kk]},
If[x > i, largest[i, nn-1, kk], {nn,x}]
]
id2combo[id_, n_, 0] := {}
id2combo[id_, n_, k_] := Module[
{val, offset},
{val, offset} = largest[id, n, k];
Append[id2combo[id-offset, n, k-1], val]
]
Update:
The order that the combinations were being generated by MBnext_multicombination wasn't matching id2combo, so i don't think they were lexicographic. The function below generates them in the same order as id2combo and matches the order of mathematica's Sort[]function on a list of lists.
void next_combo(unsigned int *ar, unsigned int n, unsigned int k)
{
unsigned int i, lowest_i;
for (i=lowest_i=0; i < k; ++i)
lowest_i = (ar[i] < ar[lowest_i]) ? i : lowest_i;
++ar[lowest_i];
i = (ar[lowest_i] >= n)
? 0 // 0 -> all combinations have been exhausted, reset to first combination.
: lowest_i+1; // _ -> base incremented. digits to the right of it are now zero.
for (; i<k; ++i)
ar[i] = 0;
}