How to calculate the summation $(\sum_{p = k}^{n} \binom{n}{p}) / 2^n$ quickly?

Solution 1:

$\newcommand{\bbx}[1]{\,\bbox[8px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ For 'large' $\ds{n, p\ \mbox{and}\ n - p}$, $\ds{{n \choose p} \sim 2^{n}\exp\pars{-\,{\bracks{p - n/2}^{2} \over n/2}}}$.

You can use the $\ds{\bbox[#dfd,5px]{\ Laplace\ Method\ for\ Sums\ }}$ ( see page $761$ in $\ds{\bbox[#fdd,5px]{\ Analytic\ Combinatorics\ }}$ by Philippe Flajolet and Robert Sedgewick, Cambridge University Press $2009$ ) \begin{align} {1 \over 2^{n}}\sum_{p = k}^{n}{n \choose p} & \sim {1 \over 2^{n}} \bracks{\int_{k}^{n}{n \choose n/2}\exp\pars{-\bracks{p - n/2}^{2} \over n/2}\,\dd p} \\[5mm] & \sim {1 \over 2^{n}}\,{n \choose n/2}{\root{2} \over 2}\,n^{1/2}\int_{\pars{k -n/2}/\root{n/2}}^{\infty}\exp\pars{-p^{2}}\,\dd p \\[5mm] & = {\root{2\pi} \over 4}\,{n \choose n/2}\,{n^{1/2} \over 2^{n}}\bracks{1 + \,\mrm{erf}\pars{n - 2k \over \root{2}\root{n}}} \quad \mbox{as}\ n \to \infty \end{align} where $\ds{\,\mrm{erfc}\pars{z} \equiv {2 \over \root{\pi}}\int_{0}^{z}\expo{-x^{2}}\,\dd x}$ is the Error Function.

Solution 2:

One method to compute the sum directly and without losing a lot of accuracy with finite precision arithmetic is to represent floating points as a pair of mantissa and exponent and to express factorials using the recurrence $$\binom{n}{r} = \binom{n}{r-1} \frac{n-r+1}{r}.$$

Here is a simple Python implementation for $\sum\limits_{p=0}^K \binom{N}{p}$:

def binomial_sum(N,K):
    current_exponent, current_mantissa = 0, 1.0
    total_exponent, total_mantissa = 0, 1.0
    for i in range(1, K+1):
            current_mantissa = (N-i+1)*current_mantissa/i
            while current_mantissa>=2:
                    current_mantissa/=2
                    current_exponent+=1
            while current_mantissa<=0.5:
                    current_mantissa*=2
                    current_exponent-=1
            total_mantissa += current_mantissa*pow(2,current_exponent-total_exponent)
            while total_mantissa>=2:
                    total_mantissa/=2
                    total_exponent+=1
    return total_exponent, total_mantissa

binomial_sum(10000000,5000000) is found to be 1.000252313246 × 2^9999999, correct to 12 decimal places.

Of course, this is not going to be as efficient as the approximations using integrals but it should work reasonably well for upto $N\approx10^8$ or so.