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

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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