Some Constants from Number theory

How to compute them ?

(Click here for a Postscript version of this page.)

1  Constants and prime series

There are some important constants related to number theory which are defined by series involving sums or products on prime numbers (like Artin's Constant, the Twin prime Constant, Mertens' Constant, ...). Usually the direct estimation of such constants is extremely slow and only provides a few digits (say around 10). Some specific methods are available which allow a much greater accuracy (hundred or even thousand of decimal places !), in this essay we describe some of those algorithms.

In the following sections the letter p will always denote a prime number and the letter n an integer. The sums with p only take the primes in account.

2  Riemann zeta prime function

2.1  Definitions

The classical Riemann Zeta function is defined as

z(s)= ¥ å n=1   1 ns        Â(s) > 1,

where the sum is taken over all the positive integers and s is a complex number. And as showed by Leonhard Euler (1707-1783), we have the famous and most important relation

z(s)= Õ p prime  æ è  1 1-p-s ö ø ,

(1)

where, this time, only the prime numbers are taken in account.

Therefore, it's natural to define a new zeta function if instead of using all integers we restrict ourself to all the prime numbers.

Definition 1 The Riemann zeta prime function is defined as

zp(s)= ¥ å p=2   1 ps        Â(s) > 1,

(2)

where the sum is over all prime numbers p.

We have the easy majoration

zp(s) < ( z(s)-1) æ è 1-  1 2s ö ø .

The direct use of definition (2) is not convenient for any accurate numerical evaluation of z_p(s), because the convergence of the sum is logarithmic and therefore extremely slow. For instance, to compute 20 digits of z_p(2), you will need to compute all prime numbers up to about 10²⁰ which is not reasonable. The two following theorems provide a much better way to achieve this task.

Theorem 2 For Â(s) > 1,

log( z(s)) = ¥ å n=1   zp(sn) n .

Proof. Just take the logarithm of both side of Euler's relation (1)

log( z(s)) = - ¥ å p=2  log æ è 1-  1 ps ö ø = ¥ å p=2  æ è ¥ å n=1   1 npsn ö ø

thanks to -log(1-x)=x+x²/2+x³/3+..., then inverting the sum

log( z(s)) = ¥ å n=1   1 n æ è ¥ å p=2   1 psn ö ø = ¥ å n=1   zp(sn) n .

 

2.2  A useful relation

Theorem 3 For  (s) > 1,

zp(s)= ¥ å n=1   m (n) n log( z (sn))

(3)

where m (n) is the Möbius function defined by :

m(1) = 1, m(n) = 0, if n has a squared factor, m(n) = (-1)k, if n is the product of k different primes.

Proof. It's a consequence of a Möbius inversion formula applied on the previous theorem [6].  

Example 4 With s=2,

zp(2)=logz(2)-  1 2 logz(4)-  1 3 logz(6)-  1 5 logz(10)+...

The convergence of the series is geometric thanks to the relation

log( z (sn)) ~  1 2sn

and for s=2

log( z (2n)) ~  1 4n

therefore only

20 log10(4) » 33

terms in the series (3) are required to compute 20 digits of z_p(2).

There is a nice way to increase the speed of such a formula by computing the M first terms of the series z_p(s). Set

zp(s)= p £ M å p=2   1 ps + å p > M   1 ps

then the above method is easily generalized (left as exercise !, [1]) to

Theorem 5 For  (s) > 1,M ³ 2 a prime

zp(s)= p £ M å p=2   1 ps + ¥ å n=1   m(n) n log( z(M,sn))

(4)

with

z(M,s)=z(s) Õ p £ M  æ è 1-  1 ps ö ø .

The bigger M is the faster the rate of convergence is. The best choice for M depends on your implementation, and values around M=50 give good results.

Example 6 With s=2,M=7

zp(2)=  18589 44100 +logz(7,2)-  1 2 logz(7,4)-  1 3 logz(7,6)-  1 5 logz(7,10)+...

The convergence of the series is geometric thanks to the relation

log( z(7,sn)) ~  1 11sn

and for s=2

log( z(7,2n)) ~  1 121n

therefore, this time, less than

20 log10(121) » 10

terms in the series (4) are required to compute 20 digits of z_p(2).

2.3  Some values of z_p(s)

Using the previous theorems with relation (3) or (4) allows an easy and very fast computation of hundred digits of z_p(s), here are the first values computed with 60 digits :

ê ê ê ê ê ê ê ê ê ê ê ê ê zp(2)=0.452247420041065498506543364832247934173231343239892421736418... zp(3)=0.174762639299443536423113314665706700975412121926149289888672... zp(4)=0.076993139764246844942619295933157870162041059714843190264938... zp(5)=0.035755017483924257132818242538855711131697276726651331690092... zp(6)=0.017070086850636512954133673266059399209585941874544244733163... zp(7)=0.008283832856133592535124138729448723089183328885307806244641... zp(8)=0.004061405366517830560523439142683080522977144512071741001032...

In 1748, Euler made an estimation of the numbers z_p(s) up to z_p(36) but only for the even values of s with 15 digits of accuracy (around 13 digits were correct) [2].

In 1881, Merrifield computed all those numbers up to z_p(35) with 15 digits of accuracy [10]. In 1948, this computation was extended up to z_p(167) with 50 digits by Liénard [8], [13]. Nowadays it's possible to compute thousand of decimal places and for example P. Sebah found more than 10000 digits for z_p(2).

2.4  Applications

The values of numbers z_p(s) is useful to compute other important constants of number theory. For example in 1737, Euler proved the beautiful result that the sum of the inverse of the prime numbers is divergent (a proof is given in [13]) :

å p prime   1 p =¥,

but the divergence grows like log(log(m)). This suggest to define, just like Euler's constant, a new constant B₁:

B1= lim m® ¥  æ è å p £ m   1 p -log( log(m)) ö ø .

Unfortunately this definition is of no practical use to make any accurate computation of Mertens' constant B₁ and unlike Euler's Constant g, Euler-Maclaurin acceleration method cannot be used here. The following theorem due to Mertens in 1874 will produce a much better algorithm to find enough digits of B₁.

Theorem 7 (Mertens' 1874, [11])

Õ p £ m  æ è 1-  1 p ö ø ~  e-g log(m)        m > 0.

(5)

Proof. A proof may be found in [6].  

Taking minus the logarithm of both sides of Mertens' formula (5) gives

- å p £ m  log æ è 1-  1 p ö ø =g+log( log(m)) +em

(e_m is the remainder tending to zero with m) hence

g = lim m® ¥  æ è å p £ m  -log æ è 1-  1 p ö ø -log( log(m)) ö ø

and using the series expansion of -log(1-x)=x+x²/2+x³/3+...

g = lim m® ¥  æ è å p £ m  æ è  1 p +  1 2p2 +  1 3p3 +... ö ø -log( log(m)) ö ø

the two new formulae (the second is given in [6]) are then

B1 = g- ¥ å n=2   zp(n) n (6) B1 = g+ å p  æ è log æ è 1-  1 p ö ø +  1 p ö ø . (7)

2.4.1  Numerical value of B₁

The formula

B1=g-  zp(2) 2 -  zp(3) 3 -  zp(4) 4 -...

is convenient to evaluate B₁ because z_p(n) ~ 1/2ⁿ tends to zero at a geometrical rate. For example the first 100 digits may be computed in a few seconds only (in agreement with the 40 digits approximation given in [7]) and 8000 digits were obtained in a few hours by Sebah [5]. Here are its first decimal places :

B1 = 0.26149721284764278375542683860869585905156664 8261199206192064213924924510897368209714142631...

In [7] and [4], the following equivalent formula is given

B1=g+ ¥ å n=2   m(n) n log( z(n))

which avoid to compute then numbers z_p(n) in (6) by replacing them thanks to formula (3). This new relation may also be used for a verification of any numerical estimation of B₁.

2.4.2  Estimation of Artin's constant

Artin's constant occurs in the properties of prime numbers and is defined by

C= Õ p ³ 2   æ è 1-  1 p(p-1) ö ø

and therefore if j and y are the two roots of p²-p-1

log(C)= ¥ å p=2  log æ è  p2-p-1 p(p-1) ö ø = ¥ å p=2  log æ è  (1-j/p)(1-y/p) (1-1/p) ö ø

and as proved in [13], follows easily the relation

log(C)=- ¥ å n=2   (un-1) n zp(n)

with u₁=1,u₂=3 and u_n+1=u_n+u_n-1

log(C)=-zp(2)-zp(3)-  3 2 zp(4)-...

and Wrench gave in 1961, 45 digits of C [14]. Here are the first 60 digits

C=0.3739558136192022880547280543464164151116292486061500420947428...

The rate of convergence is much faster if, as suggested in [4], we compute the first terms of the product, for example :

C=  779 2016 Õ p ³ 11   æ è 1-  1 p(p-1) ö ø =  779 2016 D

and with the same method

log(D)=- ¥ å n=2   (un-1) n æ è zp(n)-  1 2n -  1 3n -  1 5n -  1 7n ö ø

with a convergence similar to (2/11)ⁿ.

2.4.3  Estimation of the twin primes constant

The twin primes constant (sometimes called Shall and Wilson constant [9]) is important in the estimation of density of twin primes (primes p such as p+2 is also prime, for example the couples (3,5),(5,7),(11,13),... are twin primes).

This constant is defined as [6]

C2= Õ p ³ 3   æ è 1-  1 (p-1)2 ö ø = Õ p ³ 3   æ è  p(p-2) (p-1)2 ö ø = Õ p ³ 3   æ è  (1-2/p) (1-1/p)2 ö ø

and just like for the Artin's constant

log(C2)= ¥ å n=2   (2-2n) n æ è zp(n)-  1 2n ö ø ,

therefore the first 60 digits are (Wrench also gave 45 digits [14] and Sebah extended the computation to 5000 decimal places [5]) :

C2=0.660161815846869573927812110014555778432623360284733413319448...

This method or some improvements [1] can be applied to compute the Hardy-Littlewood constants which occur in the density of prime constellations (extension of twin primes like prime triplets, prime quadruplets, ...). H. Cohen [1] also describes some similar algorithms to compute series of the form

å p ³ 2   1 pmlog( p) ,       m ³ 1.

Recently (1999), P. Moree and G. Niklasch have evaluated many of those constants up to 1000 decimal of accuracy. Their work was about some constants defined by infinite products of the form :

C= Õ p ³ 2   æ è 1-  f(p) g(p) ö ø

with f(p) and g(p) being polynomials with integer coefficients [12].

3  Derivative of the Riemann zeta prime function

Again, using Euler's relation (1) and taking the logarithm

log( z(s)) = - ¥ å p=2  log æ è 1-  1 ps ö ø

and if we take the derivative of this relation

z¢(s) z(s) =- ¥ å p=2   log(p)p-s 1-p-s

hence follows the theorem

Theorem 8 If  (s) > 1,

¥ å p=2   log(p) ps-1 =-  z¢(s) z(s) .

(8)

We suppose that we know how to evaluate z(s) and z^¢(s) (it's a classical problem and Euler-Maclaurin formula may be used)

z¢(s)=- ¥ å n=2   log(n) ns .

3.1  Application

The following constant l occurs in the estimation of the mean of the inverse Euler's totient function f(n)

l = å p ³ 2   log(p) p2-p+1 .

The direct use of this formula imposes to compute too many prime numbers to have a good estimation of l (the difficulty is similar to the brute estimation of z_p(2)), one way to avoid such huge calculation is to use a similar method as Kummer's acceleration process (see the essay Acceleration of the convergence of series in [5]). For instance

1 p2-p+1 -  1 p2-1 =  p-2 (p2-p+1)(p2-1) =O æ è  1 p3 ö ø

and the definition formula for l becomes

l = å p ³ 2  log(p) æ è  p-2 (p2-p+1)(p2-1) +  1 p2-1 ö ø

hence using (7)

l = å p ³ 2  æ è  (p-2)log(p) (p2-p+1)(p2-1) ö ø -  z¢(2) z(2)

this series has a faster convergence. The process may be repeated a few times giving for example

l = - å p ³ 2  æ è  (p5+p3-1)(p2-1)(p-1)log(p) (p8-1)(p6-1)(p5-1) ö ø -  z¢(2) z(2) -  z¢(3) z(3) +  z¢(4) z(4) +  z¢(5) z(5) +3  z¢(6) z(6) -  z¢(8) z(8)

the term of the new series tends to zero much faster (in O(1/p¹¹)), an estimation of the constant to 60 digits is

l = 0.608381717863324722683834585815620187759148598226022521199573...

and to achieve this accuracy all the prime numbers up to 10⁶ were computed (that is the first 78498 primes) which is reasonable.

4  Riemann zeta prime modulo functions

In this section we are interested in the evaluation of the Riemann zeta prime function but this time the primes taken in account are of the form 4k+1 or of the form 4k+3.

Definition 9 The Riemann zeta prime modulo functions are defined by

zQ(s) = ¥ å q=5   1 qs =  1 5s +  1 13s +  1 17s +...,       zR(s) = ¥ å r=3   1 rs =  1 3s +  1 7s +  1 11s +...       Â(s) > 1,

where q are the primes =1 mod 4 and r are the primes =3 mod 4.

In order to evaluate those functions with many digits we need to introduce the Dirichlet L functions.

4.1  The Dirichlet L functions

It will be convenient to use the following alternating series

L(s)= ¥ å n=0   (-1)n (2n+1)s =1-  1 3s +  1 5s -  1 7s +...       Â(s) ³ 1,

and note that

L(1) =  p 4 L(2) = G L(2m+1) =  | E2m| p2m+1 4m+1(2m)!

where G is Catalan's constant and E_2m are Euler's numbers. This new function is also sometime known as the Dirichlet Beta function.

And just like the Riemann zeta function, this function can also be expressed as a prime product in the form

L(s)= ¥ å n=1   c(n) ns = Õ p prime  æ è  1 1-c(p)p-s ö ø

where the character c is such as c(n)=0 if n is even and (-1)^(n-1)/2 when n is odd (see [6] for more details). Now, if we note q all the prime numbers of the form 4k+1 and r all the prime numbers of the form 4k+3 this product becomes

L(s)= Õ q prime=1 (4)  æ è  1 1-q-s ö ø Õ r prime=3 (4)  æ è  1 1+r-s ö ø .

4.2  Fast series for z_Q and z_R

We are now ready to derive a new expression for z_Q and z_R, if we rewrite Euler's product (1) as

z(s)= æ è  1 1-2-s ö ø Õ q prime=1 (4)  æ è  1 1-q-s ö ø Õ r prime=3 (4)  æ è  1 1-r-s ö ø

and we also have

z(2s) z(s) = Õ p prime  æ è  1-p-s 1-p-2s ö ø = Õ p prime  æ è  1 1+p-s ö ø = æ è  1 1+2-s ö ø Õ q prime=1 (4)  æ è  1 1+q-s ö ø Õ r prime=3 (4)  æ è  1 1+r-s ö ø .

This suggest, in order to separate the primes q and the primes r, to form the two new functions

a(s) = (1+2-s)-1  L(s)z(s) z(2s) b(s) = (1-2-s)  z(s) L(s)

and comes from this definition

a(s) = Õ q prime=1 (4)  æ è  1+q-s 1-q-s ö ø b(s) = Õ r prime=3 (4)  æ è  1+r-s 1-r-s ö ø

Theorem 10 For  (s) > 1,

1 2 log( a(s)) = ¥ å n=0   zQ((2n+1)s) 2n+1  1 2 log( b(s)) = ¥ å n=0   zR((2n+1)s) 2n+1

Proof. The proof is the same for the two identities and similar to the proof given in the first section for the Riemann zeta function. Just take the logarithm of the two previous relations and then use the series expansion

1 2 ( log(1+x)-log(1-x)) = x+  x3 3 +  x5 5 +...

for x=q^-s and x=r^-s respectively. The identities then follow easily.  

Theorem 11 For Â(s) > 1

zQ(s) =  1 2 ¥ å n=0   m(2n+1) 2n+1 log( a( (2n+1)s) ) (9) zR(s) =  1 2 ¥ å n=0   m(2n+1) 2n+1 log( b( (2n+1)s) ) (10)

where m(m) is the Möbius function.

Proof. It's again a consequence of the Möbius inversion formula applied on the previous theorem [6].  

Notice that the first relation of this theorem was given in [4].

By adding those expressions for z_Q(s) and z_R(s), we find a new form of the expression for the Riemann zeta prime function z_P(s)=1/2^s+z_Q(s)+z_R(s).

Corollary 12 For Â(s) > 1

zP(s)=  1 2 ¥ å n=0   m(2n+1) 2n+1 log( g( (2n+1)s) )

with

g(s)=  z2(s) z(2s) .

And for example comes the series

zP(2) =  1 2 ¥ å n=0   m(2n+1) 2n+1 log æ è  2(8n+4)!B4n+22 (4n+2)!2| B8n+4| ö ø =  1 2 log æ è  5 2 ö ø -  1 6 log æ è  715 691 ö ø -  1 10 log æ è  524875 523833 ö ø -...

where the B_m are Bernoulli's numbers and the rational number inside the log tends to 1 at a geometrical rate.

4.3  Some values of z_Q and z_R

Using the series (8) and (9) allow a fast computation of a few hundred digits and more of z_Q(s) and z_R(s).

4.3.1  z_Q table

Evaluation to 60 decimal places of

zQ(s)=  1 5s +  1 13s +  1 17s +...

ê ê ê ê ê ê ê ê ê ê ê ê ê zQ(2)=0.053813763574057670280678287341536562285675501495085532293911... zQ(3)=0.008755082732970504494226765813746675051112061220425472440026... zQ(4)=0.001649584154029291598996761313638851827487909943834732147811... zQ(5)=0.000323474034221797518511908186041083977442733705799147336695... zQ(6)=0.000064250963664773791101819138043576598984545546978815052898... zQ(7)=0.000012818448599795268251026582166507935820606749563344794362... zQ(8)=0.000002561371680396469808248432312473936447260601807298870666...

4.3.2  z_R table

Evaluation to 60 decimal places of

zR(s)=  1 3s +  1 7s +  1 11s +...

ê ê ê ê ê ê ê ê ê ê ê ê ê zR(2)=0.148433656467007828225865077490711371887555841744806889442507... zR(3)=0.041007556566473031928886548851960025924300060705723817448645... zR(4)=0.012843555610217553343622534619519018334553149771008458117126... zR(5)=0.004181543449702459614306334352814627154254543020852184353396... zR(6)=0.001380835886971739163031854128015822610601396327565429680264... zR(7)=0.000458514407533797266873112147282215153362722135744461450279... zR(8)=0.000152593994837434090715190710370606586529883910264442130365...

4.4  Application

Consider the function T(n) defined as the number of primes pairs of the form (i-1)²+1 and (i+1)²+1 lesser than n. Then when n becomes large

T(n) ~ A  Ön log2(n)

and the constant A is defined by

A=  p2 2 Õ q  æ è 1-  4 q ö ø æ è  q+1 q-1 ö ø 2   ,

with q being a prime of the form 4k+1 (this constant and many similar ones are detailed in Finch's excellent site on Mathematical constants [3]).

The direct use of this definition makes it very hard to find more than a few digits of this constant, but taking the logarithm of this expression gives

log(A)=log æ è  p2 2 ö ø + å q  log æ è 1-  4 q ö ø +2 å q  log æ è  q+1 q-1 ö ø

and expanding the two logarithms

log(A)=log æ è  p2 2 ö ø - å q  æ è å s ³ 1   4s s  1 qs ö ø +4 å q  æ è å s ³ 0   1 2s+1  1 q2s+1 ö ø

finally

log(A)=log æ è  p2 2 ö ø - å s ³ 2   4s s zQ(s)+4 å s ³ 1   1 2s+1 zQ(2s+1)

and comes the 50 digits evaluation in a few seconds of computation

A=1.9504911124462870744465855658095536925267084977189...

The convergence may be accelerated a lot if we compute the first terms of the product up to the prime q_s. For example with q_s=13, the formula becomes

log(A)=log æ è  441 1040 ö ø +log æ è  p2 2 ö ø - å s ³ 2   4s s zQ > 13(s)+4 å s ³ 1   1 2s+1 zQ > 13(2s+1)

with

zQ > 13(s)=zQ(s)-  1 5s -  1 13s ~  1 17s .

4.5  Other family of primes

It's possible to find with a similar method fast expressions for the Riemann zeta prime function restricted to primes of the form 3k+1,3k+2,6k+1 and 6k+5. The Dirichlet function must here be replaced by two other similar functions.

References

[1]

H. Cohen, High precision computation of Hardy-Littlewood constants, preprint, (1991)

[2]

L. Euler, Introduction à l'analyse infinitésimale (french traduction by Labey), Barrois, ainé, Librairie, (original 1748, traduction 1796), vol. 1

[3]

S. Finch, Favorite mathematical constants. World Wide Web site at http://pauillac.inria.fr/algo/bsolve/constant/constant.html, (1995)

[4]

P. Flajolet and I. Vardi, Zeta Function Expansions of Classical Constants, (1996)

[5]

X. Gourdon and P. Sebah, Numbers, Constants and Computation, World Wide Web site at the adress : http://numbers.computation.free.fr/Constants/constants.html, (1999)

[6]

G.H. Hardy and E. M. Wright, An Introduction to the Theory of Numbers, Oxford Science Publications, (1979)

[7]

D.E. Knuth, The Art of Computer Programming, Vol. II, Seminumerical Algorithms, Addison Wesley, (1998)

[8]

R. Liénard, Tables fondamentales à 50 décimales des sommes S_n,u_n,å_n, Paris, Centre de Docum. Univ., (1948)

[9]

F. Le Lionnais, Les nombres remarquables, Paris, Hermann, (1983)

[10]

C.W. Merrifield, The sums of the series of reciprocals of the prime numbers and of their powers, Proc. Roy. Soc. London, (1881), vol. 33, p. 4-10

[11]

F. Mertens, Journal für Math., (1874), vol. 78, p. 46-62

[12]

P. Moree, Approximation of singular series and automata, Manuscripta Math., (2000),  vol. 101, p. 385-399.

[13]

P. Ribenboim, The new Book of Prime Number Records, Springer, (1996)

[14]

J.W. Wrench, Evaluation of Artin's constant and the twin prime constant, Math. Comp., (1961), vol. 15, p. 396-398.