Thomas R. Nicely

Current e-mail address

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.

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

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.

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.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.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)) ,

giving

(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.

- 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.
- 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.
- 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.
- Carl-Erik Fröberg, "On the sum of inverses of primes and twin primes," Nordisk Tidskr. Informationsbehandling (BIT) 1 (1961), 15-20.
- 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.
- Reference removed.
- J. W. L. Glaisher, "An enumeration of prime-pairs," Messenger of Mathematics 8 (1878), 28-33.
- F. J. Gruenberger and G. Armerding, "Statistics on the first six million primes," UMT 73, Math. Comp. 19 (1965), 503-505.
- 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.
- Martin Kutrib and Jörg Richstein, "Primzahlen: Zwillinge aus dem Parallelrechner," (German) Spektrum der Wissenschaft (February, 1996), 26-31.
- 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.
- 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/.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- C. R. Sexton, "Counts of twin primes less than 100000," MTAC 8 (1954), 47-49.
- 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.
- 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.
- Charles S. Sutton, "An investigation of the average distribution of twin prime numbers," Journal of Mathematics and Physics 16 (1937), 1-42.
- Sol Weintraub, "Four tables concerning the distribution of primes," UMT 38, Math. Comp. 27 (1973), 676-677.
- John W. Wrench, Jr., "Evaluation of Artin's constant and the twin-prime constant," Math. Comp. 15 (1961), 396-398. MR 23 #A1619.
- 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.

- After further consideration, the author is now (2009) of the opinion that the error estimates appearing in this paper are better characterized as representing only "at least one standard deviation," that is, a confidence interval of at least 68.27 %.
- The latest updated count of the twin primes, along with corresponding estimates for Brun's constant and a definitive error bound, may be found at http://www.trnicely.net/counts.html.
- More extensive versions of Table 1 are available on this site for the intervals 0 to 1e16 and 1e16 to 2e16.
- Also available are detailed tables of pi(x) and associated functions for the intervals 0 to 1e16 and 1e16 to 2e16.
- Complete counts and reciprocal sums of the prime constellations from Nicely's computations (1993-2009), including over two million data points from 0 to 2e16, are now available (these are very large files, over 60MB each even in zipped form).
- The reciprocal sums S_2(x) listed in Table 1 are believed to be correct, and correctly rounded, to all digits shown. The calculation of these quantities was carried out in long double precision (extended precision, 64-bit mantissa, 19 significant digits), using the Intel FPU hardware (numeric coprocessor). However, these calculations were also carried out to 53 decimal places in software, using the author's modification of the ZBIGINT ultraprecision integer package developed, and graciously placed in the public domain, by Arjen K. Lenstra, Mark Riordan, and Marc Ringuette (1988-1991). The original ZBIGINT package is available at http://www.funet.fi/pub/crypt/cryptography/rpem/rpem/ (thanks to Charles Doty for this pointer). The predecessor of this package was Lenstra's LIP (long integer package). Its descendant is NTL (Number Theory Library), maintained by Victor Shoup. However, these packages have to a large extent been rendered obsolete by GMP (GNU MP), the GNU multiple precision library, and MPFR, a C library for multiple-precision floating-point computations with correct rounding, reliable precision control, and compatibility with the ANSI/IEEE 754-1985 standard.
*This note is in response to an inquiry from Dr. Marek Wolf*. Efforts to use weighted and/or unweighted data averaging or smoothing, or linear regression techniques, in an attempt to obtain a more accurate value of Brun's constant, have not met with success. The results thus obtained, at values of x < 1.6e15, were in some cases slightly better, but also in some cases slightly worse, than the values of F_2(x). "Better" or "worse" here refers to whether or not the resulting estimate was closer to the presumably superior estimates of B_2 obtained from more extensive calculations, i.e., F_2(x) for a much larger value of x. However, I have not mounted a full-scale effort using weighted or non-linear regression; I invite third parties to carry out such an effort, using the extensive data tables I have provided. It may well be possible to produce a more rapidly converging sequence of estimates by some such method; this would amount to a successful second order extrapolation. I recommend as critical test values for this purpose the estimates for Brun's constant at x = 1e11 and x = 1e13; can your proposed extrapolation produce estimates substantially superior to F_2(1e11) and F_2(1e13)?- At the request of Professor Kenneth Hicks, Department of Physics, Ohio University, Athens, Ohio, the following table has been prepared.

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 ==================================================================

- Top of Page
- "Enumeration to 1e14 of the twin primes and Brun's constant" (paper, 1995)
- "A new error analysis for Brun's constant" (paper, 2001)
- The latest count of the twin primes and Brun's constant
- Enumeration of the twin-prime pairs to 1e16 (table)
- Enumeration of the twin-prime pairs from 1e16 to 2e16 (table)
- Complete counts and reciprocal sums of the prime constellations from Nicely's computations (tables)
- Home Page