Numerical multidimensional infinite series

Another new feature in mpmath: support for multidimensional series in nsum. Some very simple examples (finite and infinite summations can be combined in any desired order):

>>> nsum(lambda i,j: 2**(-i-j), [0,inf], [0,inf])
4.0
>>> nsum(lambda i,j,k: 2**(-i-j-k), [0,inf], [0,inf], [0,inf])
8.0
>>> nsum(lambda i,j,k,l: i/2**j*(i+l)/2**k, [1,2], [0,inf], [1,inf], [2,3])
50.0
>>> nsum((lambda i,j,k,l: 1/(2**(i**2+j**2+k**2+l**2))), *([[-inf,inf]]*4))
20.5423960756379
>>> jtheta(3,0,0.5)**4
20.5423960756379

One could of course also make nested calls to nsum, but having a direct syntax is much more convenient. Nested evaluation is also usually inefficient for convergence acceleration, so this is not what nsum does internally. Instead, it combines all the infinite summations to a single summation over growing hypercubes. The distinction is very important for conditionally convergent series.

For example, nsum can now directly evaluate the Madelung constant in three dimensions (ignore=True is to ignore the singular term at the origin):

>>> nsum(lambda i,j,k: (-1)**(i+j+k)/(i**2+j**2+k**2)**0.5,
...     [-inf,inf], [-inf,inf], [-inf,inf], ignore=True)
-1.74756459463318

While this evaluation takes several seconds, it is somewhat remarkable that a precise value can be obtained at all. A better way to compute the Madelung constant is using the following rapidly convergent 2D series:

>>> f = lambda i,j: -12*pi*sech(0.5*pi*sqrt((2*i+1)**2+(2*j+1)**2))**2
>>> nsum(f, [0,inf], [0,inf])
-1.74756459463318
>>> mp.dps = 100
>>> nsum(f, [0,inf], [0,inf])
-1.74756459463318219063621203554439740348516143662474175815282535076504062353276
117989075836269460789

Another nice application is to evaluate Eisenstein series directly from the definition:

>>> tau = 1j
>>> q = qfrom(tau=tau)
>>> nsum(lambda m,n: (m+n*tau)**(-4), [-inf,inf], [-inf,inf], ignore=True)
(3.1512120021539 + 0.0j)
>>> 0.5*sum(jtheta(n,0,q)**8 for n in [2,3,4])*(2*zeta(4))
(3.1512120021539 + 0.0j)