I've recently been very interested in the wonderfully complex world of Euler sums, i.e. (convergent) infinite sums that, roughly speaking, consist of some rational polynomial combination of generalized harmonic numbers $H_p^{(q)} := \sum\limits_{k = 1}^{p}\frac{1}{k^q}$ and the parameters itself. For example $$\sum\limits_{k = 1}^{\infty} \frac{H_k}{k^2}, \sum\limits_{k = 1}^{\infty} \frac{H_k^3H_k^{(2)}}{k^2}, \sum\limits_{k = 1}^{\infty} \frac{H_k H_k^{(2)}H_k^{(4)}}{(k+1)^2}$$

are all convergent sums of this kind that are (with some exceptions) expressible as some polynomial combination of integer $\zeta$-values. Since I can't reasonably sink time into understanding all of the details of doing this for the cases I want to generalize, I turned to Mathematica in the hope that I could somehow squeeze out closed forms only dependent on zeta-values with some educated guessing for $\mathbb{Z}$-basis elements and the function FindIntegerNullVector.

I googled around a bit on how to apply this to my problem; aside from supplying sufficiently high precision of all involved numbers (done that up to 150 digits) and sorting the numbers to be related in certain ways, it seems to be wildly unreliable in outputting simple combinations in more complicated cases that Mathematica can't evaluate directly or are dependent on more than three combinations of $\zeta$-values.

For example:

Take the sum $\sum_{k = 1}^{\infty} \frac{H_k^3}{k^2} = 10\zeta(5) + \zeta(2)\zeta(3)$ - let's pretend we don't know this yet. Having worked with some Euler sums before, you might suspect that $\{\zeta(5),\zeta(2)\zeta(3),1\}$ is indeed a reasonable combination to express this sum as a linear combination of those, and in fact

In[*] := FindIntegerNullVector[N[{NSum[HarmonicNumber[k]^3/k^2,{k,1,Infinity},WorkingPrecision->20],Zeta[5],Zeta[2]Zeta[3],1},20]]
Out[*]= {1,-10,-1,0}

corresponding to the exact value of the sum if we take 20 digits of precision. Let's double the precision:

In[*] := FindIntegerNullVector[N[{NSum[HarmonicNumber[k]^3/k^2,{k,1,Infinity},WorkingPrecision->40],Zeta[5],Zeta[2]Zeta[3],1},40]]
Out[*]= {6293450745760373317489774,-62934507457603733174897742,-6293450745760373317489777,8} 

Ouch. I know that some sort of "overfitting" can occur when using Integer Relation algorithms, but how could one prevent this from happening? To make matters worse, let's look at another example.

Take the sum $\sum\limits_{k = 1}^{\infty} \frac{H_k^3}{k^6} = \frac{521}{24}\zeta(9)+3\zeta(2)\zeta(7)-\frac{97}{8}\zeta(3)\zeta(6) - \frac{51}{4}\zeta(4)\zeta(5) + 2\zeta(3)^3$. Again, with some educated guessing you might suspect the individual irrational terms to be a combination of this sum, and running the same calculations again

In[*] := FindIntegerNullVector[N[{NSum[HarmonicNumber[k]^3/k^6,{k,1,Infinity},WorkingPrecision->20],Zeta[9],Zeta[2]*Zeta[7],Zeta[3]*Zeta[6],Zeta[4]*Zeta[5],Zeta[3]^3},20]]
Out[*]= {55,42,-72,10,60,-35} 


In[*] := FindIntegerNullVector[N[{NSum[HarmonicNumber[k]^3/k^6,{k,1,Infinity},WorkingPrecision->40],Zeta[9],Zeta[2]*Zeta[7],Zeta[3]*Zeta[6],Zeta[4]*Zeta[5],Zeta[3]^3},40]]
Out[*]= {15349009,-24000024,-3069275,2446852,-14189761,14812626} 

where coefficients of both of these results are very far from their exact value, i.e. some "simpleness" assumptions would also be fruitless in this case.

Now why do I suspect that it should be theoretically possible to find closed forms of this kind with some integer relation algorithm? I actually got my inspiration for this from Vladimir Reshetnikov over at Math.SE, who has explained his out-of-the-blue thousands-of-digits-correct conjectures of some integrals/sums based on these relation algorithms* - which makes the poor performance of this PSLQ-function in this rather simple case... bizarre.

*See here and here. I do not actually have a link to it, but I vaguely recall him using the term "PSLQ" and presenting a high-caliber approximation in the same post.

My two questions are:

  1. How can I avoid "overfitting" (unreasonably large coefficients) when choosing a precision for my problem?
  2. What is the trick -other than educated guessing of the irrationals involved- that people like Reshetnikov use to obtain these obscenely good closed forms of really, really complicated expressions? Would this also work for this problem concerning Euler sums? If not, are there alternatives to PSLQ better fit for it?
  • 2
    $\begingroup$ I’ll take a (wild) guess at the failing examples. If NSum was unable to deliver on 40 digits of precision then a PSLQ-type relation finder will likely falter due to incorrect lower digits. $\endgroup$ Feb 11, 2023 at 17:16
  • 1
    $\begingroup$ Re “What is the trick”, there might be some inspiration to be found in work by various subsets of Jon and Peter Borwein, David Bailey and Simon Plouffe. $\endgroup$ Feb 11, 2023 at 17:18
  • 1
    $\begingroup$ @DanielLichtblau Damn... I suspected some kind of numerical inaccuracy in NSum as well. It didn't fail to converge in these test examples. Weirdly enough, cranking up the WorkingPrecision far enough should fix problematic precision inaccuracies if I interpret the "failed to converge[...]with error estimate[...]" correctly (as it would imply that enough digits were correct for the wanted precision). Did I misinterpret this? Is there any way to squeeze more correct digits out of the sum? $\endgroup$
    – TheOutZ
    Feb 11, 2023 at 17:25
  • 1
    $\begingroup$ Maybe try computing NSum to 60 digits and compare to computation at 40. If, say, 32 agree, then use that amount in your null vector computation (so you would have N[…,32]). This is just a heuristic of course. $\endgroup$ Feb 11, 2023 at 17:34
  • 2
    $\begingroup$ Assuming you fiddled with NSum option settings, I for one have no other suggestions. Might want to contact Tech Support since this might indicate a bug in NSum (or else limitations in numeric evaluation of the summands). $\endgroup$ Feb 11, 2023 at 17:56

2 Answers 2


(* "13.2.1 for Mac OS X ARM (64-bit) (January 27, 2023)" *)


To get high precision with NSum, in addition to increasing the WorkingPrecision, you should use more than the default value for NSumTerms

Options[NSum, NSumTerms]

(* {NSumTerms -> 15} *)

Prepend[(table = Flatten[Outer[{#1, #2,
        NSum[HarmonicNumber[k]^3/k^2, {k, 1, Infinity},
          NSumTerms -> #1, WorkingPrecision -> #2] -
         (10 Zeta[5] + Zeta[2] Zeta[3])} &,
      {15, 30, 60}, {20, 40, 60}], 1]),
  {"#terms", "prec", "error"}] //
 Grid[#, Frame -> All] &

enter image description here

Prepend[(table = Flatten[Outer[{#1, #2,
        NSum[HarmonicNumber[k]^3/k^6, {k, 1, Infinity},
          NSumTerms -> #1, WorkingPrecision -> #2] -
         (521/24 Zeta[9] + 3 Zeta[2] Zeta[7] - 97/8 Zeta[3] Zeta[6] -
           51/4 Zeta[4] Zeta[5] + 2 Zeta[3]^3)} &,
      {15, 30, 60}, {20, 40, 60}], 1]),
  {"#terms", "prec", "error"}] //
 Grid[#, Frame -> All] &

enter image description here


Here is a way to get at least 50 decimal places. One can boost this amount modestly (I did not try to quantify the error term too closely and skip that here).

First sum 10000 terms directly.

ee = HarmonicNumber[k]^3/k^2;
ss = Series[ee, {k, Infinity, 12}];
sum10k = Sum[ee, {k, 10000}];
nsum10k = N[sum10k, 50]

(* Out[901]= 12.217613264921967590943295626274344756881149234350 *)

To estimate the rest we sum the series expansion termwise from 1 to infinity (this can be done symbolically). We'll add this to the value above. Then we brute-force sum the first 10000 terms of the series expansions and subtract that from the sum, to account for the fact that those terms were double-counted.

seriesterms = ss[[3]];
seriestermsums = 
  Table[Sum[seriesterms[[j]]/k^(1 + j), {k, 1, Infinity}], {j, 
seriestermsumstot = Total[N[seriestermsums, 50]]
(* Out[904]= 12.347098083612797705384266996819903504057847643218 *)

seriestermssumsto10k = 
  Table[Sum[seriesterms[[j]]/k^(1 + j), {k, 1, 10000}], {j, 
totalsum10k = Total[N[seriestermssumsto10k, 50]]

(* Out[908]= 12.218129446803769914816822317038781008286036215889 *)

approxsum = nsum10k + seriestermsumstot - totalsum10k

(* Out[910]= 12.346581901730995381510740306055467252652960661679 *)

Compare to the proposed value.

N[10*Zeta[5] + Zeta[2]*Zeta[3], 50]

(* Out[867]= 12.346581901730995381510740306055467252652960661679 *)

I would expect that this approach might get one 100 digits without too much trouble. Maybe more, but I think 1000 would be out of reach and require a finer analysis and different (read: better) ideas.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.