仅基于索引计算第N个有重复元素的多元组组合

Question

仅基于索引计算第N个有重复元素的多元组组合

combinationscombinatoricsfactorialbinomial-coefficients

4

如何根据索引计算第N个组合。有(n+k-1)!/(k!(n-1)!)种带重复的组合。

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}

因此，black_magic_function(3) 应该产生 {0,0,1,1,1}。

这将要传入 GPU 着色器，所以我希望每个工作组/线程能够在不必全局存储序列的情况下，找出它们的排列子集。

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}

生成它的算法可以在http://www.martinbroadhurst.com/combinatorial-algorithms.html上看到，算法名称为MBnext_multicombination。

更新：

所以我想用(n+k-1)!/(k!(n-1)!)替换帕斯卡三角形中的二项式系数，看看效果如何。

(* 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

与本文开头的n=3，k=5的控制台输出相比较：第三条对角线{3,6,10,15,21,28,36}给出了每个翻转点{0,0,0,1,1} -> {0,0,1,1,1} -> {0,1,1,1,1}的索引，等等。左边的对角线似乎显示了前一个块中包含多少值（diagonal[2][i]==diagonal[3][i]-diagonal[3][i-1]）。如果你横向阅读金字塔的第五行，你会得到递增N的最大组合数，公式为 (n+k-1)!/(k!(n-1)!)，其中K=5。

可能有一种方法可以利用这些信息来确定任意索引的确切组合，而无需枚举整个集合，但我不确定是否需要走这么远。最初的问题只是将完整的组合空间分解成相等的子集，这些子集可以在GPU上本地生成，并通过并行处理。因此，上面的三角形给出了每个块的起始索引，其组合可以轻松地派生出来，所有连续的元素都可以逐步枚举。它还告诉我们块的大小以及总组合数。现在，问题变成将不均匀大小的块放入平均工作负载的X个线程的组中的装箱问题。

- Jean-Baptiste

n和k的最大值是多少？ - Bob Bryan

K 的最大值为 6。N 的最大值不能超过 300，但通过一些预处理技巧可以将其降至约为 20、50 或 100。 - Jean-Baptiste

因此，当n = 2时，情况相当简单 - 只需从右侧开始放置由k指定的许多1，并将其余部分保留为零。对于n = 3，k = 5的情况，我得到(3+5-1)! / 5!(3-1)! = 21。但是，我不太确定组合应该是什么样子的。您能否编辑您的问题并展示n = 3和k = 5的组合？ - Bob Bryan

人们会从编辑中收到线程更新通知吗？ - Jean-Baptiste

是的，他们会。通常您需要指定个人的姓名。例如，为了确保我得到通知，您会指定@Bob Bryan。但是，除了最后一条评论外，我已经收到了您的评论。 - Bob Bryan

3个回答

0

这里是一个组合数系统的实现，它可以处理带有和不带有重复（即多重集合）的组合，并且可选择生成字典序排序。

/// Combinatorial number system encoder/decoder
/// https://en.wikipedia.org/wiki/Combinatorial_number_system
struct CNS(
    /// Type for packed representation
    P,
    /// Type for one position in unpacked representation
    U,
    /// Number of positions in unpacked representation
    size_t N,
    /// Cardinality (maximum value plus one) of one position in unpacked representation
    U unpackedCard,
    /// Produce lexicographic ordering?
    bool lexicographic,
    /// Are repetitions representable? (multiset support)
    bool multiset,
)
{
static:
    /// Cardinality (maximum value plus one) of the packed representation
    static if (multiset)
        enum P packedCard = multisetCoefficient(unpackedCard, N);
    else
        enum P packedCard = binomialCoefficient(unpackedCard, N);
    alias Index = P;

    private P summand(U value, Index i)
    {
        static if (lexicographic)
        {
            value = cast(U)(unpackedCard-1 - value);
            i = cast(Index)(N-1 - i);
        }
        static if (multiset)
            value += i;
        return binomialCoefficient(value, i + 1);
    }

    P pack(U[N] values)
    {
        P packed = 0;
        foreach (Index i, value; values)
        {
            static if (!multiset)
                assert(i == 0 || value > values[i-1]);
            else
                assert(i == 0 || value >= values[i-1]);

            packed += summand(value, i);
        }

        static if (lexicographic)
            packed = packedCard-1 - packed;
        return packed;
    }

    U[N] unpack(P packed)
    {
        static if (lexicographic)
            packed = packedCard-1 - packed;

        void unpackOne(Index i, ref U r)
        {
            bool checkValue(U value, U nextValue)
            {
                if (summand(nextValue, i) > packed)
                {
                    r = value;
                    packed -= summand(value, i);
                    return true;
                }
                return false;
            }

            // TODO optimize: (rolling product / binary search / precomputed tables)
            // TODO optimize: don't check below N-i
            static if (lexicographic)
            {
                foreach_reverse (U value; 0 .. unpackedCard)
                    if (checkValue(value, cast(U)(value - 1)))
                        break;
            }
            else
            {
                foreach         (U value; 0 .. unpackedCard)
                    if (checkValue(value, cast(U)(value + 1)))
                        break;
            }
        }

        U[N] values;
        static if (lexicographic)
            foreach         (Index i, ref r; values)
                unpackOne(i, r);
        else
            foreach_reverse (Index i, ref r; values)
                unpackOne(i, r);

        return values;
    }
}

完整代码：https://gist.github.com/CyberShadow/67da819b78c5fd16d266a1a3b4154203

- Vladimir Panteleev

-1

我对这个问题做了一些初步的分析。在谈论我找到的低效解决方案之前，让我先给您提供一篇论文的链接，该论文介绍如何将k-索引（或组合）转换为与二项式系数相关联的组合的排名或字典序索引：

http://tablizingthebinomialcoeff.wordpress.com/

我在解决这个问题的时候也是以同样的方式开始的。我想出了以下代码，它在公式(n+k-1)!/k!(n-1)!中为每个k值使用一个循环，当k = 5时。按照现有的写法，这段代码将生成n选5的所有组合：

private static void GetCombos(int nElements)
{
   // This code shows how to generate all the k-indexes or combinations for any number of elements when k = 5.
   int k1, k2, k3, k4, k5;
   int n = nElements;
   int i = 0;
   for (k5 = 0; k5 < n; k5++)
   {
      for (k4 = k5; k4 < n; k4++)
      {
         for (k3 = k4; k3 < n; k3++)
         {
            for (k2 = k3; k2 < n; k2++)
            {
               for (k1 = k2; k1 < n; k1++)
               {
                  Console.WriteLine("i = " + i.ToString() + ", " + k5.ToString() + " " + k4.ToString() + 
                     " " + k3.ToString() + " " + k2.ToString() + " " + k1.ToString() + " ");
                  i++;
               }
            }
         }
      }
   }
}

该方法的输出为：

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

这是与您在编辑答案中提供的相同值。我也尝试了4选5，看起来它也生成了正确的组合。

我用C#编写了这个代码，但是您应该能够在其他语言（如C/C++、Java或Python）中使用它而不需要太多修改。

一个有点低效的解决方案的想法是修改GetCombos以接受k作为输入。由于k限制为6，因此可以对k进行测试。因此，生成n选k情况下所有可能组合的代码将如下所示：

private static void GetCombos(int k, int nElements)
{
   // This code shows how to generate all the k-indexes or combinations for any n choose k, where k <= 6.
   //
   int k1, k2, k3, k4, k5, k6;
   int n = nElements;
   int i = 0;
   if (k == 6)
   {
      for (k6 = 0; k6 < n; k6++)
      {
         for (k5 = 0; k5 < n; k5++)
         {
            for (k4 = k5; k4 < n; k4++)
            {
               for (k3 = k4; k3 < n; k3++)
               {
                  for (k2 = k3; k2 < n; k2++)
                  {
                     for (k1 = k2; k1 < n; k1++)
                     {
                        Console.WriteLine("i = " + i.ToString() + ", " + k5.ToString() + " " + k4.ToString() +
                           " " + k3.ToString() + " " + k2.ToString() + " " + k1.ToString() + " ");
                        i++;
                     }
                  }
               }
            }
         }
      }
   }
   else if (k == 5)
   {
      for (k5 = 0; k5 < n; k5++)
      {
         for (k4 = k5; k4 < n; k4++)
         {
            for (k3 = k4; k3 < n; k3++)
            {
               for (k2 = k3; k2 < n; k2++)
               {
                  for (k1 = k2; k1 < n; k1++)
                  {
                     Console.WriteLine("i = " + i.ToString() + ", " + k5.ToString() + " " + k4.ToString() +
                        " " + k3.ToString() + " " + k2.ToString() + " " + k1.ToString() + " ");
                     i++;
                  }
               }
            }
         }
      }
   }
   else if (k == 4)
   {
      // One less loop than k = 5.
   }
   else if (k == 3)
   {
      // One less loop than k = 4.
   }
   else if (k == 2)
   {
      // One less loop than k = 3.
   }
   else
   {
      // k = 1 - error?
   }
}

现在，我们有一个可以生成所有感兴趣组合的方法。但是，问题是如何从字典顺序或排名中获取特定的组合。这可以通过简单计数并在达到指定值时返回正确的组合来实现。因此，为了适应这一点，需要向该方法添加表示等级的额外参数。因此，执行此操作的新函数如下：

private static int[] GetComboOfRank(int k, int nElements, int Rank)
{
   // Gets the combination for the rank using the formula (n+k-1)!/k!(n-1)! where k <= 6.
   int k1, k2, k3, k4, k5, k6;
   int n = nElements;
   int i = 0;
   int[] ReturnArray = new int[k];
   if (k == 6)
   {
      for (k6 = 0; k6 < n; k6++)
      {
         for (k5 = 0; k5 < n; k5++)
         {
            for (k4 = k5; k4 < n; k4++)
            {
               for (k3 = k4; k3 < n; k3++)
               {
                  for (k2 = k3; k2 < n; k2++)
                  {
                     for (k1 = k2; k1 < n; k1++)
                     {
                        if (i == Rank)
                        {
                           ReturnArray[0] = k1;
                           ReturnArray[1] = k2;
                           ReturnArray[2] = k3;
                           ReturnArray[3] = k4;
                           ReturnArray[4] = k5;
                           ReturnArray[5] = k6;
                           return ReturnArray;
                        }
                        i++;
                     }
                  }
               }
            }
         }
      }
   }
   else if (k == 5)
   {
      for (k5 = 0; k5 < n; k5++)
      {
         for (k4 = k5; k4 < n; k4++)
         {
            for (k3 = k4; k3 < n; k3++)
            {
               for (k2 = k3; k2 < n; k2++)
               {
                  for (k1 = k2; k1 < n; k1++)
                  {
                     if (i == Rank)
                     {
                        ReturnArray[0] = k1;
                        ReturnArray[1] = k2;
                        ReturnArray[2] = k3;
                        ReturnArray[3] = k4;
                        ReturnArray[4] = k5;
                        return ReturnArray;
                     }
                     i++;
                  }
               }
            }
         }
      }
   }
   else if (k == 4)
   {
      // Same code as in the other cases, but with one less loop than k = 5.
   }
   else if (k == 3)
   {
      // Same code as in the other cases, but with one less loop than k = 4.
   }
   else if (k == 2)
   {
      // Same code as in the other cases, but with one less loop than k = 3.
   }
   else
   {
      // k = 1 - error?
   }
   // Should not ever get here.  If we do - it is some sort of error.
   throw ("GetComboOfRank - did not find rank");
}

ReturnArray返回与排名相关联的组合。所以，这段代码应该适用于您。但是，如果进行表查找，可以实现比这更快的速度。300选6的问题在于：

300 choose 6 = 305! / (6!(299!) = 305*304*303*302*301*300 / 6! = 1,064,089,721,800

这可能是存储在内存中的数据量过大了。因此，如果您可以通过预处理将n降至20，则总计会是：

20 choose 6 = 25! / (6!(19!)) = 25*24*23*22*21*20 / 6! = 177,100
20 choose 5 = 24! / (5!(19!)) = 24*23*22*21,20 / 5!    =  42,504
20 choose 4 = 23! / (4!(19!)) = 23*22*21*20 / 4!       =   8,855
20 choose 3 = 22! / (3!(19!)) = 22*21*20 / 3!          =   1,540
20 choose 2 = 21! / (2!(19!)) = 22*21 / 2!             =     231
                                                         =======
                                                         230,230

如果每个组合的值都使用一个字节，则可以计算出在内存中存储表格（通过锯齿形数组或者可能是5个单独的表格）所使用的总字节数为：

177,100 * 6 = 1,062,600
 42,504 * 5 =   212,520
  8,855 * 4 =    35,420
  1,540 * 3 =     4,620
    231 * 2 =       462
              =========
              1,315,622

这取决于目标计算机和可用内存量，但是在如今许多计算机具有千兆字节可用内存的情况下，1,315,622字节并不算太多。

- Bob Bryan

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jean-Baptiste · Accepted Answer

请参考以下示例： https://en.wikipedia.org/wiki/Combinatorial_number_system#Finding_the_k-combination_for_a_given_number 只需将二项式系数替换为(n+k-1)!/(k!(n-1)!)。

假设n=3,k=5，我们想计算第19个组合（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}

我们要得到的结果是 {1,2,2,2,2}。

解析 'T2' 三角形： n=3,k=5 指向的是第三条斜线（从左到右），从上往下数第五个数，即 21。

            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

我们需要在这一行中找到不超过我们的 id=19 值的最大数（横向，而不是对角线）。因此，从 21 左移，我们到达 6（该操作由下面的 largest 函数执行）。由于 6 是此行中的第二个数字，因此它对应于 n==2（或者来自下面代码的 g[2,5] == 6）。

现在，我们已经找到了组合中的第5个数字，我们向金字塔中的上一层移动，所以 k-1=4。我们还将下面遇到的 6 从 id 中减去，因此 id=19-6=13。重复整个过程，我们发现 5（再次是n==2）是小于 13 的最大数字。

接下来：13-5=8，此行中最大的是 4 （再次是n==2）。

接下来：8-4=4，此行中最大的是 3 （再次是n==2）。

接下来：4-3=1，此行中最大的是 1（n==1）。

因此，在每个阶段收集索引，我们得到{1,2,2,2,2} 以下是 Mathematica 代码：

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]
]

更新： MBnext_multicombination 生成组合的顺序与 id2combo 不匹配，因此我认为它们不是字典序。下面的函数按照与 id2combo 相同的顺序生成它们，并且与 Mathematica 的 Sort[] 函数在列表的列表上的顺序匹配。

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