在 python 中实现嵌套求和
Implementation of nested summing in python
我想从这个 publication 中实现求和,表示为公式 29。如您所见,有 5 个嵌套求和。现在我很难实现它。根据老师的理解,我应该按照以下方式嵌套:
B=0.
for beta in range():
coeff1=...
sum1=0.
for beta_prim in range():
coeff2=...
sum2=0.
for alfa in range():
coeff3=...
sum3=0.
for alfa_prim in range():
coeff4=...
sum4=0.
for lambda in range():
coeff5=...
sum4+=coeff5
sum3+=sum4*coeff3
sum2+=sum3*coeff2
sum1+=sum2*coeff1
B+=sum1
return B
现在我用 coeff{1,2,3,4} 表示每个 sigma 符号后的那些表达式。
但是我做错了,我不知道在哪里。
你能给我一些建议吗?
祝福!
嵌套求和公式:
您的方法应该有效。这是我想出的:
from math import factorial
import scipy.misc
def comb(n,k):
return scipy.misc.comb(n,k,True)
def monepow(n):
if n % 2 == 0:
return 1
else:
return -1
def B(mu, nu, n, np, l, lp, m):
B = 0
beta1 = max( 0, mu - np - l, nu - np - l + m )
beta2 = min( n - l, nu )
for beta in range(beta1, beta2+1):
c1 = comb(n-l, beta)
beta1p = max( 0, mu - beta - l - lp, nu - beta - l - lp + m )
beta2p = min( np - lp, nu - beta )
for betap in range(beta1p, beta2p+1):
delta = min( mu + nu + l + lp, mu + n + np ) + 1
c2 = comb(np - lp, betap) * factorial(delta) / factorial (mu + beta + betap + l + lp + 1)
alpha1 = max( m, nu - beta - betap - lp + m )
for alpha in range(alpha1, l+1):
c3 = monepow(alpha) * comb(l, alpha) * comb(l, alpha - m)
alpha1p = max( m, nu - beta - betap - alpha + m )
for alphap in range(alpha1p, lp):
c4 = monepow(alphap) * comb(lp, alphap) * comb(lp, alphap-m) * comb(alpha - alphap - m, nu - beta - betap) * factorial(alpha - alphap + beta + lp) * factorial( alpha - alphap + beta + l)
for lam in range(0, mu+1):
c5 = monepow(lam) * comb( alpha - alphap + beta + lp + lam, lam ) * comb (alpha - alphap + betap + l + mu - lam, mu - lam) * comb(mu, lam)
x = c1*c2*c3*c4*c5
B += x
return B
def testAll():
for mu in range (0,10):
for nu in range (0,10):
for n in range(1,5):
for np in range(1,5):
for l in range(0,n):
for lp in range(0,np):
for m in range(0,l+1):
b = None
try:
b = B(mu, nu, n, np, l, lp, m)
except:
pass
if b is not None:
print [mu, nu, n, np, l, lp, m], " -> ", b
# print B(1, 1, 6, 6, 0, 4, 3) # should be == 0
# print B(1, 2, 6, 6, 4, 0, 0) # should be == 0
# print B(2, 2, 6, 6, 4, 0, 0) # might be == 0
# print B(3, 2, 6, 6, 4, 0, 0) # might be == 0
testAll()
try...except...
块出现在 testAll()
中,因为
我不知道 mu 和 nu 的有效范围是多少(以及 m 的范围也可能不正确),因此对于某些输入,B 将尝试计算负数的阶乘并抛出异常。
如果您发现任何错误,请告诉我...
我想从这个 publication 中实现求和,表示为公式 29。如您所见,有 5 个嵌套求和。现在我很难实现它。根据老师的理解,我应该按照以下方式嵌套:
B=0.
for beta in range():
coeff1=...
sum1=0.
for beta_prim in range():
coeff2=...
sum2=0.
for alfa in range():
coeff3=...
sum3=0.
for alfa_prim in range():
coeff4=...
sum4=0.
for lambda in range():
coeff5=...
sum4+=coeff5
sum3+=sum4*coeff3
sum2+=sum3*coeff2
sum1+=sum2*coeff1
B+=sum1
return B
现在我用 coeff{1,2,3,4} 表示每个 sigma 符号后的那些表达式。 但是我做错了,我不知道在哪里。
你能给我一些建议吗?
祝福!
嵌套求和公式:
您的方法应该有效。这是我想出的:
from math import factorial
import scipy.misc
def comb(n,k):
return scipy.misc.comb(n,k,True)
def monepow(n):
if n % 2 == 0:
return 1
else:
return -1
def B(mu, nu, n, np, l, lp, m):
B = 0
beta1 = max( 0, mu - np - l, nu - np - l + m )
beta2 = min( n - l, nu )
for beta in range(beta1, beta2+1):
c1 = comb(n-l, beta)
beta1p = max( 0, mu - beta - l - lp, nu - beta - l - lp + m )
beta2p = min( np - lp, nu - beta )
for betap in range(beta1p, beta2p+1):
delta = min( mu + nu + l + lp, mu + n + np ) + 1
c2 = comb(np - lp, betap) * factorial(delta) / factorial (mu + beta + betap + l + lp + 1)
alpha1 = max( m, nu - beta - betap - lp + m )
for alpha in range(alpha1, l+1):
c3 = monepow(alpha) * comb(l, alpha) * comb(l, alpha - m)
alpha1p = max( m, nu - beta - betap - alpha + m )
for alphap in range(alpha1p, lp):
c4 = monepow(alphap) * comb(lp, alphap) * comb(lp, alphap-m) * comb(alpha - alphap - m, nu - beta - betap) * factorial(alpha - alphap + beta + lp) * factorial( alpha - alphap + beta + l)
for lam in range(0, mu+1):
c5 = monepow(lam) * comb( alpha - alphap + beta + lp + lam, lam ) * comb (alpha - alphap + betap + l + mu - lam, mu - lam) * comb(mu, lam)
x = c1*c2*c3*c4*c5
B += x
return B
def testAll():
for mu in range (0,10):
for nu in range (0,10):
for n in range(1,5):
for np in range(1,5):
for l in range(0,n):
for lp in range(0,np):
for m in range(0,l+1):
b = None
try:
b = B(mu, nu, n, np, l, lp, m)
except:
pass
if b is not None:
print [mu, nu, n, np, l, lp, m], " -> ", b
# print B(1, 1, 6, 6, 0, 4, 3) # should be == 0
# print B(1, 2, 6, 6, 4, 0, 0) # should be == 0
# print B(2, 2, 6, 6, 4, 0, 0) # might be == 0
# print B(3, 2, 6, 6, 4, 0, 0) # might be == 0
testAll()
try...except...
块出现在 testAll()
中,因为
我不知道 mu 和 nu 的有效范围是多少(以及 m 的范围也可能不正确),因此对于某些输入,B 将尝试计算负数的阶乘并抛出异常。
如果您发现任何错误,请告诉我...