MathOverflow will be down for maintenance for approximately 3 hours, starting Monday evening (06/24/2013) at approximately 9:00 PM Eastern time (UTC-4).

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

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

-

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.

-
 Ah, thanks; I'll see whether I can find something about the implementation in mpmath. – Gottfried Helms Aug 26 at 17:14 The algorithm is described here: mpmath.googlecode.com/svn/trunk/doc/build/… Specifically, it does numerical integration. – Fredrik Johansson Aug 26 at 17:39 @Fredrik: Well, very nice, indeed! That made it immediately into my bookmarks; thank you very much! I'll see what I can take from it... – Gottfried Helms Aug 26 at 20:49 @Fredrik: I added some heuristic for the achievable precision to my original question. – Gottfried Helms Aug 26 at 22:28 @Fredrik: I've re-computed the 1000-digit version with prec 1600 and uploaded to the given link; the last entry matches now exactly your reference where the last digit is rounded upwards (you have 3 digits more). So I seem to have at least a rough estimate for the lower bound of computational requirements... – Gottfried Helms Aug 26 at 22:59