# How to calculate an infinite sum to 100 exact digits with NSum?

In the discussion https://math.stackexchange.com/a/3419778/198592 I stumbled of the question how to calculate the sum

$$s= \sum _{n=3}^{\infty } \frac{n \cot \left(\frac{\pi }{n}\right)}{4^{n-2}}$$

to high precision.

I used NSum[] and played around with WorkingPrecision and Method, and the best I could get was 29 digits:

NSum[(n Cot[\[Pi]/n])/4^(n - 2), {n, 3, \[Infinity]}, WorkingPrecision -> 29,
Method -> "WynnEpsilon"]

(* Out[123]= 0.85222988130298006255574122938 *)


Increasing WorkingPrecision did not help as can be see from

sn[w_] := NSum[(n Cot[\[Pi]/n])/4^(n - 2), {n, 3, \[Infinity]},
WorkingPrecision -> w, Method -> "WynnEpsilon"]

Table[{w, Length[RealDigits[sn1[w]][[1]]]}, {w, 24, 32}]

(* Out[126]= {{24, 24}, {25, 24}, {26, 26}, {27, 26}, {28, 27}, {29, 29}, {30, 28}, {31,
28}, {32, 28}} *)


How can I get more exact digits, say 50 or 100?

I haven't had much success with NSum, and would suggest to sum terms until you reach the desired accuracy.

If you sum terms from $$n=3$$ to $$m$$ (for some large $$m$$), then you can get an estimate of the error by summing a series-expansion of the remaining terms. For large $$n$$, your term is approximately (16*n^2/π - 16*π/3)/(4^n).

Example: sum exactly up to $$m=1000$$,

With[{m = 1000},
N[Sum[4^(2-n)*n*Cot[π/n], {n, 3, m}], 10^3]]
(*    0.85222988130298006255...    *)


The error is approximately

With[{m = 1000},
N[Sum[(16*n^2/π - 16*π/3)/(4^n), {n, m+1, ∞}], 10^3]]
(*    1.4825646546726929...*10^-596    *)


So with 1000 terms we've gotten almost 600 correct digits. That's roughly $$\ln(4)/\ln(10)=0.60...$$ digits per term.

• Thank you very much for the elegant solution. Nice measure "digits per term". – Dr. Wolfgang Hintze Nov 3 '19 at 22:24

To achieve the result with NSum, use more terms for the extrapolation to get improved accuracy, increase WorkingPrecision to handle rounding error:

n1 = NSum[(n Cot[\[Pi]/n])/4^(n - 2), {n, 3, \[Infinity]},
NSumTerms -> 150, WorkingPrecision -> 150, Method -> "WynnEpsilon"]
Precision[%]
(*
0.85222988130298006255574122919169970056195391003233107805879735766771\
1344129758954639414777305170965859236369534172722

117.849
*)

n2 = NSum[(n Cot[\[Pi]/n])/4^(n - 2), {n, 3, \[Infinity]},
NSumTerms -> 250, WorkingPrecision -> 200, Method -> "WynnEpsilon"]
Precision[%]
(*
0.85222988130298006255574122919169970056195391003233107805879735766771\
1344129758954639414777305170965859236369534172721688746076335155778061\
150336000423518013038870770528545386615821

180.328
*)


The higher-precision result agrees with the 117-digit result to all digits (when rounded to the same precision).

n1 - n2
(*  0.*10^-118  *)

• @ Michael E2 thank you for this very instructive answer to the question as requested in the OP (using NSum). I had thought about this parameter before but I did not try to play with it. – Dr. Wolfgang Hintze Nov 4 '19 at 13:54
• @Dr.WolfgangHintze You're welcome. It's less successful on the oscillatory sum of (Sin[k]/k)^7. (You probably need to go out to well over 500000 terms to get a stable result, since k^-7 equals 0.5*^-40 around k -> 570000.) – Michael E2 Nov 4 '19 at 14:07

My question has been fully answered by the solutions of @Roman and @ Michael E2.

Here I'd like to add a comment to this solution.

Sometimes we need really very many terms to reach the required number of exact digits.

Here's an example of an infinite sum which was discussed earlier (Possible bug in infinite sum Sum[(Sin[k]/k)^m,{k,0,∞}]) and can be calculated in closed form:

s=Table[Sum[(Sin[k]/k)^n, {k, 0, \[Infinity]}], {n, 7, 7}]

(* Out[200]= {1 + (-23040 + 129423 \[Pi] - 201684 \[Pi]^2 +
144060 \[Pi]^3 - 54880 \[Pi]^4 + 11760 \[Pi]^5 - 1344 \[Pi]^6 +
64 \[Pi]^7)/46080} *)


With fourty digits the result is

N[%, 40]

(* Out[201]= {1.302724350729152883809457286983788824390} *)


If we use Sum to calculate these 40 digits correctly we need no less than 340'000 terms:

N[1 + Sum[(Sin[k]/k)^7, {k, 1, 340000}], 40]

(* Out[229]= 1.302724350729152883809457286983788824390 *)


Only one term less changes the last digit from $$0$$ to $$1$$:

N[1 + Sum[(Sin[k]/k)^7, {k, 1, 340000-1}], 40]

(* Out[239]= 1.302724350729152883809457286983788824391 *)

• NSum does much better if you break (Sin[k]/k)^7 up into terms of Exp[m I k]/k^7, m = 1, 3, 5, 7. You need less than 1000 terms with Wynn-Epsilon. – Michael E2 Nov 4 '19 at 14:21
• Thank you for the interesting hint. I understand it as a statement about the series for $\operatorname{Li}(7,e^{m i})$ – Dr. Wolfgang Hintze Nov 5 '19 at 9:45
• Yes, that's a succinct way of putting it, but I was thinking that the numerical behavior of terms like Exp[3 I k]/k^7 is easier to extrapolate than the sine form. – Michael E2 Nov 5 '19 at 14:27