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 zp(s), because the convergence of the sum is logarithmic and therefore extremely slow. For instance, to compute 20 digits of zp(2), you will need to compute all prime numbers up to about 1020 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+x2/2+x3/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 zp(2).

There is a nice way to increase the speed of such a formula by computing the M first terms of the series zp(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 zp(2).

2.3  Some values of zp(s)

Using the previous theorems with relation (3) or (4) allows an easy and very fast computation of hundred digits of zp(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 zp(s) up to zp(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 zp(35) with 15 digits of accuracy [10]. In 1948, this computation was extended up to zp(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 zp(2).

2.4  Applications

The values of numbers zp(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 B1:
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 B1 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 B1.

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

(em 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+x2/2+x3/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 B1

The formula
B1=g-  zp(2)

2
-  zp(3)

3
-  zp(4)

4
-...

is convenient to evaluate B1 because zp(n) ~ 1/2n 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 zp(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 B1.

2.4.2  Estimation of Artin's constant

Artin's constant occurs in the properties of prime numbers and is defined by
C=
Õ
p ³ 
æ
è
1-  1

p(p-1)
ö
ø

and therefore if j and y are the two roots of p2-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 u1=1,u2=3 and un+1=un+un-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)n.

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 ³ 
æ
è
1-  1

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

(p-1)2
ö
ø
=
Õ
p ³ 
æ
è
 (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 ³ 
æ
è
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 zp(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/p11)), an estimation of the constant to 60 digits is
l = 0.608381717863324722683834585815620187759148598226022521199573...

and to achieve this accuracy all the prime numbers up to 106 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 E2m 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 zQ and zR

We are now ready to derive a new expression for zQ and zR, 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 zQ(s) and zR(s), we find a new form of the expression for the Riemann zeta prime function zP(s)=1/2s+zQ(s)+zR(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 Bm are Bernoulli's numbers and the rational number inside the log tends to 1 at a geometrical rate.

4.3  Some values of zQ and zR

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

4.3.1  zQ 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  zR 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)2+1 and (i+1)2+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 qs. For example with qs=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 Sn,un,å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.



Back to Numbers, Constants and Computation



File translated from TEX by TTH, version 3.01.
On 27 Nov 2001, 17:53.