Accelerating convergence for some double sums

I am interested in the following general double sums, for integers $$a\geq 1$$ and $$b\geq 2$$,

$$Z(a,b) = \sum_{k,\ell \geq 0} \frac{2k+3}{\binom{k+2}{2}^a} \frac{2\ell+3}{(\binom{k+2}{2}+\binom{\ell+2}{2})^b},$$

which are converging very slowly. For these sums, there is also an alternative expression as an iterated integral in dimension $$a+b$$, similar to multiple zeta values.

I would like more particularly to find $$Z(2,2)$$ and $$Z(1,3)$$ with precision as large as possible. Using the naive summation, I could only obtain 4 decimal digits, namely $$Z(2,2) \simeq 4.7058$$ and $$Z(1,3) \simeq 1.6470$$.

It is known that $$2 Z(2,2) + 4 Z(1,3) = 16$$.

What would be a smart way to accelerate the convergence, in general and in the special case using maybe the previous formula ?

• Please use a high-level tag like "nt.number-theory". I added this tag now. May 17 at 14:21
• $2𝑍(2,2)+4𝑍(1,3)=16$ ??? I want a combinatorial proof! May 19 at 14:45

This is easy to do with PARI/GP. Here is my code

p(n) = binomial(n+2,2);
Y(k,b) = sumnum(l=0, (2*l+3)/(p(k)+p(l))^b,sumtable);
Z(a,b) = sumnum(k=0, (2*k+3)/p(k)^a*Y(k,b));
default(realprecision,57);
sumtable = sumnuminit();
print(2*Z(2,2)+4*Z(1,3))
/* 16.0000000000000000000000000000000000000000000000000000000 */


It takes 1.361 CPU seconds for 57 decimal places. The documentation has some information about the methods being used and there are several summation functions other than sumnum. For example, I originally used sumpos but sumnum is faster. Thanks to Henri Cohen for his comment to use sumnuminit to speed up the calculation of Y(k,b). Thanks to Jorge Zuniga for his comment that replacing sumnum with summonien is much faster. And especially thanks to the developers of PARI/GP who provided these fast numerical summation functions.

Some details are in the arXiv paper Gaussian Summation: An Exponentially Converging Summation Scheme by Hartmut Monien. (published in 2010 in Mathematics of Computation).

• Thanks, very nice to see that open-source pari can do that. May 18 at 8:13
• @F.C. The recent book "Numerical algorithms for number theory - Using PARI/GP" by Belabas and Cohen discusses numerical summation in chapter 4, it also provides open-source programs. The relevant programs, besides sumpos, include SumLagrange and SumSidi. May 18 at 17:06
• @Somos: when you have double sums or double integrals, it is preferable to initialize nodes and weights once and for all for the inner sum. So if you set tab=sumnuminit() and write Y(k,b)=sumnum(l=0,(2*l+3)/(p(k)+p(l))^b,tab) the result will be obtained much faster (3 times at 57 decimals): note that the sumnuminit() command must be performed AFTER setting the precision to 57. May 19 at 12:02
• @Somos. Monien summation is much faster. p(n) = binomial(n+1,2); Y(k,b) = sumnummonien(l=1, (2*l+1)/(p(k)+p(l))^b,sumtable); Z(a,b) = sumnummonien(k=1, (2*k+1)/p(k)^aY(k,b)); default(realprecision,57); sumtable = sumnummonieninit(); print(2*Z(2,2)+4*Z(1,3)) /** 15.9999999999999999999999999999999999999999999999999999999, time = 16 ms. **/ May 19 at 16:42
• @Somos. Reference: H. Monien Math. Comp. 79 (2010), 857-869 "Gaussian quadrature for sums: A rapidly convergent summation scheme" ams.org/journals/mcom/2010-79-270/S0025-5718-09-02289-3/… May 19 at 21:53

With Mathematica I can first sum the series over $$\ell$$ to get a closed-form expression in terms of polygamma functions, $$Z(2,2)= \sum_{k,\ell \geq 0} \frac{2k+3}{\binom{k+2}{2}^2} \frac{2\ell+3}{(\binom{k+2}{2}+\binom{\ell+2}{2})^2}$$ $$\qquad=\sum_{k\geq 0}\frac{16 (2 k+3) }{(k+1)^2 (k+2)^2 \sqrt{-4 k (k+3)-7}}\left(\psi ^{(1)}\left[\frac{1}{2} \left(3-\sqrt{-4 k (k+3)-7}\right)\right)-\psi ^{(1)}\left(\frac{1}{2} \left(\sqrt{-4 k (k+3)-7}+3\right)\right)\right],$$ and then obtain a high-precision result using the Wynn epsilon method.. To fifteen decimal places I find

$$Z(2,2)=4.705905644127174\cdots$$

NSum[Z22,{k,0,Infinity},NSumTerms->2000,WorkingPrecision->250,Method->"WynnEpsilon"]


As a test, I also computed $$Z(1,3)$$ to fifteen decimal places, which gave

$$Z(1,3)=1.647047177936412\cdots.$$

And thus I find

$$2Z(2,2)+4Z(1,3)=15.999999999999996\cdots$$ correct to fifteen decimal places.

Here is a little Maple. Note that using "sum" on the inside causes it to find an algebraic expression for the sum over $$\ell$$ and using "Sum" on the outside tells it to not try to sum that algebraically over $$k$$.

I'll request 120 digits then round to 100 for confidence but I didn't need to. Maple increases the working precision automatically, and this is confirmed by the fact that exactly the same answers are obtained if 100 digits are requested directly.

The documentation says that Levin's u-transform is used, with these references:

Fessler, T.; Ford, W.F.; and Smith, D.A. "HURRY: An acceleration algorithm for scalar sequences and series." ACM Trans. Math. Software, Vol. 9, (1983): 346-354.

Levin, D. "Development of non-linear transformations for improving convergence of sequences". Internat. J. Comput. Math, Vol. B3, (1973): 371-388.

evalf[120](Sum((2*k+3)/((k+2)*(k+1)/2)^2
* sum( (2*l+3)/((k+2)*(k+1)/2 + (l+2)*(l+1)/2)^2,l=0..infinity),
k=0..infinity)):
Z22 := evalf[100](%);
Z22 := 4.7059056441271748413982167408282763463686149628226287913\
07611234496885020388225363260883986399848630 + 0. I

evalf[120](Sum((2*k+3)/((k+2)*(k+1)/2)^1
* sum( (2*l+3)/((k+2)*(k+1)/2 + (l+2)*(l+1)/2)^3,l=0..infinity),
k=0..infinity)):
Z13 := evalf[100](%);
Z13 := 1.6470471779364125793008916295858618268156925185886856043\
46194382751557489805887318369558006800075685 + 0. I

evalf[100](2*Z22 + 4*Z13);
16.0000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000 + 0. I

• Nice, thanks a lot. Too bad that I cannot reproduce the strategy with my open-source tools, not able to perform the inner symbolic summation. Would you be so kind to provide me with the value of $Z(1,2)$ too ? May 17 at 15:56
• I found maple somewhere and got Z12 := 6.133848088861232674022375960972559526576210014050489627587538049554517487250070394823631665263642367 - 0. I May 17 at 16:33
• @F.C. Yes I agree with that value. May 18 at 1:11