Sum of Reciprocals of Germain Primes

Samuel S. Wagstaff, Jr.

Center for Education and Research in Information Assurance and Security

and Department of Computer Sciences Purdue University West Lafayette, IN 47907-1398 USA

                     ssw@cerias.purdue.edu

Abstract A prime p is a Germain prime if 2p + 1 is also prime. We compute the sum of the reciprocals of the Germain primes and related sequences. Since we do not know whether there are infinitely many Germain primes, all we can do is bound the sum in an interval.

1 Introduction

In 1919, Brun [4] proved that the sum

B = \frac{1}{3} + \frac{1}{5} + \frac{1}{5} + \frac{1}{7} + \frac{1}{11} + \frac{1}{13} + \frac{1}{17} + \frac{1}{19} + \frac{1}{29} + \frac{1}{31} + ⋯

of the reciprocals of all twin primes either converges or is a finite sum. Various authors [17] have estimated Brun’s constant B. Sebah [16] computed the sum of the reciprocals of all twin primes < 10^{16} and found that B > 1.830. Klyve [8] showed that B < 2.347 as reported in Crandall and Pomerance [5]. Platt and Trudgian [14] used the number of twin primes < 4·10^{18} computed by Oliveira e Silva [13] and found that 1.840503 < B < 2.288513. Several authors [3, 2, 10, 11, 13, 14, 16, 17] agree that the most probable value of B is about 1.90216, but this estimate is not rigorous. These results show how little we know about the distribution of twin primes, even whether there are infinitely many of them. The twin primes are Sequence \underline{A001097} in the OEIS [12].

A prime p is a Germain prime if 2p + 1 is also prime. They are named after Sophie Germain, who proved nearly 200 years ago that the first case of Fermat’s Last Theorem is true when the exponent is a Germain prime. The Germain primes are Sequence \underline{A005384} in the OEIS. Computing the sum of the reciprocals of the Germain primes is similar to computing Brun’s constant B. First, it is easy to prove using a sieve that this sum either converges or is a finite sum. Second, we do not know whether there are infinitely many Germain primes, so the best lower bound for the sum is the partial sum up to the limit to which we can compute it directly. One can compute a most probable value for the sum just as for B.

Little extra work is required to study the primes p for which 2^{k}p + 1 is also prime for any fixed k ≥ 1. The same methods apply to the sum of the reciprocals of the primes p for which p + 2^{k} is also prime. See Lee and Park [9] for the case of p, p + 8 both prime, for example.

2 Notation and easy results

Throughout this work p denotes a prime. Let

c_{2} = ∏_{p>2} (1 − \frac{1}{(p − 1)^{2}} ) ≈ 0.6601618158468695739278121100

be the twin prime constant as computed by Harley (See Lee and Park [9]).

Let a and b be positive integers with gcd(a, b) = 1 and 2 | ab. Let S_{a,b} = {p : ap + b is prime\}. For real x > 0, let S_{a,b}(x) denote the sum of \frac{1}{p} over all primes p ≤ x with p ∈ S_{a,b}. Also let S′_{a,b}(x) be the sum of \frac{1}{p} + \frac{1}{ap+b} over all primes p ≤ x with p ∈ S_{a,b}. Clearly one sum is finite if and only if the other is finite. Let S_{a,b} = lim_{x→∞} S_{a,b}(x) and S′_{a,b} = lim_{x→∞} S′_{a,b}(x).

Thus S′_{1,2} = B, Brun’s constant, and

S_{2,1} = \frac{1}{2} + \frac{1}{3} + \frac{1}{5} + \frac{1}{11} + \frac{1}{23} + \frac{1}{29} + \frac{1}{41} + \frac{1}{53} + ⋯

is the sum of the reciprocals of the Germain primes. Also,

S_{4,1} = \frac{1}{3} + \frac{1}{7} + \frac{1}{13} + \frac{1}{37} + \frac{1}{43} + \frac{1}{67} + \frac{1}{73} + \frac{1}{79} + ⋯

and

S_{8,1} = \frac{1}{2} + \frac{1}{5} + \frac{1}{11} + \frac{1}{17} + \frac{1}{29} + \frac{1}{71} + \frac{1}{101} + \frac{1}{107} + · · · .

See Sequences \underline{A023212} and \underline{A023228} in the OEIS. The case of S_{1,8} was studied by Lee and Park [9].

We show first that the limits S_{a,b} and S′_{a,b} always exist. Let π_{a,b}(x) be the number of primes p ≤ x for which ap + b is also prime.

Theorem 1. Suppose gcd(a, b) = 1 and 2 | ab. Define S_{a,b}(x), S′_{a,b}(x), S_{a,b} and S′_{a,b} as above. Then the limits in the definitions of S_{a,b} and S′_{a,b} exist.

Proof. Since S_{a,b}(x) is the sum of positive terms we need only show that it is bounded. According to Theorem 3.12 of Halberstam and Richert [7], we have

π_{a,b}(x) ≤ \frac{cx}{(log x)^{2}} (1 + O (\frac{log log x}{log x})) .

for some constant c. Therefore,

S_{a,b}(x) = ∑_{p≤x,p∈S_{a,b}} \frac{1}{p} = ∑_{t=1}^{x} \frac{π_{a,b}(t) − π_{a,b}(t − 1)}{t}

≤ ∑_{t=2}^{x} (\frac{1}{t} − \frac{1}{t + 1}) π_{a,b}(t) (since π_{a,b}(1) = 0)

≤ ∑_{t=2}^{x} \frac{1}{t(t + 1)} \frac{ct}{(log t)^{2}} (1 + O (\frac{log log t}{log t})) = O(1)

as x → ∞. The proof for S′_{a,b} is identical.

Fix x_{0}. Since 1/p > 0 always, we have S_{a,b} ≥ S_{a,b}(x_{0}) for every (a, b). These are the best lower bounds for S_{a,b} that we know how to compute. Table 1 shows these lower bounds for a few values of x_{0}. We computed lower bounds for Brun’s constant for comparison, although a better lower bound is found in Platt and Trudgian [14]. The sums were computed in double precision, but only nine decimal places are shown in the table. Table 2 shows the corresponding values of π_{a,b}. The data for twin primes is exactly the same as Nicely [10] found.

The most probable value for S_{a,b} is based on the Hardy-Littlewood [6] (or see Bateman and Horn [1]) approximation to π_{a,b}(x). This heuristic estimate is

π_{a,b}(x) ≈ 2c_{2} ∫_{2}^{x} \frac{dt}{(log t)^{2}} ∏_{p>2, p|ab} \frac{p − 1}{p − 2}

for fixed coprime integers a, b. When ab is a power of 2, which we always assume in this work, the product over primes p > 2 is 1. Assuming this heuristic and using partial summation, one finds

S_{a,b} − S_{a,b}(x_{0}) ≈ 2c_{2} ∫_{x_{0}}^{∞} \frac{dt}{t(log t)^{2}} = \frac{2c_{2}}{log x_{0}} .

Thus, after computing S_{a,b}(x_{0}), the most probable value for S_{a,b} is S_{a,b}(x_{0})+2c_{2}/ log x_{0}. The constant 2c_{2} is replaced by 4c_{2} for S′_{1,2} because we add the reciprocals of two (consecutive) primes in that case. These values are shown in Table 3.

\begin{matrix}x_{0} & S′_{1,2}(x_{0}) & S_{2,1}(x_{0}) & S_{4,1}(x_{0}) & S_{8,1}(x_{0}) \\ 10^{2} & 1.330990366 & 1.268745760 & 0.654987903 & 0.898299886 \\ 10^{3} & 1.518032463 & 1.350207371 & 0.718127473 & 0.985924806 \\ 10^{4} & 1.616893558 & 1.395040698 & 0.758673876 & 1.021448995 \\ 10^{5} & 1.672799585 & 1.422243022 & 0.783466978 & 1.045080342 \\ 10^{6} & 1.710776931 & 1.440222094 & 0.800867822 & 1.061448591 \\ 10^{7} & 1.738357044 & 1.453223143 & 0.813394490 & 1.073447436 \\ 10^{8} & 1.758815621 & 1.463079361 & 0.822866962 & 1.082581870 \\ 10^{9} & 1.774735958 & 1.470771610 & 0.830302440 & 1.089777846 \\ 10^{10} & 1.787478503 & 1.476946485 & 0.836290746 & 1.095594911 \\ 10^{11} & 1.797904311 & 1.482013752 & 0.841220054 & 1.100393501 \\ 10^{12} & 1.806592419 & 1.486246659 & 0.845347507 & 1.104420686 \\ 2 · 10^{12} & 1.808931050 & 1.487387756 & 0.846461664 & 1.105509109 \\ 5 · 10^{12} & 1.811852563 & 1.488814135 & 0.847855248 & 1.106871419 \\ 10^{13} & 1.813943761 & 1.489835793 & 0.848854028 & 1.107848344 \end{matrix}

Table 1: Some values of S′_{1,2}(x) and S_{a,b}(x).

\begin{matrix}x_{0} & π_{1,2}(x_{0}) & π_{2,1}(x_{0}) & π_{4,1}(x_{0}) & π_{8,1}(x_{0}) \\ 10^{2} & 8 & 10 & 9 & 6 \\ 10^{3} & 35 & 37 & 31 & 34 \\ 10^{4} & 205 & 190 & 176 & 161 \\ 10^{5} & 1224 & 1171 & 1057 & 1019 \\ 10^{6} & 8169 & 7746 & 7422 & 7090 \\ 10^{7} & 58980 & 56032 & 53709 & 51464 \\ 10^{8} & 440312 & 423140 & 407198 & 392990 \\ 10^{9} & 3424506 & 3308859 & 3198946 & 3095744 \\ 10^{10} & 27412679 & 26569515 & 25773602 & 25030013 \\ 10^{11} & 224376048 & 218116524 & 212205881 & 206597032 \\ 10^{12} & 1870585220 & 1822848478 & 1777532673 & 1734464714 \\ 2 · 10^{12} & 3552770943 & 3464491157 & 3380477421 & 3300430590 \\ 5 · 10^{12} & 8312493003 & 8112446112 & 7921805792 & 7740062386 \\ 10^{13} & 15834664872 & 15462601989 & 15107599980 & 14768799353 \end{matrix}

Table 2: Some values of π_{a,b}(x_{0}).

\begin{matrix}x_{0} & S′_{1,2} & S_{2,1} & S_{4,1} & S_{8,1} \\ 10^{2} & 1.904399633 & 1.555450394 & 0.941692537 & 1.185004520 \\ 10^{3} & 1.900305309 & 1.541343794 & 0.909263896 & 1.177061229 \\ 10^{4} & 1.903598191 & 1.538393015 & 0.902026193 & 1.164801312 \\ 10^{5} & 1.902163292 & 1.536924875 & 0.898148831 & 1.159762196 \\ 10^{6} & 1.901913353 & 1.535790305 & 0.896436033 & 1.157016803 \\ 10^{7} & 1.902188263 & 1.535138753 & 0.895310100 & 1.155363045 \\ 10^{8} & 1.902167938 & 1.534755520 & 0.894543121 & 1.154258029 \\ 10^{9} & 1.902160239 & 1.534483751 & 0.894014581 & 1.153489986 \\ 10^{10} & 1.902160356 & 1.534287412 & 0.893631672 & 1.152935837 \\ 10^{11} & 1.902160541 & 1.534141868 & 0.893348170 & 1.152521616 \\ 10^{12} & 1.902160630 & 1.534030764 & 0.893131612 & 1.152204792 \\ 2 · 10^{12} & 1.902160522 & 1.534002492 & 0.893076401 & 1.152123846 \\ 5 · 10^{12} & 1.902160560 & 1.533968133 & 0.893009246 & 1.152025417 \\ 10^{13} & 1.902160571 & 1.533944198 & 0.892962433 & 1.151956749 \end{matrix}

Table 3: Most probable values of S′_{1,2} and S_{a,b}.

The consistency of the estimates in the columns instills confidence, although there is no proof that they are anywhere near the true values. One might consider “higher order” terms in this approximation, but Shanks and Wrench [17] tell why this is unlikely to provide a closer estimate; see Sections 3 and 5 of their paper.

3 Upper bound on the limit

The upper bound is based on Lemma 5 of Riesel and Vaughan [15], which we quote here using Inequality (3.20) from their proof, as was done by Klyve [8].

Theorem 2. Let a and b be integers with a > 0, b ≠ 0 and gcd(a, b) = 1. Let

R(x, a, b) = sup_{I} ∑_{p∈I,ap+b prime} 1,

where the supremum is taken over all intervals of length x. Suppose that L and C = C(L) are related by Table 4. Then, whenever x ≥ e^{L} we have

R(x, a, b) < (\frac{16c_{2}x}{(log x)(C + log x)} + 2\sqrt{x}) ∏_{p>2, p|ab} \frac{p − 1}{p − 2} .

\begin{matrix}L & C & L & C \\ 24 & 0.97 & 48 & 8.2054 \\ 25 & 2.31 & 60 & 8.302 \\ 26 & 3.40 & 82 & 8.3503 \\ 27 & 4.28 & 100 & 8.3708 \\ 28 & 5.00 & 127 & 8.3905 \\ 29 & 5.58 & 147 & 8.404 \\ 31 & 6.45 & 174 & 8.4102 \\ 34 & 7.24 & 214 & 8.4201 \\ 36 & 7.56 & 278 & 8.4301 \\ 42 & 8.04 & 396 & 8.44004 \\ 44 & 8.11 & 690 & 8.45001 \end{matrix}

Table 4: L and C of Theorem 2.

The product in Theorem 2 is 1 when ab is a power of 2, as we always assume. Klyve [8] found it more convenient not to have the 2\sqrt{x} term in the theorem. He gave an alternate version, which we quote here.

Theorem 3. Let a and b be coprime positive integers with ab a power of 2. Suppose that L and D = D(L) are related by Table 5. Then, whenever x ≥ e^{L} we have

π_{a,b}(x) < \frac{16c_{2}x}{(log x)(D + log x)} .

The same bound applies when π_{a,b}(x) is replaced by the count over any other interval of length x ≥ e^{L}.

In Table 5 we correct a tiny error in Klyve [8]. He had D = 2.30 when L = 25; it should be D = 2.296.

\begin{matrix}L & D & L & D \\ 24 & 0.95 & 48 & 8.20 \\ 25 & 2.296 & 60 & 8.30 \\ 26 & 3.39 & 82 & 8.35 \\ 27 & 4.27 & 100 & 8.37 \\ 28 & 4.99 & 127 & 8.39 \\ 29 & 5.57 & 147 & 8.40 \\ 31 & 6.44 & 174 & 8.41 \\ 34 & 7.23 & 214 & 8.42 \\ 36 & 7.55 & 278 & 8.43 \\ 42 & 8.03 & 396 & 8.44 \\ 44 & 8.10 & 690 & 8.45 \end{matrix}

Table 5: L and D of Theorem 3.

Note that the upper bound in Theorem 3 is independent of a and b, so long as ab is a power of 2.

Theorem 4. We have 1.4898 < S_{2,1} < 1.8027, 0.8488 < S_{4,1} < 1.1617 and 1.1078 < S_{8,1} < 1.4208.

Proof. The lower bounds come from Table 1 with x_{0} = 10^{13}.

Let a and b be positive integers with gcd(a, b) = 1 and 2 | ab. Let 0 < M < N. Stieltjes integration by parts yields

∑_{t=M}^{N} \frac{π_{a,b}(t) − π_{a,b}(t − 1)}{t} = \frac{π_{a,b}(N)}{N} − \frac{π_{a,b}(M)}{M} + ∫_{M}^{N} \frac{π_{a,b}(t)}{t^{2}}dt .  (1)

Now

S_{a,b} = S_{a,b}(x_{0}) + ∑_{t=x_{0}}^{∞} \frac{π_{a,b}(t) − π_{a,b}(t − 1)}{t} .

We will divide the interval [x_{0}, ∞) into segments with boundaries at e^{L} for L in Table 5 to take advantage of the constants D(L) in that table. We will use Stieltjes integration to bound S_{a,b}(N) − S_{a,b}(M) on each interval [M, N).

Use Theorem 3 to bound the integral above. If L, L′ are consecutive entries in Table 5, then the integral over [M, N) = [e^{L}, e^{L′}) becomes

∫_{e^{L}}^{e^{L′}} \frac{π_{a,b}(t)}{t^{2}}dt ≤ ∫_{e^{L}}^{e^{L′}} \frac{16c_{2}t}{t^{2}(log t)(D(L) + log t)} dt

= 16c_{2} ∫_{L}^{L′} \frac{ds}{s(s + D(L))} (s = log t)

= \frac{16c_{2}}{D(L)}(log s − log(s + D(L))|_{L}^{L′}

= \frac{16c_{2}}{D(L)} log (\frac{L′(L + D(L))}{L(L′ + D(L))}) .

Note that log x_{0} ≈ 29.993606. If L ≤ log x_{0} < L′ in Table 5, then L = 29, L′ = 31 and the first integral is bounded by

\frac{16c_{2}}{D(L)} log (\frac{L′(log(x_{0}) + D(L))}{log(x_{0})(L′ + D(L))}) ≈ 0.010262 .

The last integral is bounded by

\frac{16c_{2}}{8.45} log (\frac{8.45 + 690}{690}) ≈ 0.015216 .

The values of these upper bounds are shown in Table 6.

M N  Upper Bound  M N  Upper Bound
x_{0} e^{31} 0.0102619142161059 e^{100} e^{127} 0.0208935182022260
e^{31} e^{34} 0.0250838160282704 e^{127} e^{147} 0.0106594237333466
e^{34} e^{36} 0.0143023911370928 e^{147} e^{174} 0.0105916018816358
e^{36} e^{42} 0.0350847550483840 e^{174} e^{214} 0.0108704126795364
e^{42} e^{44} 0.0096317734983264 e^{214} e^{278} 0.0109807270912696
e^{44} e^{48} 0.0170051336041258 e^{278} e^{396} 0.0110369750969163
e^{48} e^{60} 0.0381486495943858 e^{396} e^{690} 0.0111777020604581
e^{60} e^{82} 0.0421828971523529 e^{690} ∞ 0.0152151240070501
e^{82} e^{100} 0.0212203397984129
Table 6: M, N and the upper bound for ∫_{M}^{N} π_{a,b}(t)t^{−2}dt.

The sum of all upper bounds gives

∫_{x_{0}}^{∞} \frac{π_{a,b}(t)}{t^{2}}dt < 0.3143472 ,

which is the same for every (a, b) with ab a power of 2. The sum of the terms π_{a,b}(N)/N − π_{a,b}(M)/M in Equation (1) telescopes to give −π_{a,b}(x_{0})/x_{0}. Adding these values to the lower bounds from Table 1 gives

S_{2,1} < 1.4898358 − 0.0015463 + 0.3143472 = 1.8026367,
S_{4,1} < 0.8488541 − 0.0015108 + 0.3143472 = 1.1616905, and
S_{8,1} < 1.1078484 − 0.0014769 + 0.3143472 = 1.4207187.

Better upper and lower bounds would result from using a larger x_{0}. One probably can derive slightly smaller upper bounds (with the same x_{0}) by assuming the ERH as Klyve [8] did for S′_{1,2}.

4 Acknowledgments

I thank Harold Diamond and Carl Pomerance for valuable correspondence concerning this research. I am grateful to Nick Egbert for noticing the typo in Klyve’s thesis [8].

References

[1] P. T. Bateman and R. A. Horn, A heuristic asymptotic formula concerning the distribution of prime numbers, Math. Comp. 16 (1962), 363–367.

[2] J. Bowman, Some computational results regarding the prime numbers below 2,000,000,000, Nordisk Tidskr. Informationsbehanding (BIT) 13 (1973), 242–244.

[3] R. Brent, Irregularities in the distribution of primes and twin primes, Math. Comp. 29 (1975), 43–56. Corrigendum, ibid. 30 (1976), 1183.

[4] V. Brun, La série 1/5 + 1/7 + 1/11 + 1/13 + 1/17 + 1/19 + 1/29 + 1/31 + 1/41 + 1/43 + 1/59 + 1/61 + ⋯, ou les dénominateurs sont ‘nombres premieres jumeaux,’ est convergente ou finie, Bull. Sci. Math. 43 (1919), 100–104, 124–128.

[5] R. Crandall and C. Pomerance, Prime Numbers: A Computational Perspective, Springer-Verlag, 2001.

[6] G. H. Hardy and J. E. Littlewood, Some problems of ‘partitio numerorum,’ III: on the expression of a number as a sum of primes, Acta Math. 44 (1923), 1–70. Reprinted in Collected Papers of G. H. Hardy, Clarendon Press, Vol. 3, 1966, pp. 561–630.

[7] H. Halberstam and H.-E. Richert, Sieve Methods, London Mathematical Society Monographs, No. 4. Academic Press, 1974.

[8] D. Klyve, Explicit bounds on twin primes and Brun’s constant, Ph.D. thesis at Dartmouth College, Hanover, NH, 2007. Available at https://collections.dartmouth.

edu/archive/object/dcdis/dcdis-klyve2007.

[9] H. S. Lee and Y. Y. Park, On the primes with p_{n+1} − p_{n} = 8 and the sum of their reciprocals, J. Appl. Math. Comput. 22 (2006), 441–452.

[10] T. R. Nicely, Enumeration to 10^{14} of the twin primes and Brun’s constant, Virginia J. Sci. 46 (1996), 195–204.

[11] T. R. Nicely, A new error analysis for Brun’s constant, Virginia J. Sci. 52 (2001), 45–55.

[12] OEIS Foundation Inc., The On-Line Encyclopedia of Integer Sequences, 2021. Available

at https://oeis.org.

[13] T. Oliveira e Silva, Tables of values of pi(x) and of pi2(x), 2015. Available at http://

sweet.ua.pt/tos/primes.html.

[14] D. Platt and T. Trudgian, Improved bounds on Brun’s constant. Springer Proceedings in Mathematics and Statistics 313 (2020), 395–406. https://doi.org/10.

1007/978-3-030-36568-4_25. Also available at http://arxiv.org/abs/1803.01925,

March 5 2018.

[15] H. Riesel and R. C. Vaughan, On sums of primes, Ark. Mat. 21 (1983), 46–74.

[16] P. Sebah and X. Gourdon, Introduction to twin primes and Brun’s constant com-

    putation, 2002. Computation reported at http://numbers.computation.free.fr/
    Constants/constants.html.

[17] D. Shanks and J. W. Wrench, Jr., Brun’s constant, Math. Comp. 28 (1974), 293–299. Corrigendum, ibid., 28 (1974), 1183.

2010 Mathematics Subject Classification: Primary 11N05. Secondary 11L20, 11N37, 11N55.

Keywords: Brun’s constant, Germain prime, twin prime, sum of reciprocals.

(Concerned with sequences \underline{A001097}, \underline{A005384}, \underline{A023212}, and \underline{A023228}.)

Received August 16 2021; revised version received October 2 2021. Published in Journal of Integer Sequences, October 2 2021.

Return to Journal of Integer Sequences home page.