How to check numerical precision of my computation of Stieltjes-constants?

Question

In a thread in MSE I proposed an older routine of mine for the efficient computation of coefficients; I use a very similar routine for the quick&dirty computation of the Stieltjes-constants.

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.

The only freely available table around seems to be that of S. Plouffe (linked via wikipedia) but they display only the first 78 numbers to 256 digits precision.

Update2: This is the effective formula to which the Pari/GP code reduces:

Let $ \qquad h_c = {1\over c!} \sum_{k=0}^\infty (-1)^k {\ln(1+k)^c\over1+k}$ This is done using the sumalt-procedure.

Next let $ \qquad r_c = - {\ln(2)^{c-1}\over c!} b_c$ where $b_c$ are the bernoulli numbers

Then $ \qquad \gamma_c = c! \sum_{d=0}^{c+1} h_d \cdot r_{c+1-d} $

So my question:

how could I possibly get an educated guess for the number of correct digits based on my Pari/GP-routine?*

Alternatively:

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

(here is the table with my current computations of 512 coeffs by 1000 digits)

---

Update1:
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.
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.
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...

---

Here is my routine. I reduced the precision-parameter so that this can just be copied & 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

Prepare computations with parameters for precision of computation

termsforseries = 32
digitstocompute = 200; digitstoshow = 12;
default(realprecision,digitstocompute)
default(format,Str("g0.",digitstoshow))
default(seriesprecision,termsforseries)

Compute the coefficients of the Laurent-expansion of the zeta by conversion from the same series-type of the eta-function (the alternating zeta)

\\ ========= 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

Answer 1

You could compare with the output from mpmath:

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

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.

Answer 2

I have recently computed a large table of rigorous values of the Stieltjes constants. Thanks to some coding by Jon Bober, the table can be browsed using a web interface on LMFDB.org (the L-functions and modular forms database):

http://beta.lmfdb.org/riemann/stieltjes/

For any $n \le 10^5$, the web interface allows printing $\gamma_n$ to at least 10000 digits (over 30000 digits for small $n$). The raw data (huge files) can also be downloaded.