65.9K
CodeProject 正在变化。 阅读更多。
Home

二项式系数的部分和

Tasos Kipriotis
starIconstarIconstarIconstarIconstarIcon

5.00/5 (1投票)

2023 年 10 月 23 日

CPOL

9分钟阅读

viewsIcon

4807

downloadIcon

79

一种更快地计算二项式系数部分和的方法

二项式系数

二项式系数[1]是二项式展开式中系数 \( a_k \) 处的正整数。

$ (x+y)^n = a_0 \ x^n + a_1 \ x^{n-1} \ y \ + \ldots \ + \ a_{n-1} \ x \ y^{n-1} \ + \ a_n \ y^n= \sum_{k=0}^n { a_k \ x^{n-k} \ y^k} $

系数 \( a_k \)，用 \( \dbinom{n}{k} \) 表示，读作“\( n \) choose \( k \)”，表示从 \( n \) 个不同元素中选取 \( k \) 个元素的组合数，其计算公式为：

$ { n \choose k } \ = \ \frac {n!} {k!\ (n-k)!} \tag{1} $

其中包含的阶乘[4]意味着需要进行大量的乘法运算。

通过将分子和分母同时除以 \( (n-k)! \)，公式 (1) 可以写成：

$ { n \choose k } \ = \ \frac {(n-k+1) \cdot \ldots \cdot (n-1) \cdot n} {k!} \ = \ \frac { \prod \limits _{j=n-k+1} ^{n} j } {k!} \tag{2} $

此公式计算速度更快，因为它比公式 (1) 使用的乘法次数更少。

二项式系数之和

二项式系数之和由下式给出：

$\sum _{k=0} ^n = 2^n \tag{3} $

但是，对于这些系数的部分和，即对于下式，并没有一个通用的公式：

$ \sum _{j=a} ^b { n \choose j } = { n \choose a } + { n \choose a + 1 } \ + \ldots \ + \ { n \choose b - 1 } \ + \ { n \choose b } \tag{4} $

其中 \( 0 \leqslant a \leqslant b \leqslant n \)。这意味着必须计算从 \( \dbinom{n}{a} \) 到 \( \dbinom{n}{b} \)（包括两端）的所有系数，同时累加到一个临时结果中。这显然是一个缓慢的过程。

一举两得

使用杨辉三角法则：

$ { n \choose k } + { n \choose { k + 1 } } = { { n + 1 } \choose { k + 1 } } $

公式 (4) 的求和可以改写为：

$ \begin{align} \sum _{j=a} ^b { n \choose j } &= { n \choose a } + { n \choose a + 1 } + { n \choose a + 2 } + { n \choose a + 3 } + \ldots \\ &= { n + 1 \choose a + 1 } + { n + 1 \choose a + 3 } + \ldots \end{align} \tag{5} $

这样可以使计算次数减半。如果系数的数量是奇数，则会剩下一个单独的系数需要计算。

$ \begin{align} \sum _{j=a} ^b { n \choose j } &= { n \choose a } + { n \choose a + 1 } + { n \choose a + 2 } + { n \choose a + 3 } + { n \choose a + 4 } + \ldots \\ &= { n \choose a } + { n + 1 \choose a + 2 } + { n + 1 \choose a + 4 } + \ldots \end{align} \tag{6} $

令 \( s=(b-a+1) \ mod\ 2 \) 和 \( t=\lfloor(b-a+1)/2\rfloor \)，公式 (5) 和 (6) 可以写成：

$ \sum _{j=a} ^b { n \choose j } = s \cdot { n \choose a } \ + \ \sum _{r=0} ^t { n+1 \choose a + 1 + 2 \cdot r + s } \tag{7} $

前两个系数和的计算

再次使用杨辉三角法则：

$ { n \choose k } + { n \choose {k+1} } = { {n+1} \choose {k+1} } = \frac {(n+1)!} {(k+1)!\ \cdot\ (n-k)!} = \frac {\prod \limits _{j=n-k+1} ^{n+1} j} {(k+1)!} $

计算接下来的两个系数之和

一旦我们计算出单个系数，再进行一些数学推导，可以得到：

$ \begin{align} {n \choose {i+1}} + {n \choose {i+2}} &= {{n+1} \choose {i+2}}\\ \\ &=\frac {(n+1)!} {(i+2)!\cdot (n-i-1)!}\\ \\ &= \frac {n! \cdot (n+1)} {i!\cdot (i+1) \cdot (i+2)\cdot \large {\frac {(n-i)!} {n-i}}} \\ \\ &= {n \choose i} \cdot \frac {(n-i) \cdot (n+1)} {(i+1) \cdot (i+2)} \end{align} \tag{8} $

同样，在计算两个系数之和后，经过一些数学推导，我们得到：

$ \begin{align} {n \choose {i+2}} + {n \choose {i+3}} &= {{n+1} \choose {i+3}}\\ \\ &= \frac {(n+1)!} {(i+3)! \cdot (n-i-2)!}\\ \\ &= \frac {(n+1)!} {(i+1)! \cdot (i+2) \cdot (i+3) \cdot \large {\frac {(n-i)!} {(n-i-1) \cdot (n-i)}}}\\ \\ &= {{n+1} \choose {i+1}} \cdot \frac {(n-i-1) \cdot (n-i)} {(i+2) \cdot (i+3)} \end{align} \tag{9} $

无论哪种情况，接下来的两个系数之和都可以快速计算。

一些简单的情况

当 \( k>n \) 时，\( \dbinom{n}{k} =0 \)，因此：

$ \sum _{i=a} ^b {n \choose i} = \sum _{i=a} ^n {n \choose i} \textsf{ 当 } b>n \textsf{ 时} \tag{10} $

以下关系：

$ \binom{n}{0}=\binom{n}{n}=1,\ \binom{n}{1}=\binom{n}{n-1}=n\textsf{ 和 } \binom{n}{2}=\binom{n}{n-2}= \frac{n!}{2!\cdot(n-2)}=\frac{(n-1)\cdot n}{2} $

导致：

$ \sum _{i=0} ^0 \binom{n}{i}=1 \tag{11} $
  
$ \sum _{i=0} ^1 \binom{n}{i}=1+n \tag{12} $
$ \sum _{i=1} ^2 \binom{n}{i}=n+\frac{n \cdot (n-1)}{2}=\frac{n \cdot(n+1)}{2} \tag{13} $
$ \sum _{i=0} ^{n-1} \binom{n}{i}=2^n-1 \tag{14} $
  
$ \sum _{i=0} ^{n-2} \binom{n}{i}=2^n-n-1 \tag{15} $
$ \sum _{i=0} ^2 \binom{n}{i}= \sum _{i=n-2} ^n \binom{n}{i}=1+n+ \frac {n \cdot(n+1)}{2}=1+\frac {n \cdot (n+1)}{2} \tag{16} $

显然：

$ \sum _{i=a} ^a \binom{n}{i} = \binom{n}{i} \tag{17} $

恒等式：

$ \binom {n}{k}=\binom {n}{n-k} \textsf{ 当 } 0 \leqslant k \leqslant n \textsf{ 时} $

导致：

$ \begin{align} \sum _{i=a} ^b \binom {n}{i} &= \binom {n}{a} + \binom {n}{a+1} + \ldots \ + \binom {n}{b-1} + \binom {n}{b}\\ \\ &= \binom {n}{n-a} + \binom {n}{n-a-1} + \ldots \ + \binom {n}{n-b+1} + \binom {n}{n-b}\\ \\ &= \binom {n}{n-b} + \binom {n}{n-b+1} + \ldots \ + \binom {n}{n-a-1} + \binom {n}{n-a}\\ \\ &= \sum _{j=n-b} ^{n-a} \binom{n}{j} \textsf{ 当 } 0 \leqslant a \leqslant b \leqslant n \textsf{ 时} \end{align} \tag{18} $

这在 \( a > n - b \) 的情况下，允许我们“镜像”求和的范围，从而在计算中使用较小的数字。

代码

以上所有内容都实现在了 `Partial_Sum` 函数中。它以 x86-64 汇编语言编写，分别使用了 GNU As（基于 System V AMD64 ABI for Linux，文件名：Partial_Sum.s）和 MASM（基于 Windows 64 位调用约定，文件名：Partial_Sum.Asm）。

汇编文件由于操作系统的约定，寄存器使用方式有所不同，但通过恰当的寄存器选择，在某一点后它们是相同的。

它也以标准 C 语言（文件名：Partial_Sum.c）实现。汇编文件将 C 语句作为注释包含在内，因此这里只显示 C 代码。

该函数接收三个 `unsigned int` 类型的输入参数：

  • n：二项式的阶数
  • a：要相加的第一个系数的索引
  • b：要相加的最后一个系数的索引

返回值类型为 `unsigned long long`。

typedef unsigned int       DWORD;
typedef unsigned long long QWORD;

#define MAXULLEXP 63  /* Maximum ...

QWORD Partial_Sum( DWORD n, DWORD a, DWORD b ) {

    DWORD i,
          nCoeff;
    QWORD Result,
          Coeff,
          x,
          y,
          z;

首先，检查参数的有效性。如果 `a > n` 或 `a > b`，则返回 `0`。如果 `b > n`，则将 `b` 设置为 `n`，因为超出此范围的任何系数都为 `0` - 参见公式 10。

    if ( ( a > n ) | a > b )
        return 0;

    if ( b > n )
        b = n;

接下来，如果需要，将求和范围“镜像”。这不仅有助于处理较小的数字，而且还能简化简单情况的测试，后续将进一步说明。

    i = n - b;
    if ( a > i ) {
        b = n - a;
        a = i;
    }

然后，检查 `a == 0` 的情况，如果满足，则检查 `b` 的值：

  1. 如果 `b` 为 `0`，则只有一个系数，结果为 1，公式 12；
  2. 如果 `b` 为 `1`，则结果为 n + 1，即前两个系数的和，公式 13；
  3. 如果 `b` 为 `n`，则结果为所有系数之和，即 `2^n`；
  4. 如果 `b` 为 `n - 1`，则结果为 `2^n - 1`，公式 14；
  5. 如果 `b` 为 `n - 2`，则结果为 `2^n – n - 1`，公式 15。

如果求和范围已被“镜像”，则情况 2 也涵盖了 `a = n - 1`、`b = n` 的情况，情况 4 则涵盖了 `a = 1`、`b = n` 的情况。

此处，也检查 `n > 63` 的情况，因为这会导致溢出，因为 `unsigned long long` 的最大值为 `2^64 - 1`，此时返回 `0`。

    if ( !a )
        if ( !b )               /* sum _0 ^0 */  /* sum _n ^n */
            return 1;
        else if ( b == 1 )      /* sum _0 ^1 */  /* sum _n-1 ^n */
            return n + 1;
        else if ( b == n )      /* sum _0 ^n */
            if ( n > MAXULLEXP )
                return 0;
            else
                return ( 1 << n );
        else if ( b == n - 1 )  /* sum _0 ^(n-1) */
            if ( n > MAXULLEXP )
                return 0;
            else
                return ( ( 1 << n ) - 1 );
        else if ( b == n - 2 )  /* sum _0 ^(n-2) */
            if ( n > MAXULLEXP )
                return 0;
            else
                return ( ( 1 << n ) - n - 1 );

最后一个要检查的简单情况是第一个和最后一个索引是否相同，即 `a == b`，这意味着结果是单个系数的值，该值将被计算出来。这里也测试了一些简单情况——`a == 0`、`a == n` 和 `a == 1`，这也涵盖了 `a = n – 1`（因此 `b = n – 1`）的情况，因为在这种情况下，求和范围已被“镜像”。

    if ( a == b ) {
        if ( ( a == 0 ) | ( a == n ) )
            return 1;
        if ( a == 1 )  /* no need to check a == n – 1 */
            return n;
        Coeff = 1;
        for ( i = n - a + 1; i <= n ; i++ )
            Coeff *= i;
        x = 1;
        for ( i = 2; i <= a; i++ )
            x *= i;
        Coeff /= x;
        return Coeff;
    }

如果到达此处，则意味着必须按照公式 8 进行求和。首先，在汇编代码中，将要使用的非易失性寄存器保存起来。正确选择寄存器即可生成适用于 Linux/Unix 和 Windows 的相同代码。

计算要相加的系数的数量 `nCoeff`，该数量至少为 2，并设置两个辅助变量 `x` 和 `y`。

    nCoeff = b - a + 1;
    x = n - a;
    y = a;

然后，测试要相加的系数的数量是奇数还是偶数，代码将相应地进行处理。

如果要相加的系数的数量 `nCoeff` 是奇数，则测试另一个简单情况：第一个系数的索引是否为 `0`。如果是，则使用公式 11 和公式 13 计算值。

    if ( nCoeff & 1 ) {
        if ( !a ) {
            Coeff = n * ( n + 1 ) / 2;
            Result = Coeff + 1;
        }

否则，使用公式 2 计算第一个单独的系数，并将和初始化为该值。然后，使用公式 8 计算下一个两个系数的和，并将其加到总和中。

无论哪种情况，起始索引 `a` 都会增加，以补偿单个系数。

        else {
            if ( a == 1 )
                Coeff = n;
            else {
                z = 1;
                for ( i = 2; i <= y; i++ )
                    z *= i;
                Coeff = 1;
                for ( i = x + 1; i <= n ; i++ )
                    Coeff *= i;
                Coeff /= z;
            }
            Result = Coeff;
            Coeff *= ( ( n - a ) * ( n + 1 ) );
            Coeff /= ( ( a + 1 ) * ( a + 2 ) );
            Result += Coeff;
        }
        a++;
    }

如果 `nCoeff` 是偶数，再次测试 `a == 0` 的情况，如果是，则使用公式 12 计算第一项。此外，还测试 `a == 1` 的情况，如果为 `true`，则使用公式 13 计算；否则，使用公式 2 计算。最后，将和初始化为该值。

    else {
        if ( !a )
            Coeff = n + 1;
        else if ( a == 1 )
            Coeff = n * ( n + 1 ) / 2;
        else {
            y++;
            z = 1;
            for ( i = 2; i <= y; i++ )
                z *= i;
            Coeff = x + 1;
            for ( i = x + 2; i <= n + 1; i++ )
                Coeff *= i;
            Coeff /= z;
        }
        Result = Coeff;
    }

此时，剩余系数的数量是偶数。`nCoeff` 除以 2，只保留整数部分，并减 1 以补偿已计算的项。然后，使用公式 9 计算任何剩余偶数对系数的和，并将其加到结果中。

最后，在汇编代码中，恢复保存的寄存器，并将结果返回给调用例程。

    if ( nCoeff /= 2 )
        if ( nCoeff-- ) {
            x = n - a;
            y = a + 2;
            for ( i = 0; i < nCoeff; i++ ) {
                z = y++;
                z *= y++;
                Coeff *= x--;
                Coeff *= x--;
                Coeff /= z;
                Result += Coeff;
            }
        }

    return Result;
}

演示

文件 `Demo.c` 是一个用于测试该函数的驱动程序。顶部的 `#define` 设置了函数参数的值。

汇编/编译代码

Linux

在终端中，可以使用以下命令汇编函数：

as --64 -o Partial_Sum.o Partial_Sum.s

可以使用以下命令构建并运行演示：

gcc -o bps-a Demo.c Partial_Sum.s -z noexecstack
./Demo

Windows

在 Visual Studio Developer Command Prompt 中，可以使用以下命令汇编函数：

ml64 -c Partial_Sum.Asm

可以使用以下命令构建并运行演示：

cl /c Demo.c
ml64 -c Partial_Sum.Asm
link Demo.obj Partial_Sum.obj
Demo

当然，在任何情况下，您都可以添加任何所需的选项，例如优化或列表...

参考文献

历史

  • 2023 年 10 月 23 日：初版
© . All rights reserved.