Are the Stieltjes constants irrational? Some computer experiments
AAre the Stieltjes constants irrational?Some computer experiments
Krzysztof D. Ma´slanka , Marek Wolf Institute for the History of Science of the Polish Academy of Sciences,ul. Nowy ´Swiat 72, pok. 9, 00-330 Warsaw, e-mail: [email protected] Cardinal Stefan Wyszynski University, Faculty of Mathematics and Natural Sciences.College of Sciences,ul. W´oycickiego 1/3, PL-01-938 Warsaw, Poland, e-mail: [email protected]
Abstract
Khnichin’s theorem is a surprising and still relatively little known result.It can be used as a specific criterion for determining whether or not anygiven number is irrational. In this paper we apply this theorem aswell as the Gauss–Kuzmin theorem to several thousand high precision(up to more than 53000 significant digits) initial Stieltjes constants γ n , n = 0 , , ..., a r X i v : . [ m a t h . N T ] S e p Introduction
The famous zeta function ζ ( s ) discovered by L. Euler in 1737 and published in 1744 [3] asa function of real variable was investigated by G. F. B. Riemann in the complex domainin his famous memoir submitted in 1859 to the Prussian Academy [14]. It is defined as: ζ ( s ) = ∞ (cid:88) n =1 n s (cid:60) ( s ) > . (1)It is divergent in the most interesting area of the complex plane, i.e., in the so calledcritical strip 0 ≤ (cid:60) ( s ) ≤ s = 1 by means of the following contour integral: ζ ( s ) = Γ( − s )2 πi (cid:90) P ( − x ) s e x − dxx , (2)where the integration is performed along the following path P : (cid:54) (cid:45) (cid:39)(cid:38)(cid:24)(cid:25) (cid:45)(cid:27) .Till now dozens of integrals and series representing the ζ ( s ) function are known, forcollection of such formulas see for example the entry Riemann Zeta Function in [17] andreferences cited therein and [12].Another representation of this function is given by a power series where appear certainconstants γ n . These constants are essentially coefficients of the Laurent series expansionof the zeta function around its only simple pole at s = 1: ζ ( s ) = 1 s − ∞ (cid:88) n =0 ( − n n ! γ n ( s − n (3)Primary definition of these fundamental constants was found by Th. J. Stieltjes andpresented in a letter to Ch. Hermite dated June 23, 1885 [5, letters no. 71–74] γ n = lim m →∞ (cid:34)(cid:32) m (cid:88) k =1 (ln k ) n k (cid:33) − (ln m ) n +1 n + 1 (cid:35) . (4)(When n = 0 the numerator of the fraction in the first summand in (4) is formally 0 which is taken to be 1.)Effective numerical computing of the constants γ n is quite a challenge because theformulas (4) converge extremely slowly. Even when n = 0, which corresponds to the well-known Euler-Mascheroni constant γ , in order to obtain just 10 accurate digits one has2o sum up exactly 12366 terms whereas in order to obtain 10000 digits (which is indeedrequired in some applications) one would have to sum up unrealistically large numberof terms: nearly 5 · which is of course far beyond capabilities of the present daycomputers. However, various fast algorithms were found to efficiently compute specificvalue of the zeroth Stieltjes constant γ , i.e. the fundamental Mascheroni-Euler constant,see e.g. [16], [2]. For n > Mathematica . An efficient but rather complicatedmethod based on Newton-Cotes quadrature has been proposed by R. Kreminski in 2003[9]. Quite recently F. Johansson presented particularly efficient method [6].In the Appendix at the end of the present paper yet another method of computingStieltjes constants will be described which is perhaps not as efficient as Johansson’s ap-proach, yet it is by far more simple and it may be easily and quickly used in practicalcalculations for obtaining γ n up to n ∼ ∼ γ n with accuracies ranging from about 53000 significant digits ( γ ) to about24000 digits ( γ ). Having these numbers we intend to provide an argument in favor oftheir irrationality. Then we consider the question of their normality, as real expansionsin the base equal 10. Finally, in Sect. 3, we develop γ n ’s into continuous fractions andnext use the remarkable theorems due to Khinchin, L´evy and Gauss–Kuzmin. Obtainedresults support the common opinion that γ n are indeed irrational. Let us recall that a number r is normal in base b if each finite string of k consecutive digitsappears in this expansion with asymptotic frequency b − k . In the usual decimal base wehave that each digit 0 , , , . . . , r with limitingfrequency 0.1, each 2–digits string 00 , , . . . ,
99 appears with density 0.01. Having thefirst 5000 Stieltjes constants with accuracies as described earlier we checked that eachdigit 0 , , , . . . , ×
10 data points in one plot. In the Fig. 1 we employed the following artifice: thefrequency h n (0) of appearance of digit 0 in the Stieltjes constant γ n is plotted at x –axisvalue n with the y value 0 . − h n (0), i.e. the distance from the expected value 0.1, whichin this case of a = 0 should be around 0.1. In general, the frequency h n ( a ) of appearanceof digit a in the Stieltjes constant γ n is plotted with the y value a × . . − h n ( a )).We calculated also density of 100 strings of two digits 00 , , . . . ,
99 for all 5000 Stieltjesconstants γ n . Now the result consisted of half a million points, what is impossible torepresent on the plot. Instead, in the Table I we present for each pattern of digits ab themaximal difference between calculated frequency of appearance and the expected value of0.01 and the number n of the Stieltjes constant γ n for which this discrepancy appeared.The difference between the actual computed value of the frequency of two digits patternsand the expected value 0.01 was typically of a few percents.3ig.1 The plot of the differences between 0.1 and actual frequencies of digits 0 , , . . . , a is plotted at y value a × . Continued fractions often reveal various profound and unexpected properties of irrationalnumbers that are normally hidden in their traditional decimal (or other basis) notation,see e.g. [11].In this Section we are going to exploit three facts about the continued fractions:the existence of the Khinchin constant, Khinchin–L´evy constant and the Gauss–Kuzmindistribution, see e.g. [8, chapter III, § § § γ n . The paper [1] presents the regular continued fraction for theEuler’s–Mascheroni constant γ . Let r = [ a ( r ); a ( r ) , a ( r ) , a ( r ) , . . . ] = a ( r ) + 1 a ( r ) + 1 a ( r ) + 1 a ( r ) + . . . (5)be the continued fraction expansion of the real number r , where a ( r ) is an integer and alldenominators a k ( r ) (“partial quotients”) with k ≥ a k . Khinchin has proved [8], see4lso [15], that limits of geometrical means of a k ( r ) are the same for almost all real r :lim l →∞ (cid:0) a ( r ) . . . a l ( r ) (cid:1) l = ∞ (cid:89) m =1 (cid:26) m ( m + 2) (cid:27) log m ≡ K = 2 . . . . . (6)The Lebesgue measure of (all) the exceptions is zero and include rational numbers ,quadratic irrationals and some irrational numbers too, like for example the Euler con-stant e = 2 . . . . for which the limit (6) is infinity. Table 1
In the columns A, C, E and G the two digits patterns are given, in the columns B, D, Fand H the maximal differences between 0.01 and the frequency that a given pattern ab , a, b = 0 , , . . . , γ n , n = 1 , , . . . , A B C D E F G H00 2 . × −
25 1 . × −
50 2 . × −
75 1 . × −
01 2 . × −
26 2 . × −
51 2 . × −
76 2 . × −
02 2 . × −
27 2 . × −
52 2 . × −
77 2 . × −
03 2 . × −
28 2 . × −
53 2 . × −
78 2 . × −
04 1 . × −
29 2 . × −
54 1 . × −
79 2 . × −
05 1 . × −
30 1 . × −
55 2 . × −
80 1 . × −
06 1 . × −
31 2 . × −
56 2 . × −
81 2 . × −
07 2 . × −
32 2 . × −
57 1 . × −
82 1 . × −
08 2 . × −
33 2 . × −
58 1 . × −
83 2 . × −
09 2 . × −
34 1 . × −
59 2 . × −
84 2 . × −
10 2 . × −
35 2 . × −
60 1 . × −
85 2 . × −
11 2 . × −
36 2 . × −
61 1 . × −
86 1 . × −
12 1 . × −
37 2 . × −
62 1 . × −
87 2 . × −
13 1 . × −
38 1 . × −
63 2 . × −
88 2 . × −
14 2 . × −
39 2 . × −
64 2 . × −
89 2 . × −
15 2 . × −
40 2 . × −
65 1 . × −
90 1 . × −
16 1 . × −
41 2 . × −
66 2 . × −
91 2 . × −
17 1 . × −
42 2 . × −
67 2 . × −
92 2 . × −
18 2 . × −
43 1 . × −
68 2 . × −
93 1 . × −
19 2 . × −
44 2 . × −
69 1 . × −
94 1 . × −
20 1 . × −
45 1 . × −
70 2 . × −
95 2 . × −
21 2 . × −
46 1 . × −
71 1 . × −
96 1 . × −
22 2 . × −
47 2 . × −
72 1 . × −
97 2 . × −
23 2 . × −
48 1 . × −
73 1 . × −
98 1 . × −
24 2 . × −
49 1 . × −
74 1 . × −
99 2 . × − The constant K is called the Khinchin constant, see e.g. [4, § K ( r ; l ) = (cid:0) a ( r ) a ( r ) . . . a l ( r ) (cid:1) l (7)for a given number r are close to K we can regard it as an indication that r is irrational.We developed the fractional parts of Stieltjes constants (in Sect.2, investigating thenormality, we used the whole number, e.g. γ = 111670 . . . . and we5ig.2 The plot of maximal a k ( n ) for n = 1 , , , . . . , contfrac ( r, { nmax } ) which creates the row vector a ( r ) whose components are the de-nominators a k ( r ) of the continued fraction expansion of r , i.e. a = [ a ( r ); a ( r ) , . . . , a l ( r )]means that r ≈ a ( r ) + 1 a ( r ) + 1 a ( r ) + 1. . . 1 a l ( r ) (8)The parameter nmax limits the number of terms a nmax ( r ); if it is omitted the expansionstops with a declared precision of representation of real number r at the last significantpartial quotient: the values of the convergents P k ( r ) /Q k ( r ) P k ( r ) Q k ( r ) = a ( r ) + 1 a ( r ) + 1 a ( r ) + 1 a ( r ) + . . . + 1 a k (9)approximate the value of r with accuracy at least 1 /Q k [8, Theorem 9, p.9]: (cid:12)(cid:12)(cid:12)(cid:12) r − P k ( r ) Q k ( r ) (cid:12)(cid:12)(cid:12)(cid:12) < Q k ( r ) , (10)6ence when 1 /Q k is smaller than the accuracy of the number r the process stops.We checked that the PARI precision set to \ p 120000 digits is sufficient in the sensethat scripts with larger precision generated exactly the same results: the rows a ( γ n ) ob-tained with accuracy 140000 digits were the same for all n as those obtained for accuracy120000 and the continued fractions with accuracy set to 100000 digits had different de-nominators a k ( γ n ). The number of partial quotients a k varied from over 110000 for initialStieltjes constants to 48027 for γ , i.e. the value of l ( n ) was roughly 2 times the numberof digits in the expansion of γ n . However, there have been cases of extremely large valuesof partial quotients. The largest was a = 17399017050 for γ , marked by the redarrow at the top in Fig. 2.Fig.3 The plot of K l ( n ( l )) for n = 1 , , , . . . , K than 0 .
001 and 30 points closer to K than 0 . | K − K n ( l ( n )) | is 4 . × − and it occurred for the Stieltjes constant number n = 3235 (marked withthe red arrow), the smallest value of | K − K n ( l ( n )) | is 1 . × − and it occurred for γ .With the precision set to 120000 digits we have expanded each γ n , n = 1 , , . . . . = means “approximately equal”) γ n . = [ a ( n ); a ( n ) , a ( n ) , a ( n ) , . . . , a l ( n ) ( n )] ≡ a ( n ) (11)without specifying the parameter nmax , thus the length of the vector a ( n ) depended on γ n and it turns out that the number l ( n ) of denominators was contained between 53000for Stieltjes constants with index around 5000 and 110000 for gammas with smallest index n . The value of the product a a . . . a l ( n ) was typically of the order 10 for beginningStieltjes constants to 10 for the last γ n ’s. It means that, if these Stieltjes constants7re rational numbers P/Q then Q are larger then those big numbers, for justification seee.g. [8, Theorems 16, 17]. Next for each n we have calculated the geometrical means: K n ( l ( n )) = l ( n ) (cid:89) k =1 a k ( n ) /l ( n ) . (12)The results are presented in the Fig.3. Values of K n ( l ( n )) are scattered around the redline representing K . To gain some insight into the rate of convergence of K n ( l ( n )) wehave plotted in the Fig. 4 the number of sign changes S K ( n ) of K n ( m ) − K for each n when m = 100 , , . . . l ( n ), i.e. S K ( n ) = number of such m that ( K n ( m + 1) − K )( K n ( m ) − K ) < . (13)The largest S K ( n ) was 961 and it occurred for the γ and for 124 gammas there wereno sign changes at all. It is well known that the convergence to Khinchin’s constant isvery slow. In the Fig.4 for each γ n we present the closest to the Khnichin constant K value of the “running” geometrical means K n ( m ) = (cid:32) m (cid:89) k =1 a k ( n ) (cid:33) /m , m = 100 , , . . . , l ( n ) . (14)Fig.4 The number of sign changes S K ( n ) for each n , i.e. the number of such m that( K n ( m + 1) − K )( K n ( m ) − K ) < m were skipped—sign changes were detected for m = 100 , , . . . l ( n ))..8ig.5 The plot of the closest to the Khnichin constant K values of the “running”geometrical means K n ( m ).Let the rational P k /Q k be the n -th partial convergent of the continued fraction: P k Q k = [ a ; a , a , a , . . . , a k ] . (15)For almost all real numbers r the denominators of the finite continued fraction approxi-mations fulfill [8, chapter III, § k →∞ (cid:0) Q k ( r ) (cid:1) /k = e π /
12 ln 2 ≡ L = 3 . . . . (16)where L is called the Khinchin—L´evy’s constant [4, § P l ( n ) ( γ n ) /Q l ( n ) ( γ n ) be the l -th partial convergent of the continuedfractions (11) of γ n : P l ( n ) ( γ n ) Q l ( n ) ( γ n ) = a ( n ) . = γ n . (17)For each Stieltjes constant γ n we calculated the partial convergents P l ( n ) ( γ n ) /Q l ( n ) ( γ n )using the recurrence: P = a , Q = 1 , P = 1 + a a , Q = a P k = a k P k − + P k − , Q k = a k Q k − + Q k − , k ≥ . (18)9ext from these denominators Q l ( n ) ( γ n ) we have calculated the quantities L n ( l ( n )): L n ( l ( n )) = (cid:0) Q l ( n ) (cid:1) /l ( n ) , n = 1 , , . . . , . (19)The obtained values of L n ( l ( n )) are presented in the Fig.6. These values scatter around thered line representing the Khinchin—L´evy’s constant L and are contained in the interval( L − . , L + 0 . γ ( n ) for which there appear sign changes of L − L n ( m ) , m = 1 , , . . . , l ( n ). As inthe case of K n ( m ) Fig.7 presents the number of sign changes of the difference L n ( m ) − L of the denominator of the m -th convergent P m /Q m S L ( n ) = number of such m that ( L n ( m + 1) − L )( L n ( m ) − L ) < . (20)The maximal number of sign changes was 922 and appeared for the Stieltjes constant γ and there were 117 gammas without sign changes.Finally we looked into the distribution of the values of partial quotients a l ( n ). TheGauss–Kuzmin theorem [8, chapter III, §
15] asserts that the density d ( k ) of the denomi-nators a m , m = 1 , , . . . l , with the value k is given bylim l →∞ (cid:93) { m : a m = k } l = log (cid:16) k k (cid:17) (21)for almost all real numbers. In the Fig. 9 the results are presented for the Stieltjesconstants. In 1997 it was shown by one of the authors of the present note [10] (M.K.) that theRiemann zeta function may be expressed as ζ ( s ) = 1 s − (cid:20) A + (cid:16) − s (cid:17) A + (cid:16) − s (cid:17) (cid:16) − s (cid:17) A
2! + ... (cid:21) = (22)= 1 s − ∞ (cid:88) k =0 A k k ! k (cid:89) i =1 (cid:16) i − s (cid:17) = (23)= 1 s − ∞ (cid:88) k =0 Γ (cid:0) k + 1 − s (cid:1) Γ (cid:0) − s (cid:1) A k k ! s ∈ C \{ } (24)where A k = k (cid:88) j =0 ( − j (cid:18) kj (cid:19) (2 j + 1) ζ (2 j + 2) = (25)= 12 k (cid:88) j =0 (cid:18) kj (cid:19) (2 j + 1) (2 π ) j +2 B j +2 (2 j + 2)! (26)10ig.6 The plot of L n ( l ( n )) for n = 1 , , , . . . , L than 0 .
001 and 38 closer to L than 0 . | L − L n ( l ( n )) | is4 . × − and it occurred for the Stieltjes constant number l = 3235 (marked withthe red arrow), the smallest value of | L − L n ( l ( n )) | is 2 . × − and it occurred forthe Stieltjes constant number l = 3226.Here B n denotes the n th Bernoulli numbers. However, the particular choice of nodes in s = 2 , , , ... , albeit the most natural, is by no means the only one. One only requires thatthe prescribed points be strictly equally spaced. For the purpose of present calculationswe choose the following sequence of points:1 + ε, ε, ε, ... where ε is certain real, not necessarily small number.More precisely, define certain entire function ϕ as: ϕ ( s ) := ( s − ζ ( s ) s (cid:54) = 1together with ϕ (1) = 1 which stems from the appropriate limit. Then, instead of (22),we have ϕ ( s ) = ∞ (cid:88) k =0 Γ (cid:0) k − s − ε (cid:1) Γ (cid:0) − s − ε (cid:1) α k k !with α k = k (cid:88) j =0 ( − j (cid:18) kj (cid:19) ϕ (1 + jε ) (27)Note that coefficients α k depend on ε but we shall for simplicity drop temporarily thisdependence in notation. 11ig.7 The number of sign changes S L ( n ) for each n , i.e. the number of such m that( L n ( m + 1) − L )( L n ( m ) − L ) < m were skipped— signchanges were detected for m = 100 , , . . . l ( n )).As mentioned in the Introduction the Stieltjes constants are essentially coefficients ofthe Laurent series expansion of the zeta function around its only simple pole at s = 1: ζ ( s ) = 1 s − ∞ (cid:88) n =0 ( − n n ! γ n ( s − n (28)Now directly from (3) we have: γ n = ( − n n + 1 d n +1 ds n +1 ϕ ( s ) (cid:12)(cid:12)(cid:12)(cid:12) s − . Then, after some elementary calculations, we get the following useful result: γ n = ( − n n ! ε n +1 ∞ (cid:88) k = n +1 ( − k k ! α k S ( k, n + 1) (29)where S ( k, i ) are signed Stirling numbers of the first kind. Note that in the literature thereare different conventions concerning denotation and indices of Stirling numbers which canbe confusing. Here we shall adopt the following convention involving the Pochhammersymbol: ( x ) k ≡ Γ( k + x )Γ( x ) = k − (cid:89) i =0 ( x + i ) = ( − k k (cid:88) i =0 ( − i S ( k, i ) x i Denoting β nk ≡ ( − n + k n ! k ! S ( k, n + 1) ε n +1 L values of the “running”values of m (cid:112) Q n ( m ) , n = 0 , , , . . . , γ n = ∞ (cid:88) k = n +1 β nk α k (30)The summation over k starts from n + 1 since β nk ≡ k ≤ n . Accuracy of α isequal to accuracy of precomputed values of ϕ ( s ) in equidistant nodes. When k grows theaccuracy of consecutive α k quickly tends do zero. Thus there always exists certain cut-offvalue of k = k . Therefore the summation in (30) may be performed to this value: γ n = k (cid:88) k = n +1 β nk α k (31)(Numerical experiment confirm that adding more terms do not affect the value of thesum (31).) As pointed earlier ε need not to be small, however, choosing smaller ε greatlyaccelerates convergence of the series. However, it also turns out that smaller ε impliessmaller k . What is really important: All significant digits of γ n obtained from the finitesum (31) are correct.Of course, γ n eventually does not depend on ε although α k as well as the rate ofconvergence of (29) does. In fact series (29) converges for any value of ε > ε (cid:29)
1. On the other hand, the smaller ε the faster the rate of convergence. However, since α k also depends on ε , choosing smaller13ig.9 The plot of the density of partial quotients a l equal to k = 1 , , . . . ,
10 from top tobottom for first 5000 Stieltjes constants. In red are the values of (21) plotted. The y axis is logarithmic to move the plots apart.value for ε requires higher accuracy of precalculated values of ϕ ( s ) which in turn may bevery time consuming. Hence, an appropriate compromise in choosing ε is needed.Formula (29) is particularly suited for numerical calculations. As already pointedabove, one has to choose parameter ε in order to optimally perform the calculations.Typically the algorithm has three simple steps:1. Tabulating ϕ (1 + jε ) , j = 0 , , , ... This requires appropriate choosing of parameter ε (see below) and is most time consuming. The most convenient for this seems smallbut extremely efficient program PARI/GP which has implemented particularly optimalzeta procedure. The first of authors used Cyfronet ZEUS computer in Cracow, wherecalculating single value of ϕ ( s ) with 51000 significant digits requires about 13 minutes.Since this procedure may easily be parallelized therefore in order to compute 10000 valuesof ϕ
20 independent routines were performed (each calculating 500 values of ϕ ) which tooknearly one week.2. Calculating α k using (27) and the precomputed values.3. Calculating Stieltjes constants using (29).(Contrary to the above step 1 which requires a powerful computer, steps 2 and 3 canbe quickly performed on a typical PC.) Several properties concerning accuracies may beobtained experimentally. It should be stressed out that given α k calculating single γ n withaccuracy of about 50000 digits requires several minutes on a very modest PC machine. Acknowledgement : One of the authors (KM) would like to express his gratitude to14he Academic Computer Center Cyfronet, AGH, Cracow, for the computational grant of1000 hours under the PL-Grid project (Polish Infrastructure for Supporting Computa-tional Science in the European Research Space).
References [1] R. P. Brent. Computation of the regular continued fraction for Euler’s constant.
Mathematics of Computation , 31(139):771–777, Jul 1977.[2] R. P. Brent and E. M. McMillan. Some New Algorithms for High-Precision Compu-tation of Euler’s Constant.
Mathematics of Computation , 34(149):305–312, 1980.[3] L. Euler. Variae observationes circa series infinitas.
Commentarii academiae scien-tiarum Petropolitanae , 9:160–188, 1744.[4] S. R. Finch.
Mathematical Constants . Cambridge University Press, 2003.[5] C. Hermite and T. J. Stieltjes. Correspondance d’Hermite et de Stieltjes. vol. 1, (8novembre 1882 - 22 juillet 1889), 1905.[6] F. Johansson. Rigorous high-precision computation of the Hurwitz zeta function andits derivatives.
Numerical Algorithms , 69(2):253–270, Jul 2014.[7] J. B. Keiper. Power series expansions of Riemann’s ξ function. Mathematics ofComputation , 58(198):765–765, May 1992.[8] A. Y. Khinchin.
Continued Fractions . Dover Publications, New York, 1997.[9] R. Kreminski. Newton-Cotes integration for approximating Stieltjes (generalizedEuler) constants.
Mathematics of Computation , 72(243), Dec 2002.[10] K. Ma´slanka. A hypergeometric-like Representation of Zeta-function of Riemann,Cracow Observatory preprint no. 1997/60, 1997. posted at arXiv: math-ph/0105007 ,2001. " http://xxx.lanl.gov/abs/math/0105007”.[11] G. Martin. The unreasonable effectualness of continued function expansions. Journalof the Australian Mathematical Society , 77(3):305–320, 2004.[12] M. Milgram. Integral and Series Representations of Riemann’s Zeta Function andDirichlet’s Eta Function and a Medley of Related Results.
Journal of Mathematics ,2013:Article ID 181724, 2013.[13] PARI/GP version , 64 bits, 2019. available from http://pari.math.u-bordeaux.fr/ .[14] B. Riemann. ¨Ueber die Anzahl der Primzahlen unter einer gegebenenGr¨osse.
Monatsberichte der K¨oniglich Preußischen Akademie der Wissenschaftenzu Berlin.
Studia Mathematica , 12:74–79, 1951.[16] D. W. Sweeney. On the Computation of Euler’s Constant.
Mathematics of Compu-tation , 17(82):170, Apr 1963.[17] E. W. Weisstein.