Implementation of nested summing in python

Question

I want to implement summation from this publication which is denoted as formula 29. As you may see there are 5 nested summations. Now I struggle to implement that. According to what I understand from my teacher, I should nest it in following way:

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

Now I by coeff{1,2,3,4} I mean those expressions after each sigma sign. I do it wrong however and I cannot tell where.

Could you give me a tip on that?

Best wishes!

The nested summation formula:

Answer 1

You approach should work. Here is what I came up with:

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

The try...except... block is present in testAll() because I don't know what the valid ranges are for mu and nu (and the range for m may also not be correct), so for some inputs B will attempt to compute the factorial of a negative number and throw an exception.

Let me know if you spot any errors...