在Node.js中高效计算n个数中选k个数的组合数

Question

我在一个Node.js服务器上有一些性能敏感的代码需要计算组合数。从这个SO答案中，我使用了这个简单的递归函数来计算n个选择k个：

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

因为我们都知道迭代几乎总是比递归更快，所以我根据乘法公式编写了这个函数：

function choosei(n,k){
    var result = 1;
    for(var i=1; i <= k; i++){
        result *= (n+1-i)/i;
    }
    return result;
}

我在自己的机器上运行了几个 基准测试。以下是其中一个测试的结果：

Recursive x 178,836 ops/sec ±7.03% (60 runs sampled)
Iterative x 550,284 ops/sec ±5.10% (51 runs sampled)
Fastest is Iterative

结果始终表明，在Node.js中，迭代方法确实比递归方法快3到4倍左右（至少在我的机器上）。

这对我的需求来说可能足够快了，但有没有办法让它更快呢？我的代码需要频繁调用这个函数，有时候n和k的值相当大，所以速度越快越好。

编辑

在使用le_m和Mike的解决方案进行了几次测试后，结果显示虽然两种方法都比我提出的迭代方法快得多，但是使用Pascal三角形的Mike方法似乎比le_m的对数表方法稍微快一些。

Recursive x 189,036 ops/sec ±8.83% (58 runs sampled)
Iterative x 538,655 ops/sec ±6.08% (51 runs sampled)
LogLUT x 14,048,513 ops/sec ±9.03% (50 runs sampled)
PascalsLUT x 26,538,429 ops/sec ±5.83% (62 runs sampled)
Fastest is PascalsLUT

在我的测试中，使用对数查找方法比迭代方法快了26-28倍，而使用Pascal三角形的方法比对数查找方法快1.3到1.8倍。

请注意，我按照le_m的建议使用mathjs预先计算了具有更高精度的对数，然后将它们转换回常规JavaScript Numbers（始终为双精度64位浮点数）。

- NanoWizard

---

记忆化可以是最快的选择，如果空间不是您关心的问题。 - Godfather

@Godfather @NanoWizard 如果你沒有空間來優化choosei的結果，你可以通過優化階乘來獲得良好的性能提升。你甚至可以使用動態編程技術來計算階乘fact(n) = max_known_fact_value(n) * [i for evey int for max_known_fact_int(n) to n]，這將節省大量時間。 - gbtimmon

@Mike'Pomax'Kamermans 为什么要构建三角形，当你只需要(log)阶乘来获得O(1)的解决方案呢？ - le_m

@Mike'Pomax'Kamermans 我觉得我没有理解错误：帕斯卡三角形LUT（查找表）需要O(n²)的空间，这可能会有问题。 - le_m

虽然没有实际用途，但使用UINT16，您可以在380字节+数组开销中存储LUT，最多可达n=19。 UINT32允许在略大于2kb的情况下最多达到n=35。 UINT64允许最多达到n=68，现在我们已经超越了“您真的认为您需要这么高的二项式系数吗？”的值，略大于18kb。如果已经致力于进行二项式计算，则这些是任何理智人都不会关心的mallocs。特别是在Node.js中，这些LUT与即使是基本Node的内存占用相比微不足道。 - Mike 'Pomax' Kamermans

Answer 1

以下算法在使用线性查找表的对数阶乘时，具有时间复杂度 O(1)，空间复杂度为 O(n)。

由于binomial(1000,500)已经接近Number.MAX_VALUE，因此将n和k限制在范围[0,1000]内是有意义的。我们需要一个大小为1000的查找表。

在现代JavaScript引擎中，n个数字的紧凑数组大小为n*8字节。因此，完整的查找表需要8千字节的内存。如果我们将输入限制在范围[0,100]内，则表格只占用800字节。

var logf = [0, 0, 0.6931471805599453, 1.791759469228055, 3.1780538303479458, 4.787491742782046, 6.579251212010101, 8.525161361065415, 10.60460290274525, 12.801827480081469, 15.104412573075516, 17.502307845873887, 19.987214495661885, 22.552163853123425, 25.19122118273868, 27.89927138384089, 30.671860106080672, 33.50507345013689, 36.39544520803305, 39.339884187199495, 42.335616460753485, 45.38013889847691, 48.47118135183523, 51.60667556776438, 54.78472939811232, 58.00360522298052, 61.261701761002, 64.55753862700634, 67.88974313718154, 71.25703896716801, 74.65823634883016, 78.0922235533153, 81.55795945611504, 85.05446701758152, 88.58082754219768, 92.1361756036871, 95.7196945421432, 99.33061245478743, 102.96819861451381, 106.63176026064346, 110.32063971475739, 114.0342117814617, 117.77188139974507, 121.53308151543864, 125.3172711493569, 129.12393363912722, 132.95257503561632, 136.80272263732635, 140.67392364823425, 144.5657439463449, 148.47776695177302, 152.40959258449735, 156.3608363030788, 160.3311282166309, 164.32011226319517, 168.32744544842765,  172.3527971391628, 176.39584840699735, 180.45629141754378, 184.53382886144948, 188.6281734236716, 192.7390472878449, 196.86618167289, 201.00931639928152, 205.1681994826412, 209.34258675253685, 213.53224149456327, 217.73693411395422, 221.95644181913033, 226.1905483237276, 230.43904356577696, 234.70172344281826, 238.97838956183432, 243.2688490029827, 247.57291409618688, 251.8904022097232, 256.22113555000954, 260.5649409718632, 264.9216497985528, 269.2910976510198, 273.6731242856937, 278.0675734403661, 282.4742926876304, 286.893133295427, 291.3239500942703, 295.76660135076065, 300.22094864701415, 304.6868567656687, 309.1641935801469, 313.65282994987905, 318.1526396202093, 322.66349912672615, 327.1852877037752, 331.7178871969285, 336.26118197919845, 340.815058870799, 345.37940706226686, 349.95411804077025, 354.5390855194408, 359.1342053695754, 363.73937555556347];

function binomial(n, k) {
    return Math.exp(logf[n] - logf[n-k] - logf[k]);
}

console.log(binomial(5, 3));

说明

从原始的迭代算法开始，我们首先将乘积替换为对数之和：

function binomial(n, k) {
    var logresult = 0;
    for (var i = 1; i <= k; i++) {
        logresult += Math.log(n + 1 - i) - Math.log(i);
    }
    return Math.exp(logresult);
}

我们的循环现在对k个术语求和。如果我们重新排列总和，我们可以轻松地看到我们对连续的对数 log(1) + log(2) + ... + log(k) 等进行求和，我们可以用一个 sum_of_logs(k) 来代替它，实际上它与 log(k!) 相同。预先计算这些值并将它们存储在我们的查找表 logf 中，然后就可以得到上面提到的单行算法。

计算查找表:

我建议使用更高的精度预先计算查找表，并将结果元素转换为64位浮点数。如果您不需要额外的那一点精度或想在客户端运行此代码，请使用以下代码：

var size = 1000, logf = new Array(size);
logf[0] = 0;
for (var i = 1; i <= size; ++i) logf[i] = logf[i-1] + Math.log(i);

数值精度：

通过使用对数阶乘，我们避免了存储原始阶乘时固有的精度问题。

我们甚至可以使用斯特林公式来计算 log(n!)，而不需要查找表，在运行时间和空间复杂度O(1)的情况下仍可得到12位有效数字：

function logf(n) {
    return n === 0 ? 0 : (n + .5) * Math.log(n) - n + 0.9189385332046728 + 0.08333333333333333 / n - 0.002777777777777778 * Math.pow(n, -3);
}

function binomial(n , k) {
    return Math.exp(logf(n) - logf(n - k) - logf(k));
}

console.log(binomial(1000, 500)); // 2.7028824094539536e+299

- le_m

---

算术运算是在64位浮点值上执行的，因此您需要在打印之前通过截断无关小数位来格式化结果。 - le_m

我认为这是一个打字错误 - log(1) + log(2) + ... + log(k) 应该是 log(k!) 而不是 log(n!)。 - Boatmarker

Answer 2

永远不要计算阶乘，因为它们增长得太快了。相反，计算你想要的结果。在这种情况下，你想要的是二项式数，它们有一个非常简单的几何构造：你可以按照需要构建pascal's triangle，并使用纯算术来完成。

从[1]和[1,1]开始。下一行在开头有[1]，中间有[1+1]，结尾有[1]：[1,2,1]。下一行：在开头有[1]，在第2个位置上是前两个数字之和，在第三个位置上是接下来两个数字之和，结尾有[1]：[1,3,3,1]。下一行：[1]，然后是1+3=4，在第二个位置上是3+3=6，在第三个位置上是3+1=4，最后是[1]，依此类推。正如你所看到的，没有阶乘、对数或甚至乘法：只有干净整数数字的超快速加法。如此简单，你可以手工构建一个大型查找表。

而且你应该这样做。

不要在代码中计算可以手动计算并作为常量进行即时查找的内容。在这种情况下，编写n=20左右的表格绝对是微不足道的，然后您可以将其用作“起始LUT”，可能甚至从未访问高行。

但是，如果您需要它们或更多，则因为无法构建无限的查找表而进行妥协：您从预先指定的LUT开始，并使用可以将其“填满”到尚未包含所需项的函数：

// step 1: a basic LUT with a few steps of Pascal's triangle
const binomials = [
  [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],
  //  ...
];

// step 2: a function that builds out the LUT if it needs to.
module.exports = function binomial(n,k) {
  while(n >= binomials.length) {
    let s = binomials.length;
    let nextRow = [];
    nextRow[0] = 1;
    for(let i=1, prev=s-1; i<s; i++) {
      nextRow[i] = binomials[prev][i-1] + binomials[prev][i];
    }
    nextRow[s] = 1;
    binomials.push(nextRow);
  }
  return binomials[n][k];
};

由于这是一个整数数组，所以它的内存占用非常小。对于大量涉及二项式的工作，我们实际上甚至不需要每个整数超过两个字节，使这成为一个微不足道的查找表：在n = 19之前，我们不需要超过2个字节，而完整的查找表最多到n = 19只占用微不足道的380字节。与您的程序的其余部分相比，这算不了什么。即使我们允许32位整数，我们也可以在仅2380字节的情况下达到n = 35。

因此，查找速度很快：如果是先前计算过的值，则为O（常数），如果我们没有LUT，则为（n *（n + 1））/ 2步（在大O符号中，这将是O（n²），但大O符号几乎从来不是正确的复杂度度量），对于我们尚未在LUT中找到的术语，则介于两者之间。为您的应用程序运行一些基准测试，这将告诉您初始LUT应该有多大，只需硬编码即可（严肃点。这些是常量，它们是应该硬编码的确切值）。保留生成器以防万一。

然而，请记住您处于JavaScript领域，并受到JavaScript数值类型的限制：整数仅限于2 ^ 53, 超出此范围，整数属性（每个n都有一个不同的m = n + 1 ，使得m-n = 1 ）不能保证。虽然这几乎永远不会成为问题：一旦我们达到了这个限制，我们就正在处理您甚至不应该使用的二项式系数。

- Mike 'Pomax' Kamermans

---

lut数组的最后一行存在错误。对于任何复制粘贴的人，请将54替换为56或删除该行。懒人证明：https://www.google.com/search?q=8+choose+3 - Basic Block

我没想到你会看到这个评论，但很好你修复了它。 - Basic Block

对于那些喜欢一点JS疯狂的人：这是while循环中let s = ...之后的部分，可以作为单行代码：binomials.push(new Array(s+1).fill(0).map((_, i) => i == 0 || i == s ? 1 : binomials[s-1][i-1] + binomials[s-1][i])); - sigalor

@sigalor 不需要使用 new Array(s+1).fill(0).map()，只需使用现代 JS：[...Array(s+1)].map()。如果你试图进行微观优化，而这种优化实际上并不会使执行速度更快 =) - Mike 'Pomax' Kamermans

作为几年后的评论，我们现在有了Uint8Array和Uint16Array，它们是大小固定的数组，因此您现在可以使用new Uint8Array(rowsize)或（完全过度）new Uint16Array(rowsize)来_强制_内存占用。 - Mike 'Pomax' Kamermans

Answer 3

使用帕斯卡三角形是计算n个中选k个的快速方法。

我知道的最快的方法是利用“计算阶乘的复杂度”中的结果。只需计算所有3个阶乘，然后执行两个除法运算，每个运算的复杂度为M(n logn)。

- Quantum Guy 123

---

在“只是计算阶乘”中没有“只是”：计算n的阶乘需要n-1次乘法运算，相对于加法来说，这是CPU非常擅长的领域。你引用的论文只讨论了关于n的复杂度，而没有涉及到所谓的CPU操作，这显然与问题本身要求的高效计算有所偏差。问题是关于高效计算的。 - Mike 'Pomax' Kamermans