How to check numerical precision of my computation of Stieltjes-constants? - MathOverflow most recent 30 from http://mathoverflow.net 2013-05-25T02:06:58Z http://mathoverflow.net/feeds/question/105535 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/105535/how-to-check-numerical-precision-of-my-computation-of-stieltjes-constants How to check numerical precision of my computation of Stieltjes-constants? Gottfried Helms 2012-08-26T14:31:37Z 2012-08-27T06:22:01Z <p>In a thread in <a href="http://math.stackexchange.com/questions/186292/derivatives-of-the-riemann-zeta-function-at-s-0" rel="nofollow">MSE</a> I proposed an older routine of mine for the efficient computation of coefficients; I use a very similar routine for the quick&amp;dirty computation of the Stieltjes-constants. </p> <p>This motivated me to try to improve my earlier toy-computations to calculate now the first 512 Stieltjes to 1000 dec digits precision. I'm unable to estimate the number of correct digits by analytical arguments; at least wolframalpha allowed me to display StieltjesGamma[511] to 400 digits, which met my own computations. </p> <p>The only freely available table around seems to be that of S. Plouffe (linked via <a href="http://en.wikipedia.org/wiki/Stieltjes_constants" rel="nofollow">wikipedia</a>) but they display only the first 78 numbers to 256 digits precision. </p> <p><strong>Update2:</strong> This is the effective formula to which the Pari/GP code reduces: </p> <p>Let $\qquad h_c = {1\over c!} \sum_{k=0}^\infty (-1)^k {\ln(1+k)^c\over1+k}$ This is done using the <em>sumalt</em>-procedure. </p> <p>Next let $\qquad r_c = - {\ln(2)^{c-1}\over c!} b_c$ where $b_c$ are the bernoulli numbers </p> <p>Then $\qquad \gamma_c = c! \sum_{d=0}^{c+1} h_d \cdot r_{c+1-d}$ </p> <p><em>So my question:</em></p> <blockquote> <p>how could I possibly get an educated guess for the number of correct digits based on my Pari/GP-routine?* </p> </blockquote> <p>Alternatively: </p> <blockquote> <p>is there some table with comparable precision around such that I can at least check the match for the first m digits (where m should optimally go to 1000)?</p> </blockquote> <p><em>(here is the table with my current computations of <a href="http://go.helms-net.de/math/tables/stieltjes_512x1000.zip" rel="nofollow">512 coeffs by 1000 digits</a>)</em><br> <hr> <strong>Update1:</strong><br> Heuristically I find, that beginning with some precision, say $300$ dec digits at the first $\gamma_0$ , I simply lose one digit precision per step in the index, so in $\gamma_k$ are roughly $300-k$ digits correct, maybe a handful less.<br> For this I used differences when computed with precision $200,300,400,500,600,700$ from that with precision $800$, $\gamma_0$ had just nearly all leading digits constant, when precision was increased, so that was always correct to the full precision.<br> That would mean, that if I want $1000$ correct digits for $\gamma_{511}$ I need dec precision of (at least) $1550$ . Simple, if that is true...</p> <hr> <p>Here is my routine. I reduced the precision-parameter so that this can just be copied &amp; pasted to a Pari/GP-environment. For precision of 1000 dec digits and 512 coefficients this must be optimized due to exorbitant increase of stack and computation-time otherwise</p> <p>Prepare computations with parameters for precision of computation</p> <pre><code>termsforseries = 32 digitstocompute = 200; digitstoshow = 12; default(realprecision,digitstocompute) default(format,Str("g0.",digitstoshow)) default(seriesprecision,termsforseries) </code></pre> <p>Compute the coefficients of the Laurent-expansion of the zeta by conversion from the same series-type of the eta-function (the alternating zeta) </p> <pre><code>\\ ========= Zeta Laurent-expansion providing Stieltjes-coefficients ==== ps_eta = sumalt(k=0,taylor((-1)^k/(1+k)^(1-x),x)) tmp = Vec(1-2*2^(-(1-x))); tmp[1]=0; \\ make the first zero exact. this step is needed for \\ allowing the reciprocal of the powerseries ps_etatozeta=1/Ser(tmp) ps_zeta = ps_eta * ps_etatozeta \\ contains now the Stieltjes-coefficients tmp=Vec(ps_zeta);tmp=vector(#tmp-1,c,tmp[1+c]) \\ remove the first coefficient (at 1/x) sti = vector(#tmp,r,tmp[r]*(r-1)!) \\ extract Stieltjes-constants by mult with factorials </code></pre> http://mathoverflow.net/questions/105535/how-to-check-numerical-precision-of-my-computation-of-stieltjes-constants/105545#105545 Answer by Fredrik Johansson for How to check numerical precision of my computation of Stieltjes-constants? Fredrik Johansson 2012-08-26T15:46:08Z 2012-08-26T15:46:08Z <p>You could compare with the output from mpmath:</p> <pre><code>sage: import mpmath sage: mpmath.mp.dps = 1000 sage: %time mpmath.stieltjes(511) CPU times: user 123.17 s, sys: 0.02 s, total: 123.19 s Wall time: 123.40 s mpf('673581492593841075447052270498937988033439947306384442967711559788996269245614412865378751092398327114199475672304543519558074203937367354475627304841991475249868411079091195038704370379319922304314968920977080.4218186954910530966341150821211999689800345913062006500416130863993252444286525401536530609127800808358611180051913954061786113778487954768827917318185861285728540852470806490244553130800206629709991267757983837666355484638397085316115099902138453930569718675294835237821298508690226519561229169443578986238598614523990440226172962706436119188515904391443174279895106345752233034115379099381680958168062786627389335290431416199037643058641914376639305675292168558263653044141610653456719446309980037732502545489019580865593535176949757824659484296855986638635532332512794555243036229273521585906314889067495562018805980518215400448131311489588531760771126389926309367463577942595344292677230759234541824332220012416082001221662802813469321335808232095303910714771240349667445255785796410716571') </code></pre> <p>This only seems to agree with your result up to about 715 digits. Mathematica 8 agrees with mpmath, so presumably you will need to increase the precision in your algorithm.</p>