Как вычислить дорогостоящую сумму высокой точности в python?

Question

Моя проблема очень проста. Я хотел бы вычислить следующую сумму.

from __future__ import division 
from scipy.misc import comb 
import math 

for n in xrange(2,1000,10): 
    m = 2.2*n/math.log(n) 
    print sum(sum(comb(n,a) * comb(n-a,b) * (comb(a+b,a)*2**(-a-b))**m 
        for b in xrange(n+1)) 
       for a in xrange(1,n+1)) 

Однако питон дает RuntimeWarning: overflow encountered in multiply и nan как выход, и это также очень и очень медленно.

Есть ли разумный способ сделать это?

Answer 1

Причина, почему вы получаете NaNs это вы в конечном итоге оценки числа как

comb(600 + 600, 600) == 3.96509646226102e+359 

Это слишком большой, чтобы поместиться в число с плавающей точкой:

>>> numpy.finfo(float).max 
1.7976931348623157e+308 

прологарифмируем, чтобы избежать этого:

from __future__ import division, absolute_import, print_function 
from scipy.special import betaln 
from scipy.misc import logsumexp 
import numpy as np 


def binomln(n, k): 
    # Assumes binom(n, k) >= 0 
    return -betaln(1 + n - k, 1 + k) - np.log(n + 1) 


for n in range(2, 1000, 10): 
    m = 2.2*n/np.log(n) 

    a = np.arange(1, n + 1)[np.newaxis,:] 
    b = np.arange(n + 1)[:,np.newaxis] 

    v = (binomln(n, a) 
     + binomln(n - a, b) 
     + m*binomln(a + b, a) 
     - m*(a+b) * np.log(2)) 

    term = np.exp(logsumexp(v)) 
    print(term) 

Answer 2

Используйте шаблон Memoize. При этом, переопределять гребенку:

@memoized 
def newcomb(a, b): 
    return comb(a, b) 

И заменить все вызовы comb с newcomb. Кроме того, для незначительного улучшения удалите скобки. Если вы делаете явные списки, вы тратите время на их создание. Если вы удалите их, вы эффективно используете generator expressions.

Update:

Это не решит проблему nan, но делает это намного быстрее.

Для всех, кто не видит это быстрее, применяете ли вы декоратор memoize? На моей машине исходная функция занимает 29,7 с до 200, но только 3,8 с записью в память.

Что memoize делает, просто хранить все ваши вызовы comb в таблице поиска. Поэтому, если на более поздней итерации вы вызываете comb с теми же аргументами, которые были у вас в какой-то момент в прошлом, она не пересчитывает их - она ​​просто выглядит в таблице поиска.