计算n选k的最佳方法是什么?

44

什么是评估"nk"值的最有效方法?我认为暴力方法是通过单独计算每个阶乘来找到n! / k! / (n-k)!.

更好的策略可能是根据这个递归公式使用DP,nCk == (n-1)C(k-1) + (n-1)C(k)。有没有其他更好的方法来评估nk,以复杂度和避免溢出的风险为考虑?


阶乘计算在空间和时间方面都比您的递归替代方案更有效率。 - SomeWittyUsername
2
首先,你可以用 n*(n-1)*(n-2)*...*(k+1) 代替 n!/k!。 许多因子都会被消除,完全计算 n!k! 没有意义。 - Tim Goodman
你考虑的n的范围是什么? - Andrew Morton
@AndrewMorton 我需要计算 n 选 k,其中 n < 1000000,k < 1000000。 - Nikunj Banka
6
@M42:这个问题不是你链接到的那个问题的副本。那个问题要求从n个元素中选择k个元素的所有组合,而这个问题只想知道这样的组合数量。 - Luke Woodward
13个回答

56

以下是我的版本,它完全使用整数(k的除法始终产生整数商),速度为O(k):

function choose(n, k)
    if k == 0 return 1
    return (n * choose(n - 1, k - 1)) / k

我递归写它是因为它非常简单而且漂亮,但如果你喜欢,可以把它转换成迭代解决方案。


1
这不是O(k); k 严格小于 n,因此您不能忽略 n 对运行时间的贡献。最好的情况下,您可以说它是O(k M(n)),其中M(n)是您的乘法算法的速度。 - chepner
6
正确,但技术性强。上述函数会进行O(k)次乘除运算。我忽略了操作本身的位复杂度。 - user448810
6
@icepack:不,它并不是这样的。分子的范围从n到n-k+1,分母的范围从k到1。因此,choose(9,4) = (9876) / (4321) = 126,这是正确的。相比之下,9!/4! = 362880 / 24 = 15120。 - user448810
1
这是递归形式的乘法方法。它确实是O(k),除非使用Stirling近似估计k!足够好(http://en.wikipedia.org/wiki/Stirling%27s_approximation),否则这已经是最快的方法了。还有一种分治版本的阶乘算法*可能*也会有所帮助(http://gmplib.org/manual/Factorial-Algorithm.html)。 - Pedrom
我认为在C++中这很容易溢出,会得到各种奇怪的结果和/或错误... 我尝试计算一些范围内的n选择n-k的总和:http://pastebin.com/jhNAp8KD - shinzou
显示剩余2条评论

41

19
当 n-k < k 时,我们可以通过计算 (n 选 n-k) 而不是 (n 选 k) 来加快速度。 - Nikunj Banka
3
请注意,(n - (k-i)) / i可能不是整数。 - pomber
1
@pomber 个别的 factor (n - (k-i)) / i 可能实际上不是整数,但是从 x=1y 的因子积将总是可被 y 整除的(因为恰好有 y 个连续的整数)。 - Lambder
2
稍微更新的公式:https://wikimedia.org/api/rest_v1/media/math/render/svg/652661edd20c8121e58a2b26844ce46c24180f0f - trig-ger
这些阶乘的计算成本相当高。使用此处所见方法可以加速:https://cs.stackexchange.com/a/14476 。 - Quantum Guy 123

9
使用Pascal三角形计算二项式系数(n choose k)是避免溢出最简单的方法。不需要分数或乘法。(n choose k)。Pascal三角形的第n行和第k个条目给出该值。 看看这个页面。这是一个O(n^2)操作,只包含加法,可以用动态规划解决。对于能够适合64位整数的任何数字,它都会非常快速。

4
和 O(n^2) 的额外空间 - SomeWittyUsername
@rici 正确,但那只是对所提出方法的优化。 - SomeWittyUsername

6
如果你要计算很多这样的组合,那么计算帕斯卡三角形肯定是最佳选择。因为你已经知道了递归公式,所以我想在这里贴一些代码:
MAX_N = 100
MAX_K = 100

C = [[1] + [0]*MAX_K for i in range(MAX_N+1)]

for i in range(1, MAX_N+1):
    for j in range(1, MAX_K+1):
        C[i][j] = C[i-1][j-1] + C[i-1][j];

print C[10][2]
print C[10][8]
print C[10][3]

3

在遇到类似问题后,我决定编译最佳解决方案并对每个方案进行了一些不同值的n和k的简单测试。起初我有大约10个函数,逐渐淘汰出那些明显错误或在特定值处停止工作的函数。在所有的解决方案中,用户user448810提供的答案是最简洁、最易于实现的,我非常喜欢它。下面是包括我运行的每个测试的次数、每个函数的代码、函数的输出以及得到该输出所花费的时间的代码。我只运行了20000次,在重新运行测试时仍存在时间波动,但您应该可以获得每个函数的整体表现。

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    //EXPECTED VALUES
    //x choose x = 1
    //9 choose 4 =126
    //52 choose 5 = 2598960;
    //64 choose 33 = 1.777090076065542336E18;

    //# of runs for each test: 20000
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//https://dev59.com/3WnWa4cB1Zd3GeqP4dXd#12983878
public static double combination(long n, long k) {
    double sum=0;
    for(long i=0;i<k;i++) {
        sum+=Math.Log10(n-i);
        sum-=Math.Log10(i+1);
    }
    return Math.Pow(10, sum);
}
/*
    10 choose 10
    0.9999999999999992
    Elapsed=00:00:00.0000015
    9 choose 4
    126.00000000000001
    Elapsed=00:00:00.0000009
    52 choose 5
    2598959.9999999944
    Elapsed=00:00:00.0000013
    64 choose 33
    1.7770900760655124E+18
    Elapsed=00:00:00.0000058
*/


//........................................................
//https://dev59.com/3WnWa4cB1Zd3GeqP4dXd#19125294
public static double BinomCoefficient(long n, long k) {
    if (k > n) return 0;
    if (n == k) return 1; // only one way to chose when n == k
    if (k > n-k) k = n-k; // Everything is symmetric around n-k, so it is quicker to iterate over a smaller k than a larger one.
    double c = 1;
    for (long i = 1; i <= k; i++) {
        c *= n--;
        c /= i;
    }
    return c;
}
/*
    10 choose 10
    1
    Elapsed=00:00:00
    9 choose 4
    126
    Elapsed=00:00:00.0000001
    52 choose 5
    2598960
    Elapsed=00:00:00.0000001
    64 choose 33
    1.7770900760655432E+18
    Elapsed=00:00:00.0000006
*/

//........................................................
//https://dev59.com/UWUp5IYBdhLWcg3wY29V#15302448
public static double choose(long n, long k) {
    if (k == 0) return 1;
    return (n * choose(n-1, k-1)) / k;
}
/*
    10 choose 10
    1
    Elapsed=00:00:00.0000002
    9 choose 4
    126
    Elapsed=00:00:00.0000003
    52 choose 5
    2598960
    Elapsed=00:00:00.0000004
    64 choose 33
    1.777090076065543E+18
    Elapsed=00:00:00.0000008
*/

//........................................................
//My own version which is just a mix of the two best above.
public static double binomialCoeff(int n, int k) {
    if (k > n) return 0;
    if (k > n-k) k = n-k; // Everything is symmetric around n-k, so it is quicker to iterate over a smaller k than a larger one.
    double recusion(long n, long k) {
        if (k == 0) return 1; // only one way to chose when n == k
        return (n * recusion(n-1, k-1)) / k;
    }
    return recusion(n,k);
}
/*
    10 choose 10
    1
    Elapsed=00:00:00
    9 choose 4
    126
    Elapsed=00:00:00.0000001
    52 choose 5
    2598960
    Elapsed=00:00:00.0000002
    64 choose 33
    1.777090076065543E+18
    Elapsed=00:00:00.0000007
*/

//........................................................
//https://en.wikipedia.org/wiki/Binomial_coefficient
public static double binomial(long n, long k) {
    // take advantage of symmetry
    if (k > n-k)  k = n-k;
    
    long c = 1;
    for (long i = 1; i <= k; i++, n--) {
        // return 0 on potential overflow
        if (c/i >= long.MaxValue/n) return 0;
        // split c * n / i    into (c / i * i) + (c % i * n / i)
        c = (c / i * n) + (c % i * n / i); 
    }
    
    return c;
}

/*
    10 choose 10
    1
    Elapsed=00:00:00.0000006
    9 choose 4
    126
    Elapsed=00:00:00.0000002
    52 choose 5
    2598960
    Elapsed=00:00:00.0000003
    64 choose 33
    1.7770900760655424E+18
    Elapsed=00:00:00.0000029
*/

这是一种非常优秀的方法,支持更大的数字,但会以略微不精确的结果为代价,这些结果可能无法通过整数转换来修复。 - Ṃųỻịgǻňạcểơửṩ

2
如果您有一个阶乘查找表,那么计算C(n,k)将非常快速。

对于大量的n和k,查找表可能会受到限制。此外,应该提供一个选项来处理超出该表范围的值。 - Pedrom
@Pedrom 在问题中没有提到数字大小的限制。它被标记为“语言无关”和“算法”。 - Andrew Morton

2
n!/k!(n-k)!方法的问题不在于成本,而是在于!增长得非常迅速,以至于即使对于nCk的值来说,这些中间计算也不在64位整数的范围内。如果您不喜欢kainaw的递归加法方法,可以尝试乘法方法: nCk == product(i=1..k) (n-(k-i))/i 其中product(i=1..k)表示当i取值1,2,...,k时所有项的乘积。

你说得没错,溢出的可能性会在最终答案无法适应机器字之前破坏一切,但我不喜欢你的解决方案:它将为某些因子产生分数值,例如当n=4和k=3时i=2的因子。(当然,这些因子最终将相乘得到一个整数,但是你的方法意味着中间结果需要以浮点数形式存储 - yuck!) - j_random_hacker

2

最快的方法可能是使用公式,而不是帕斯卡三角形。当我们知道后面要除以相同的数字时,让我们先不进行乘法运算。 如果 k < n/2,则让 k = n - k。我们知道 C(n,k) = C(n,n-k) 现在:

n! / (k! x (n-k)!) = (product of numbers between (k+1) and n) / (n-k)!

至少采用这种技术,您永远不会将一个数除以您之前乘过的数字。您有(n-k)次乘法和(n-k)次除法。

我正在考虑一种避免所有除法的方法,即通过找到我们必须相乘的数字和我们必须相除的数字之间的最大公约数来实现。稍后我会尝试编辑。


找到最大公约数肯定会减少操作量。不幸的是,找到最大公约数本身将是一个更繁重的任务。 - SomeWittyUsername
是的,我很担心这个。但是当乘法有一个大数时,GCD将在小数上计算。实际上,我不确定GCD是否比除法更难。 - bruce_ricard
我倾向于持怀疑态度,但看到结果会很有趣 :) - SomeWittyUsername
除法和最大公约数对于大小为n的2个数字都具有O(n^2)的复杂度。在这里,我们将计算一个大数和一个小数的除法,而GCD将在两个小数上进行计算,但我们需要对所有数字进行计算。如果我必须手动完成它,我认为我会尝试找到至少明显的倍数和GCD,以避免进行无用的除法。 - bruce_ricard
@double_squeeze: 如果你有C(n,k)的质因数分解,你可以仅通过乘法计算出C(n,k),这难道不是你的目标吗? - rici
显示剩余2条评论

2

如果有帮助到其他人的话,我希望这个答案能够解决问题。


我之前没有看到过这个答案。

def choose(n, k):
   nominator = n
   for i in range(1,k):
      nominator *= (n-i)
      k *= i 
   return nominator/k

1

计算“部分阶乘”,即给定数字范围内的乘积的复杂度如何? - Ṃųỻịgǻňạcểơửṩ

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接