Enumeration to 1.6e15 of the twin primes and Brun's constant

Thomas R. Nicely

Current e-mail address

Freeware copyright (c) 2011 Thomas R. Nicely. Released into the public domain by the author, who disclaims any legal liability arising from its use.

Last revised 0800 GMT 05 February 2011.

Date of first release 16 July 1999.

The content of this document (other than the addendum, which was not part of the submission for publication) is essentially that of the original release, except that information rendered obsolete by subsequent events has been removed or modified, in both the main document and the addendum (this liberty is taken in view of the fact that the paper was never accepted for publication). There may also be differences in formatting, and in minor details and corrections, including updated URLs.


Enumeration of the twin primes, and the sum of their reciprocals, is extended to 1.6e15. An improved estimate is obtained for Brun's constant, B_2 = 1.90216 05824 +/- 0.00000 00030. Error analysis is presented to support the contention that the stated error bound represents a 99 % confidence level.

Mathematics Subject Classification 2000 (MSC2000)

Primary: 11A41.
Secondary: 11-04, 11Y70, 11Y60.

Key Words and Phrases

Twin primes, Brun's constant, prime numbers.

1. Introduction

The set K_2 = {(3, 5), (5, 7), (11, 13), (17, 19), ...} of twin-prime pairs (q, q+2) has been studied by many investigators, including Glaisher [7], Brun [3], Hardy and Littlewood [9], Sutton [22], Selmer [18], Sexton [19], Lehmer [11], Fröberg [4], Gruenberger and Armerding [8], Weintraub [23], Bohman [1], Shanks and Wrench [20], Brent [2], Nicely [13], and Kutrib and Richstein [10, 25]. Computations similar to those described here have also been carried out by Sebah [17] and initiated by Fry, Nesheiwat, and Szymanski [5].

The present study results from the continuation of a project initiated in 1993, with results to 1e14 described in the author's 1995 paper [13]. A detailed description of the general problem, the computational methods employed, and the incidental discovery of the Pentium FDIV flaw may be found there, with additional details given in [14] and [15]; only a brief summary will be included here.

The prime numbers themselves continue to retain most of their secrets, but still less is known about the twin primes. A matter as fundamental as the infinitude of K_2 remains undecided---the famous "twin-primes conjecture," a topic of discussion even in a recent popular motion picture [21]. Nonetheless, Brun [3] proved in 1919 that in any event the sum of the reciprocals,

(1.1)        B_2 = (1/3 + 1/5) + (1/5 + 1/7) + (1/11 + 1/13) + (1/17 + 1/19) + ... ,

is convergent, in contrast to the known divergence of the sum of the reciprocals of all the primes. The limit of this sum, styled Brun's sum or Brun's constant, is usually denoted by simply B (although Fröberg [4] used S), but herein the author will use B_2. In this instance, as in a number of others noted in this paper, identifiers have been changed from those in [2] and [13], in anticipation of the need for analogous symbols to be used in the study of prime constellations other than the twins.

The twin-primes conjecture is a consequence of a much stronger result proposed by Hardy and Littlewood [9] in 1923,

(1.2)        pi_2(x) ~ L_2(x) = 2c_2*int(1/(ln(t))^2, t, 2, x)

where pi_2(x) represents the count of twin-prime pairs (q, q+2) such that q <= x, and c_2 denotes the "twin-primes constant," computed to 42D by Wrench [24] in 1961,

(1.3)        c_2 = 0.66016 18158 46869 57392 78121 10014 55577 84326 23...

The validity of the conjecture (1.2), sometimes titled the Hardy-Littlewood approximation, is central to the estimation of Brun's constant and the error bounds in this paper. The Hardy-Littlewood approximation is itself a special case of the yet more general "prime k-tuples conjecture," also set forth in [9]. See Riesel [16, pp. 60-83] for an illuminating exposition of these concepts.

Although (1.1) is convergent, the monotonically increasing partial sums approach the limit with agonizing slowness; summing the first thousand million reciprocals is still insufficient to bring us within five percent of the estimated value of the limit. However, assuming the validity of the Hardy-Littlewood approximation (1.2), a first order extrapolation was derived by Fröberg [4] and further studied by Brent [2],

(1.4)        B_2 = S_2(x) + 4c_2/ln(x) + O(1/(sqrt(x)*ln(x))

with an accelerated rate of convergence O(sqrt(x)) faster than (1.1). Here S_2(x) is the partial sum

(1.5)        S_2(x) = sum(1/q + 1/(q+2), q, q <= x) .

of the reciprocals of all the twin-prime pairs (q, q+2) for which q <= x. Note that S_2(x) is written as B(x) in [2] and [13]. The first order extrapolation of S_2(x) to approximate B_2 consists of the first two terms of the right hand side of (1.4); this was indicated as B*(x) in [2] and [13], but we write it here as F_2(x):

(1.6)        F_2(x) = S_2(x) + 4c_2/ln(x) .

The final term in (1.4) is the author's conjectured error or remainder term, inspired by Brent's probabilistic analysis [2]. As discussed by Shanks and Wrench [20, p. 298], no effective second-order extrapolation is known.

2. Computational Technique

These calculations were carried out as part of a more comprehensive project, including in addition enumerations of other prime constellations and prime gaps. Computations began in 1993 and have since proceeded almost without interruption, although several months' work was lost early on due to the Pentium FDIV flaw. The calculations were distributed across a number (varying from a few to more than two dozen) of personal computers, using Intel processors (mostly classic Pentiums), extended DOS and Windows 3.x operating systems, and code written in C. The algorithm employed the classic sieve of Eratosthenes to carry out an exhaustive generation and enumeration of the primes. To guard against errors, all calculations were performed in duplicate on separate systems; in addition, the count pi(x) of primes was maintained and checked periodically against known values. The values obtained for the count pi_2(x) of twin-prime pairs agreed to 1e11 with those of Brent [2], and to 1e14 with those of Kutrib and Richstein [10, 25]. Additional values of pi_2(x) to 8e14 were checked with preliminary values obtained by Fry, Nesheiwat, and Szymanski [5]. In addition to the errors generated by the Pentium FDIV flaw, more than fifty other errors have been detected, most apparently the result of faults---some transient, some permanent---in RAM chips.

The reciprocal sums were calculated in both long double precision (64 bit mantissa, 19S) and ultraprecision (53D), the former in hardware using the Intel numeric coprocessors, the latter in software using a modification of a package due largely to Arjen Lenstra. Results were written to a binary file at intervals of 1e10.

As mentioned previously, additional details regarding the computational technique, and the Pentium FDIV affair, are available in the cited references [13, 14, 15], and also at the author's URL.

3. Computational Results

Table 1 contains a brief summary of the computational results, including the counts pi_2(x) of twin-prime pairs; the values of the discrepancy between pi_2(x) and the Hardy-Littlewood approximation, indicated here by delta_2(x),

(3.1)        delta_2(x) = L_2(x) - pi_2(x) ;

the partial sums S_2(x) of the reciprocals of the twins; and the first order extrapolations F_2(x) of S_2(x) to the limit, according to (1.6), members of a sequence believed to be converging to Brun's constant B_2. Note that the discrepancy delta_2(x) was written in [2] and [13] as r_3(x); Brent also rounded this value to the nearest integer.

Table 1 includes results previously published for powers of ten to 1e14, in addition to new results from the present study for x = (1.5e14)(5e13)(1.6e15). Expanded versions of Table 1 are available at the author's URL, as are extensive tables of the values of pi(x) recorded in this project, at much finer granularity than those commonly available; or see this site for complete results of the computations.

Table 1. Counts of twin-prime pairs and estimates of Brun's constant.
       x        pi_2(x)   delta_2(x)       S_2(x)           F_2(x)
    1e01              2         2.84  0.8761904761905  2.0230090113326
    1e02              8         5.54  1.3309903657191  1.9043996332901
    1e03             35        10.80  1.5180324635596  1.9003053086070
    1e04            205         9.21  1.6168935574322  1.9035981912177
    1e05           1224        24.71  1.6727995848277  1.9021632918562
    1e06           8169        79.03  1.7107769308042  1.9019133533279
    1e07          58980      -226.18  1.7383570439173  1.9021882632233
    1e08         440312        55.79  1.7588156210680  1.9021679379607
    1e09        3424506       802.16  1.7747359576385  1.9021602393210
    1e10       27412679     -1262.47  1.7874785027192  1.9021603562335
    1e11      224376048     -7183.32  1.7979043109551  1.9021605414226
    1e12     1870585220    -25353.18  1.8065924191759  1.9021606304377
    1e13    15834664872    -66566.94  1.8139437606846  1.9021605710802
    1e14   135780321665    -56770.51  1.8202449681303  1.9021605777833
 1.50e14   198466021764    -41424.45  1.8212624999372  1.9021605778478
 2.00e14   259858400254   -286596.19  1.8219692563019  1.9021605806674
 2.50e14   320314774824   -251822.16  1.8225090097373  1.9021605803406
 3.00e14   380041003032   -386165.49  1.8229446574499  1.9021605813179
 3.50e14   439169847802   -377896.90  1.8233092906771  1.9021605812270
 4.00e14   497794845572   -687458.42  1.8236224494488  1.9021605828234
 4.50e14   555983876545   -538028.53  1.8238966154614  1.9021605820874
 5.00e14   613790177314   -495402.94  1.8241402488571  1.9021605819011
 5.50e14   671255245533   -412495.10  1.8243593389029  1.9021605816063
 6.00e14   728412916123   -399030.90  1.8245582810368  1.9021605816028
 6.50e14   785290803354   -219998.15  1.8247403931446  1.9021605810260
 7.00e14   841912734248   -330271.47  1.8249082431040  1.9021605813540
 7.50e14   898298225772   -362599.07  1.8250638545912  1.9021605814440
 8.00e14   954464283498   -207253.20  1.8252088524970  1.9021605810407
 8.50e14  1010426459312   -534774.22  1.8253445623449  1.9021605818329
 9.00e14  1066196920739   -459168.78  1.8254720744001  1.9021605816527
 9.50e14  1121787648852   -422043.55  1.8255923015512  1.9021605815759
 1.00e15  1177209242304   -750443.32  1.8257060132403  1.9021605822498
 1.05e15  1232470313145   -856839.63  1.8258138623433  1.9021605824537
 1.10e15  1287579137984   -732612.87  1.8259164099410  1.9021605822159
 1.15e15  1342543717345   -879045.52  1.8260141417697  1.9021605824791
 1.20e15  1397370335220   -761338.54  1.8261074785719  1.9021605822802
 1.25e15  1452065612798   -774787.26  1.8261967900949  1.9021605822983
 1.30e15  1506635099560   -762644.45  1.8262824008978  1.9021605822837
 1.35e15  1561084329425  -1010551.62  1.8263645987282  1.9021605826562
 1.40e15  1615417411648   -785068.05  1.8264436378766  1.9021605823288
 1.45e15  1669639673260   -931388.83  1.8265197475581  1.9021605825390
 1.50e15  1723754585354   -761213.67  1.8265931311402  1.9021605823084
 1.55e15  1777766583441   -839920.90  1.8266639732589  1.9021605824132
 1.60e15  1831678961614   -851925.37  1.8267324395006  1.9021605824283

4. Brun's Constant and the Error Analysis

The first order extrapolation F_2(1.6e15) is believed to yield the most accurate value known to this point for Brun's constant,

(4.1)        B_2 = 1.90216 05824 +/- 0.00000 00030 .

The error estimate is believed to represent a 99 % confidence level, although this is a subjective opinion of the author for which definitive proof remains lacking. More precisely, the author believes that for at least 99 % of the integers in any sufficiently large set of distinct integers x > 1, Brun's constant B_2 would lie between F_2(x) - E_2(x) and F_2(x) + E_2(x). Here E_2(x) is the error bound function stated in (4.6) below, from which the error bound in (4.1) will be obtained. The error estimate is arrived at as follows.

First, it is assumed that the deviations F_2(x) - B_2 change sign infinitely often; more precisely, given any x_0, however large, there will exist integers x_1 > x_0 and x_2 > x_0 such that F_2(x_1) < B_2 and F_2(x_2) > B_2. We will refer to this as hypothesis [H2]. Note that although [H2] is neither necessary nor sufficient for the Hardy-Littlewood approximation (1.2), it will almost certainly fail (along with (4.1) in its entirety) if Hardy-Littlewood is false. In support of [H2], we note that if our best estimate F_2(1.6e15) is used for B_2, then F_2(x) - B_2 is observed to undergo 872 sign changes over the 160081 data points recorded in the present study, with 677 sign changes occurring beyond 1e15. Brent [2] presented a heuristic argument for a much stronger characterization of the deviations, one which would imply [H2], namely that the values of the scaled deviations

(4.2)        P_2(x) = sqrt(x)*ln(x)*[F_2(x) - B_2]

(indicated in [13] by P(x)) are asymptotically normally distributed with mean o(1) and standard deviation sqrt(8c_2)=2.2981. A statistical analysis of P_2(x), using F_2(1.6e15) for B_2 and the 160081 data points from the present computation, indicates a distribution that is only roughly normal, with mean -0.5294, standard deviation 0.7489, and widely distributed sign changes. As in a similar analysis performed by the author in [13], the standard deviation was significantly less than Brent's prediction; however, this statistic is quite sensitive to the value presumed for B_2.

Given [H2], we then look for the maximum computable value of the function

(4.3)        N_2(x_1, x_2) = sqrt(x_1)*ln(x_1)*|F_2(x_1) - F_2(x_2)|

where x_1 and x_2 are integers such that x_2 > x_1 > 1. N_2 may be considered as a "normalized" measure of the amplitude of the oscillations in F_2(x), or alternately as a measure of the scaled deviation of a "predictor" value F_2(x_1) from a "terminal approximation" F_2(x_2). In contrast to P_2, N_2 has the virtue that it is independent of the uncertain value of B_2. If N_2 has a global maximum, or even an upper bound, then given [H2] this must also represent an upper bound on the absolute value of P_2, thus producing an unconditional error bound; for any specific value |P_2(x_3)| will be exceeded by N_2(x_3, x_4), where x_4 is chosen so that x_4 > x_3 while F_2(x_3) - B_2 and F_2(x_4) - B_2 are of opposite sign ([H2] implies that there is an infinite sequence of such integers x_4 for any given integer x_3). In practice, we are unable to prove that N_2 has a global maximum, although (1.4) would imply that P_2 is bounded. Indeed, it is computationally impractical even to find the absolute maximum of N_2 over all integers in the range under investigation, 1 < x_1 < x_2 <= 1.6e15; this would involve the comparison of more than 1e30 data pairs (F_2(x_1), F_2(x_2)), a computation which might bring an entire Web of teraflop computers to its knees. What has been determined is that the absolute maximum of N_2 over all of the (more than 1e10 such) data pairs recorded in this study is

(4.4)        N_2(1000000, 9000000) = 3.805437 .

However, additional computations reveal an even larger value of N_2 in the vicinity of this point,

(4.5)        N_2(872140, 4265159) = 4.14083 ,

where the values in both (4.4) and (4.5) have been rounded up. Although the value 4.14083 is still not established as the sought after upper bound on N_2 (and thus on |P_2|), the numerical evidence (zero exceptions among more than 1e10 data pairs extending over fifteen orders of magnitude) indicates that if any larger values of N_2 (and thus |P_2|) exist, they must be relatively rare. Our intuitive conclusion is that for the great majority (99 percent or more?) of integers x > 1, it will be true that |P_2(x)| < 4.14083, producing the error bound formula E_2(x),

(4.6)        |F_2(x) - B_2| < E_2(x) = 4.14083/(sqrt(x)*ln(x)) .

In arriving at our error estimate in (4.1), we are assuming that (4.6) does in fact hold specifically for x = 1.6e15, so that

(4.7)        |F_2(1.6e15) - B_2| < E_2(1.6e15) = 4.14083/(sqrt(1.6e15)*ln(1.6e15)) ,


(4.8)        |F_2(1.6e15) - B_2| < 0.00000 00029 57 ,

which, with rounding up for good measure, produces the error estimate given in (4.1).

The error bound formula E_2(x) in (4.6) is similar to the one arrived at (using different reasoning) by Brent [2]. Brent obtains the constant 3.5 in place of 4.14083---partly because he used an estimate for B_2 in place of F_2(x_2), partly because he did not consider it important to perform a higher resolution search for a larger magnitude of the scaled deviation, and partly because of the different interpretation he placed on the resulting error bound; Brent regarded the value 3.5 as representing approximately 1.523 of his standard deviations of P_2(x) (that is, 1.523*sqrt(8c_2)), yielding a probability (based on the presumed normal distribution) of approximately 88 % that B_2 would fall within the resulting error bound. By contrast, the computed standard deviation of P_2(x) here is only 0.7489; even given that this is artificially low due to the use of F_2(1.6e15) as B_2, and taking into account the non-zero mean, this implies that the error bound E_2(x) represents a variation of nearly five standard deviations from the mean, additional support for the opinion of a 99 % (or better) confidence level. On the other hand, the error bound given by the author in [13] was specifically stated as representing one computed standard deviation at 1e14, so it is not especially surprising that the present value for B_2 exceeds the previous estimate by slightly more than two standard deviations (as computed at 1e14).

It should be noted that the scaling factor sqrt(x)*ln(x) was chosen for two reasons---firstly, the persuasiveness of Brent's argument in [2]; secondly, analysis of the data shows that the resulting mean absolute value of P_2(x) remains of the same order of magnitude (typically a few tenths) over most intervals of 1e12 from 0 through 1.6e15, as would be expected with a valid scaling factor. However, other scaling factors could be employed, and the author in fact investigated several alternatives. No other yielded superior results; either the resulting error bound was inferior, or a "near-maximum" value of the corresponding function N_2 (meaning one sufficiently large to provide an error bound of equal or better confidence level) could not be located without greatly escalating the computational effort.

In summary, the error bound E_2(x) set forth here is more conservative than that of Brent [2] and much more conservative than that given in the author's previous work [13]. More than ten thousand million pairs of first order extrapolations have agreed with each other within this error bound, without exception. Nonetheless, it carries no unconditional guarantee; overwhelming numerical evidence is occasionally a subtle trap, as the conjecture Li(x) > pi(x) so memorably illustrated. However, it would appear that a more rigorous error bound must await a deeper understanding of the subtleties in the distribution of the twin primes---perhaps the discovery of a "Riemann zeta function" for the twins.

5. Acknowledgments

The author was greatly aided by Eberhard Mattes' outstanding freeware package emTeX [12], which includes an extended DOS implementation of the LaTeX2e software system for composing and typesetting documents such as this one.


  1. Jan Bohman, "Some computational results regarding the prime numbers below 2,000,000,000," Nordisk Tidskr. Informationsbehandling (BIT) 13 (1973), 242-244; errata, ibid. 14 (1974), 127. MR 48 #217.
  2. Richard P. Brent, "Irregularities in the distribution of primes and twin primes," Math. Comp. 29:129 (1975), 43-56, MR 51 #5522. Corrigendum ibid. 30:133 (1976), 198, MR 53 #302. Addendum reviewed ibid. 30 (1976), 379.
  3. Viggo 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 + ..., où les dénominateurs sont 'nombres premieres jumeaux' est convergente ou finie," (French) Bulletin des sciences mathématiques 43 (1919), 100-104, 124-128.
  4. Carl-Erik Fröberg, "On the sum of inverses of primes and twin primes," Nordisk Tidskr. Informationsbehandling (BIT) 1 (1961), 15-20.
  5. Patrick Fry, Jeffrey Nesheiwat, and Boleslaw Szymanski, "Rensselaer's twin primes computing effort" (ca. 1997-2004); also, related e-mail communications (April-June, 1998). No final results of this work have been found; links to some preliminary and partial results have been collected at http://www.trnicely.net/twins/twins4.html#RTPCE.
  6. Reference removed.
  7. J. W. L. Glaisher, "An enumeration of prime-pairs," Messenger of Mathematics 8 (1878), 28-33.
  8. F. J. Gruenberger and G. Armerding, "Statistics on the first six million primes," UMT 73, Math. Comp. 19 (1965), 503-505.
  9. 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. See also "Collected papers of G. H. Hardy," Clarendon Press, Oxford, 1966, Vol. I, 561-630.
  10. Martin Kutrib and Jörg Richstein, "Primzahlen: Zwillinge aus dem Parallelrechner," (German) Spektrum der Wissenschaft (February, 1996), 26-31.
  11. Derrick Henry Lehmer, "Tables concerning the distribution of primes up to 37 millions," 1957. Copy deposited in the UMT file and reviewed in MTAC 13 (1959), 56-57.
  12. Eberhard Mattes, emTeX package available on the TeX Users Group (TUG) "TeX Live 4" CD (March, 1999), or at (July, 1999) http://www.tex.ac.uk/tex-archive/systems/msdos/.
  13. Thomas R. Nicely, "Enumeration to 1e14 of the twin primes and Brun's constant," Virginia Journal of Science 46:3 (Fall, 1995), 195-204. MR1401560 (97e:11014). Electronic reprint available at http://www.trnicely.net/twins/twins.html.
  14. Thomas R. Nicely, "New maximal prime gaps and first occurrences," Math. Comp. 68:227 (July, 1999), 1311-1315. MR1627813 (99i:11004). Electronic reprint available at http://www.trnicely.net/gaps/gaps.html.
  15. Thomas R. Nicely and Bertil Nyman, "First occurrence of a prime gap of 1000 or greater," unpublished (21 May 1999), available electronically at http://www.trnicely.net/gaps/gaps2.html.
  16. Hans Riesel, "Prime numbers and computer methods for factorization," 2nd ed., Birkhäuser Boston, 1994. 464 pp. + xvi. ISBN 0-8176-3743-5. MR 95h:11142.
  17. Pascal Sebah, unpublished work, e-mail communications (7 June 1999, 25 August 1999). Sebah's results are available at http://numbers.computation.free.fr. Sebah has carried his computations to 1e16, with our results agreeing within rounding error.
  18. Ernst S. Selmer, "A special summation method in the theory of prime numbers and its application to 'Brun's Sum', (Norwegian) " Nordisk Mat. Tidskr. 24 (1942), 74-81. MR 8,316g.
  19. C. R. Sexton, "Counts of twin primes less than 100000," MTAC 8 (1954), 47-49.
  20. Daniel Shanks and John W. Wrench, Jr., "Brun's constant," Math. Comp. 28:125 (January, 1974), 293-299; corrigenda ibid. 28 (1974), 1183. MR 50 #4510.
  21. Barbra Streisand, "The mirror has two faces," TriStar Pictures, 1996. Director, Barbra Streisand. Screen story and screenplay, Richard LaGravenese, based on the motion picture "Le miroir a deux faces" (1959). Cinematographers, Dante Spinotti and Andrzej Bartkowiak. Music, Marvin Hamlisch and Barbra Streisand. Producers, Barbra Streisand and Arnon Milchan. Starring Barbra Streisand, Jeff Bridges, Lauren Bacall, Mimi Rogers, George Segal, Pierce Brosnan, Brenda Vacarro, Austin Pendleton, and Elle Macpherson. See the first date dinner scene.
  22. Charles S. Sutton, "An investigation of the average distribution of twin prime numbers," Journal of Mathematics and Physics 16 (1937), 1-42.
  23. Sol Weintraub, "Four tables concerning the distribution of primes," UMT 38, Math. Comp. 27 (1973), 676-677.
  24. John W. Wrench, Jr., "Evaluation of Artin's constant and the twin-prime constant," Math. Comp. 15 (1961), 396-398. MR 23 #A1619.
  25. Jörg Richstein, "Unendliche Geschwisterliebe unter Zahlen oder Die Suche nach den Primzahlzwillingen," (German) Überblicke Mathematik 1996/1997, 135-143, Überbl. Math., Vieweg, Braunschweig, 1997. MR 99f:11117.


NOTE: This addendum was not part of the original submission for publication.

                             Table 2.
 Counts of primes and twin-prime pairs at primorial values of x.
    x                    x              pi(x)             pi_2(x)
    2#                   2                  1                   0
    3#                   6                  3                   2
    5#                  30                 10                   5
    7#                 210                 46                  15
   11#                2310                343                  70
   13#               30030               3248                 468
   17#              510510              42331                4636
   19#             9699690             646029               57453
   23#           223092870           12283531              896062
   29#          6469693230          300369796            18463713
   31#        200560490130         8028643010           425177757
   37#       7420738134810       259488750744         11997649372
   41#     304250263527210      9414916809095        385088898632
   43#   13082761331670030    362597750396740      13280323588034