1

I would like to compute the following sums. The problem is that the binomial coefficients are too large I think so it fails.

from __future__ import division  
import numpy as np 
from scipy.special import binom 
print [sum(binom(n,2*k)*np.sqrt(np.pi*k)**(-n/10) for k in xrange(1,int(n/2)+1)) for n in xrange(100000)]

Is there some way to approximate the answer?

Simd
  • 19,447
  • 42
  • 136
  • 271

2 Answers2

2

The problem you're facing is that scipy.special.binom is approximating the answer. You could try computing them exactly using scipy.misc.comb with the optional argument exact=True.

With exact=False it uses a gamma function for fast computation, but you can force it to calculate the coefficient explicitly with exact=True. In that case, it returns a python long.

For example:

In [1]: from scipy.misc import comb
In [2]: comb(1100, 600, exact=1)
Out[2]: 3460566959226705345639127495806232085745599377428662585566293887742644983083368677353972462238094509711079840182716572056521046152741092473183810039372681921994584724384022883591903620756613168264181145704714086085028150718406438428295606240034677372942820551517227766024953527980780035209056864110017856973033878393954438656320L

Additionally, you could try using something other than scipy: gmpy here and here.

Community
  • 1
  • 1
wflynny
  • 18,065
  • 5
  • 46
  • 67
1

Well, binom tops out pretty early:

from scipy.special import binom

binom(1019, 509)    # => 1.40313388415e+305
binom(1020, 510)    # => inf

What exactly is the calculation you are trying to perform?

Here's an reformulated version which shifts the values around a bit; we can find successive values of the binomial series for each n instead of recomputing from scratch each time, and I've pushed the power-of-a-sqrt into a single operation.

from math import pi

for n in xrange(100000):
    total = 0.

    binom = 1
    binom_mul = n
    binom_div = 1

    power = -0.05 * n
    for k in xrange(1, n // 2 + 1):
        binom = binom * binom_mul / binom_div
        binom_mul -= 1
        binom_div += 1

        total += binom * (pi * k) ** power

    print(total)
Hugh Bothwell
  • 55,315
  • 8
  • 84
  • 99
  • Mathematically it is just what I have in the question. The values should eventually tend to zero. It looks like I need some way to avoid calculating these massive binomial coefficients explicitly. – Simd Apr 09 '14 at 15:00