Download code and demo source - 6.6 KB

Binomial Coefficients

Binomial coefficients^[1] are the positive integers that occur as coefficients \( a_k \) in the expansion^[2] of the binomial:

$ (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} $

The coefficient \( a_k \), symbolized by \( \dbinom{n}{k} \) and read as "\( n \) choose \( k \)", is the number of combinations^[3] of \( n \) items taken \( k \) each time and is calculated by:

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

The presence of the factorials^[4] in it means many multiplications.

Dividing by \( (n-k)! \) both numerator and denominator, (1) is written as:

$ { 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} $

which is calculated faster as it uses less multiplications than (1).

Sum of Binomial Coefficients

The sum of binomial coefficients is given by:

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

but there is no formula^[5] for a partial sum of those, i.e., there is no formula for:

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

with \( 0 \leqslant a \leqslant b \leqslant n \). This means that one has to calculate all the coefficients from \( \dbinom{n}{a} \) to \( \dbinom{n}{b} \) inclusive, while adding on the way to a temporary result. Obviously, a slow process.

Two in One

Using Pascal's rule:

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

the sum (4) becomes:

$ \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} $

thus halving the number of calculations. In the case of odd number of coefficients, there remains a single coefficient to be calculated.

$ \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} $

Setting \( s=(b-a+1) \ mod\ 2 \) and \( t=\lfloor(b-a+1)/2\rfloor \) (5) and (6) can be written as:

$ \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} $

Calculation of the Sum of the First Two Coefficients

Using Pascal's rule, again:

$ { 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)!} $

Calculation of the Sum of the Next Two Coefficients

Once we have calculated the single coefficient, doing a little math gives:

$ \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} $

Also, after calculation of the sum of two coefficients, again with a little math, we get:

$ \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} $

In either case, the sum of the next two coefficients can be calculated quickly.

Some Simple Cases

For \( k>n \) it is \( \dbinom{n}{k} =0 \) thus

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

The relations:

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

lead to:

$ \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} $

Also, obviously:

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

The identity:

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

leads to:

$ \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{ for } 0 \leqslant a \leqslant b \leqslant n \end{align} \tag{18} $

This, in case \( a > n - b \), allows to "mirror" the range of the sum and thus deal with smaller numbers in the calculations.

The Code

All the above are realized in the function Partial_Sum. It is written in x86-64 Assembly in both GNU As, using the System V AMD64 ABI (Application Binary Interface) for Linux (Partial_Sum.s), and MASM using the Windows 64bit calling convention (Partial_Sum.Asm).

The Assembly files differ in register usage due to the OS conventions, but, with proper register selection, after a point, they are the same.

It is also written in standard C (Partial_Sum.c). The Assembly files contain the C statements as comments, so only the C code will be shown here.

The function takes three input parameters, of type unsigned int:

n the order of the binomial
a index of the first coefficient to be summed
b index of the last coefficient to be summed

The return value is of type 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;

First, the validity of parameters is checked. If a > n or a > b, 0 is returned. If b > n, then b = n, as any coefficients beyond that are 0 - eq. 10.

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

    if ( b > n )
        b = n;

Following this, the range of the sum is "mirrored", if needed. This, not only helps in working with smaller numbers, but also simplifies the tests for simple cases, as will be seen later.

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

Next, the case of a == 0 is checked and, if it is, b is checked for:

0, in which case there is only the first coefficient and the result is 1, eq. 12,
1, in which case the result is n + 1, the sum of the first two coefficients, eq. 13,
n, in which case the result is the sum of all coefficients, 2<sup>n</sup>, ,
n - 1, in which case the result is 2<sup>n</sup> - 1, eq. 14,
n - 2, in which case the result is 2<sup>n</sup> – n - 1, eq. 15.

If the range of sum has been "mirrored", case 2 also covers the case a = n - 1, b = n and case 4 covers the case a = 1, b = n.

Here, also the case of n > 63 is checked, as it causes overflow, because the maximum value of unsigned long long is 2<sup>64</sup> - 1 and 0 is returned.

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

Last of the simple cases to check is if the first and last indices are the same, a == b, which means the result is the value of a single coefficient, which is calculated. Here, too, some simple cases are tested – a == 0, a == n and a == 1 which also covers the case of a = n – 1 (and hence b = n – 1), as in this case, the range of the sum has been "mirrored".

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

Reaching here means that the sum must be calculated according to eq. 8. First, in Assembly, the non-volatile registers that will be used are saved. Proper selection of the registers results in the same code for Linux/Unix and Windows.

The number of the coefficients to be added, nCoeff, is found, which is at least 2, and 2 auxiliary variables, x and y, are setup.

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

Then the number of the coefficients to be added is tested for odd or even and the code advances accordingly.

If the number of the coefficients to be added, nCoeff, is odd, another simple case is tested: if the index of the first coefficient is 0. If so, the values are calculated using eq. 11 and 13.

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

Else, the first, single, coefficient is calculated using eq. 2 and the sum initialized to its value. Then, using eq. 8, the sum of the next two coefficients is calculated and added to the sum.

In either case, the starting index, a, is advanced, to compensate for the single coefficient.

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

If nCoeff is even, again the case a == 0 is tested, and, if so, the first term is calculated using eq. 12. Also, the case a == 1 is tested, and, if true, it is calculated using eq. 13. Else, it is calculated using eq. 2. Finally, the sum is initialized to this value.

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

At this point, the number of any remaining coefficients is even. nCoeff is divided by 2, keeping only the integer part, and decremented to comnpensate for the term that has already been calculated. Then, the sums of any remaining pairs of coefficients are calculated using eq. 9 and added to the result.

Finally, in the assembly code, the saved registers are restored, and the result returned to the calling routine.

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

A Demo

The file Demo.c is a driver program for testing the function. The #defines at the top set the values of the function parameters.

Assembling/Compiling the Code

Linux

At a terminal, the function can be assembled with the command:

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

The demo can be built and run with the commands:

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

Windows

At a Visual Studio Developer Command Prompt, the function can be assembled with the command:

ml64 -c Partial_Sum.Asm

The demo can be built and run with the commands:

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

Of course, in any case, you may add any option you want, such as optimization or listing ...

References

History

23^rd October, 2023: First version

Partial Sum of Binomial Coefficients