A new error analysis for Brun's constant

A new error analysis for Brun's constant
Thomas R. Nicely

http://www.trnicely.net
Current e-mail address

Document History

Last modified....................0700 GMT 06 October 2011.
Mathematical Reviews.............MR1853722 (2003d:11184).
Journal citation.................Virginia Journal of Science 52:1 (Spring, 2001) 45-55.
Printing date....................April 2001.
Accepted for publication.........26 January 2001.
Accepted for consideration.......09 November 2000.
Original submission..............06 November 2000 (Virginia Journal of Science).

Except for the addendum (not part of the submission for publication), the content of this document is essentially that of the journal article as published; it may contain minor revisions and corrections (such as updated URLs), and differ in format and detail.

Abstract

Enumeration of the twin primes, and the sum of their reciprocals, is extended to 3e15, yielding the count pi_2(3e15) = 3310517800844. A more accurate estimate is obtained for Brun's constant, B_2 = 1.90216 05823 +/- 0.00000 00008. Error analysis is presented to support the contention that this estimate produces a 95 % confidence interval for B_2. In addition, published values of the count pi(x) of primes, obtained previously by indirect means, are verified by direct count to x = 3e15.

Mathematics Subject Classification 2000 (MSC2000)

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

Key Words and Phrases

Twin primes, Brun's constant, prime numbers.

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 (1878), Brun (1919), Hardy and Littlewood (1922/23), Sutton (1937), Selmer (1942), Sexton (1954), Lehmer (1957), Fröberg (1961), Gruenberger and Armerding (1965), Weintraub (1973), Bohman (1973), Shanks and Wrench (1974), Brent (1975), and Nicely (1995).

The present study results from the continuation of a project initiated in 1993, with results to 1e14 previously published (Nicely, 1995). 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 (Nicely, 1999); 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." Nonetheless, Brun (1919) proved that in any event the sum of the reciprocals,

(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 (Brun actually omitted the first term in parentheses, which of course does not affect the convergence). The limit of this sum, styled Brun's sum or Brun's constant, is often denoted as simply B, but henceforth 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 (Brent, 1975) and (Nicely, 1995), 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, an asymptotic relationship conjectured by Hardy and Littlewood (1922/23, pp. 42-44):

(2) pi_2(x) ~ L_2(x) = 2c_2*integral[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 (1961),

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

The validity of the conjecture (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 consequence of the yet more general "prime k-tuples conjecture," also set forth in their 1922/23 work. See Riesel (1994, pp. 60-83) for an illuminating exposition of these concepts.

Although (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 (2), a first-order extrapolation was derived by Fröberg (1961) and further studied by Brent (1975),

(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). Here S_2(x) is the partial sum

(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 (Brent, 1975) and (Nicely, 1995). The first-order extrapolation of S_2(x) to approximate B_2 consists of the first two terms of the right hand side of (4); this was indicated as B*(x) in (Brent, 1975) and (Nicely, 1995), but we write it here as F_2(x):

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

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

Computational Technique

These calculations were carried out as part of a more comprehensive project, including in addition the tabulation of prime gaps and other prime constellations. 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 asynchronously across several (varying from a few to more than two dozen) personal computers, using Intel processors (mostly classic Pentiums), extended DOS and Windows 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, such as those published by Riesel (1994, pp. 380-383). The values obtained for the count pi_2(x) of twin-prime pairs agreed to 1e11 with those of Brent (1975), and to 1e14 with those of Kutrib and Richstein (1996, 1997). Excluding software bugs and the Pentium FDIV flaw, more than forty instances of machine errors were detected and corrected, most apparently the result of transient bit errors in memory (DRAM) chips. One of these instances contained at least 364 individual errors.

As mentioned previously, additional details regarding the computational technique, and the Pentium FDIV affair, are available in (Nicely, 1995, 1999), and also at the author's URL.

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, denoted here by delta_2(x), between pi_2(x) and the Hardy-Littlewood approximation:

(7) 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 (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 (Brent, 1975) and (Nicely, 1995) as r_3(x); Brent also rounded this value to the nearest integer.

   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)
=============================================================================
 1.0e+01              2         2.84  0.8761904761904761905  2.0230090113326
 1.0e+02              8         5.54  1.3309903657190867570  1.9043996332901
 1.0e+03             35        10.80  1.5180324635595909885  1.9003053086070
 1.0e+04            205         9.21  1.6168935574322006462  1.9035981912177
 1.0e+05           1224        24.71  1.6727995848277415480  1.9021632918562
 1.0e+06           8169        79.03  1.7107769308042211063  1.9019133533279
 1.0e+07          58980      -226.18  1.7383570439172709388  1.9021882632233
 1.0e+08         440312        55.79  1.7588156210679749679  1.9021679379607
 1.0e+09        3424506       802.16  1.7747359576385368007  1.9021602393210
 1.0e+10       27412679     -1262.47  1.7874785027192415475  1.9021603562335
 1.0e+11      224376048     -7183.32  1.7979043109551191615  1.9021605414226
 1.0e+12     1870585220    -25353.18  1.8065924191758825917  1.9021606304377
 1.0e+13    15834664872    -66566.94  1.8139437606846070596  1.9021605710802
 1.0e+14   135780321665    -56770.51  1.8202449681302705289  1.9021605777833
 2.0e+14   259858400254   -286596.19  1.8219692563019236634  1.9021605806674
 3.0e+14   380041003032   -386165.49  1.8229446574498899187  1.9021605813179
 4.0e+14   497794845572   -687458.42  1.8236224494488219106  1.9021605828234
 5.0e+14   613790177314   -495402.94  1.8241402488570614635  1.9021605819011
 6.0e+14   728412916123   -399030.90  1.8245582810368460212  1.9021605816028
 7.0e+14   841912734248   -330271.47  1.8249082431039834264  1.9021605813540
 8.0e+14   954464283498   -207253.20  1.8252088524969516994  1.9021605810407
 9.0e+14  1066196920739   -459168.78  1.8254720744000806297  1.9021605816527
 1.0e+15  1177209242304   -750443.32  1.8257060132402797152  1.9021605822498
 1.1e+15  1287579137984   -732612.87  1.8259164099409972759  1.9021605822159
 1.2e+15  1397370335220   -761338.54  1.8261074785718993129  1.9021605822802
 1.3e+15  1506635099560   -762644.45  1.8262824008978027694  1.9021605822837
 1.4e+15  1615417411648   -785068.05  1.8264436378766369280  1.9021605823288
 1.5e+15  1723754585354   -761213.67  1.8265931311402050729  1.9021605823084
 1.6e+15  1831678961614   -851925.37  1.8267324395006005931  1.9021605824283
 1.7e+15  1939218595600  -1129122.83  1.8268628327687977085  1.9021605827604
 1.8e+15  2046397121805   -678331.73  1.8269853577548725890  1.9021605822393
 1.9e+15  2153237307407   -562823.58  1.8271008903959923363  1.9021605821153
 2.0e+15  2259758303674   -612652.24  1.8272101680098151140  1.9021605821628
 2.1e+15  2365977242191   -653062.89  1.8273138179643056714  1.9021605822014
 2.2e+15  2471909670028   -643465.53  1.8274123785364204712  1.9021605821937
 2.3e+15  2577569863563   -750111.35  1.8275063150448871463  1.9021605822851
 2.4e+15  2682970233099   -552427.29  1.8275960317894826243  1.9021605821145
 2.5e+15  2788122612616   -168258.89  1.8276818830618905359  1.9021605818032
 2.6e+15  2893038573759   -430246.96  1.8277641812367275962  1.9021605820124
 2.7e+15  2997726948096   -292107.29  1.8278432012390461693  1.9021605819106
 2.8e+15  3102197972961   -876051.32  1.8279191890118998763  1.9021605823359
 2.9e+15  3206458423771   -521046.38  1.8279923621701145073  1.9021605820865
 3.0e+15  3310517800844   -897422.15  1.8280629180352850193  1.9021605823404
=============================================================================

Table 1 includes results previously published for powers of ten to 1e14, in addition to new results from the present study at additional increments of 1e14, ending with the results for the present upper bound of computation, x_0 = 3e15. Updated and more extensive versions of Table 1 are available at the author's URL, as are extensive tables of the associated values obtained for the count of primes, pi(x). Indeed, the direct enumeration of the primes has been extended to a new upper bound by this project, culminating in the value pi(3e15) = 86688602810119. This result, as well as a number of other previously published values (Riesel, 1994, pp. 380-382) which were known only through indirect calculations, is now confirmed by direct count in the present study.

Brun's Constant and the Error Analysis

The first-order extrapolation F_2(x_0) = F_2(3e15) is believed to yield the most accurate value known to date for Brun's constant,

(8) B_2 = 1.90216 05823 +/- 0.00000 00008 .

The error estimate is believed to define a 95 % confidence interval for the value of B_2. I have no rigorous proof of this assertion regarding the error estimate; rather it is an inference from the analysis (presented below) of the available numerical data. The notion of a "95 % confidence interval" is to be interpreted as follows. Based on the available numerical data, the author believes that whenever the technique used for this error analysis is applied to a sufficiently numerous sample of distinct integers x > 1, Brun's constant B_2 will lie between F_2(x) - E_2(x) and F_2(x) + E_2(x) for at least 95 % of the integers in the sample. Here E_2(x) is the error bound function stated in (11) below; the error estimate given in (8) is a special case of this error bound function, namely E_2(x_0). More precisely, given any set Z_1 of distinct integers x > 1, there will always exist a superset Z_2 of distinct integers x > 1, Z_1 \subseteq Z_2, such that F_2(x) - E_2(x) <= B_2 <= F_2(x) + E_2(x) for at least 95 % of the integers in Z_2.

The algorithm for obtaining and validating this error bound function will now be explained. Discussion and justification of certain details of the procedure will be deferred until a later point in this paper.

(A) A set S of sample test points is chosen from the available numerical data; this set should be a reasonably large subset of all the available data points, avoiding any known bias in the associated values of S_2 or F_2. Indeed, S might be chosen as the entire set T of all recorded data points, up to and including the current upper bound x_0 = 3e15 of computation; there are 300081 points in T, consisting of the lattice (1e10)(1e10)(x_0) together with the "decade values" x = k*10^n, (k=1...9, n=1...9). However, the calculations to be carried out in the error analysis then become excessive. We choose instead for S the lattice (1e12)(1e12)(x_0), consisting of 3000 equally spaced data points, extending to the current upper limit of computation, the increment being one (U. S.) trillion.

(B) For each x in S, we obtain an error bound on F_2(x), presumably representing a 95 % confidence interval, by determining the value of a parameter K_95(x) such that, for at least 95 % of the points in the set U = {t: t in T, t <= x/2},

(9) |F_2(x) - F_2(t)| < K_95(x) /(sqrt(t)*ln(t)) .

Here the form of the "scaling factor" in the denominator is inferred from the remainder term conjectured in (4). The data points t > x/2 are excluded from U to minimize any artificial reduction in the error estimate resulting from the implicit bias of F_2(t) toward F_2(x) as t approaches x.

(C) We now reason as follows. Since for each x in S, at least 95 % of the (relevant) preceding extrapolations F_2(t), t in U, agree with F_2(x) within the bound in (9), we assume that this property will remain valid for arbitrarily large values of x as well. We now interchange x and t in (9), as well as the order of the resulting terms on the left hand side, and take the limit as t approaches plus infinity.

(10) |F_2(x) - lim(F_2(t), t, +infinity)| <= lim(K_95(t), t, +infinity)/(sqrt(x)*ln(x)) .

The numerical evidence indicates that the positive function K_95(x) is either roughly constant, or exhibits an overall decreasing trend masked by small scale variations (see Table 2). Thus we can obtain an approximate upper bound on the error by using K_95(x) in place of the (unknown) limit of K_95(t) in (10). This produces the desired error bound function E_2(x):

(11) |F_2(x) - B_2| <= E_2(x) = K_95(x)/(sqrt(x)*ln(x)) .

Determination of the error bound at any specific x then becomes a matter of calculating K_95(x) and substituting into (11).

Analysis of the data yields the value K_95(x_0) = 1.380. Substitution into (11) then gives

(12) E_2(x_0) = 1.380/(sqrt(x_0)*ln(x_0)) = 0.00000 00007 06989 .

Rounding up produces the error estimate stated in (8).

Validation of the Error Analysis

Since the ad hoc error analysis algorithm described and employed clearly lacks a rigorous analytical foundation, additional examination of the empirical evidence was undertaken, in an effort to find supporting evidence, or lack thereof.

The validation process consisted of comparing the confidence intervals obtained for B_2 at each x in the "lower half" S' = {x: x in S, x <= x_0/2} of S (the values near x_0 being excluded for reasons similar to those given for set U) with the (presumably) best value obtained at x_0. Simply put, the issue is this: what percentage of the confidence intervals obtained for each x in S' actually contain the best known point estimate for B_2, given in (8) (and to greater precision, if not accuracy, in the last entry of Table 1)? For example, applying our error analysis technique to the data for x <= 1e14 yields K_95(1e14) = 1.758, and substitution into (11) then produces the confidence interval B_2 = 1.90216 05777 83 +/- 0.00000 00054 53. Since our best estimate for B_2 lies within this interval, we consider the error estimate algorithm to be a success at x = 1e14. On the other hand, applying the algorithm to the data for x <= x_1 = 8.13e14, we obtain K_95(x_1) = 1.306, with the resulting confidence interval B_2 = 1.90216 05809 53 +/- 0.00000 00013 34, which constitutes a failure.

A survey of all the points x in S' reveals that 96.93 % (1454 of 1500) produce confidence intervals containing our current best point estimate for B_2. These calculations, briefly summarized in Table 2, also show the trends in the values of K_95(x), E_2(x), and the cumulative percentage of successful (in the sense described above) error estimates generated by the algorithm. The available data thus indicates that our algorithm has been successful (actually performing beyond expectation) in producing valid 95 % confidence intervals for the estimates of B_2. Therefore we anticipate that the error bounds thus obtained for larger values of x, including our current upper bound of computation x_0 = 3e15, will also yield valid 95 % confidence intervals for the value of B_2.

 TABLE 2. Performance data for the error analysis algorithm.
===============================================================
 x/1e12          K_95(x)       E_2(x)*1e10        Success %
===============================================================
     1            2.074           750.61           100.00
    10            2.218           234.32            80.00
   100            1.758            54.53            94.00
   200            1.602            34.40            83.50
   300            1.582            27.40            89.00
   400            2.047            30.44            91.75
   500            1.688            22.30            93.40
   600            1.564            18.76            94.50
   700            1.451            16.04            94.86
   800            1.320            13.60            95.25
   900            1.487            14.39            94.89
  1000            1.658            15.18            95.40
  1100            1.623            14.13            95.82
  1200            1.626            13.52            96.17
  1300            1.608            12.82            96.46
  1400            1.606            12.31            96.71
  1500            1.584            11.70            96.93
  1600            1.612            11.51            97.12
  1700            1.747            12.08            97.29
  1800            1.511            10.14            97.44
  1900            1.455             9.49            97.58
  2000            1.454             9.23            97.70
  2100            1.450             8.97            97.81
  2200            1.434             8.65            97.91
  2300            1.449             8.54            98.00
  2400            1.381             7.96            98.08
  2500            1.269             7.16            98.16
  2600            1.321             7.30            98.23
  2700            1.277             6.92            98.30
  2800            1.399             7.43            98.36
  2900            1.307             6.82            98.41
  3000            1.380             7.07            98.47
===============================================================

Critique and Further Remarks

The results of additional analysis of the data, conducted in order to address various weaknesses of the error analysis algorithm described, are now summarized.

Further reduction of the "cutoff" fraction for the selection of sample points in sets U and S' (for example, restricting these sets to the smallest quarter, rather than the smaller half of the eligible values) had no significant effect on the results. Of course, if the restriction is relaxed or eliminated, the effect is to artificially inflate the success percentage of the algorithm. This may be observed in the entries of the last column of Table 2, for values of x > 1.5e15.
Increasing the density of the sample sets S and S' in T (for example, reducing the increment to 1e11 rather than 1e12) had no significant effect on the results.
Replacing the presumed best estimate F_2(x_0) for B_2 by another value within the specified confidence interval (8) (both endpoint values were tested) had no significant effect on the conclusions.
Replacing Brent's (1975) scaling factor sqrt(x)*ln(x) (corresponding to the denominator of the remainder term in (4)) with other plausible possibilities had no significant effect on the results. Among the candidates checked were sqrt(x)*ln(x)*ln(ln(x)), sqrt(x)*(ln(x))^2, sqrt(x)*ln(x)*(ln(ln(x)))^2, sqrt(x)*ln x* ln(ln(x))*ln(ln(ln(x))), and sqrt(x). Results produced by each of these scaling factors are summarized in Table 3; note that the values for K_95 and the error are calculated at x = x_0 = 3e15, while the success percentages are evaluated at x = x_0/2 = 1.5e15, as in our principal error analysis; furthermore, the values of the error are in units of 1e(-10). The available numerical data is seen to be insufficient to either confirm or reject the error term conjectured by the author in (4), or any of the alternatives. On the other hand, since the use of these alternatives had little impact on the final results of the error analysis, the validity of the algorithm appears to be relatively insensitive to the precise nature of the scaling factor (remainder term). Let it be noted that one could make a case, based on the results in Table 3, for a more aggressive error estimate of 0.00000 00006 58 in (8); the author prefers the more conservative value previously stated.

     TABLE 3. Impact of various scaling factors on the error analysis.
============================================================================
 Scaling factor                          K_95         Error       Success %
============================================================================
 sqrt(x)*ln(x)                           1.380         7.07         96.93

 sqrt(x)*ln(x)*ln(ln(x))                 4.809         6.89         96.27

 sqrt(x)*[ln(x)]^2                       45.40         6.53         94.80

 sqrt(x)*ln(x)*[(ln(ln(x))]^2            16.80         6.74         95.67

 sqrt(x)*ln(x)*ln(ln(x))*ln(ln(ln(x)))   6.016         6.77         95.87

 sqrt(x)                                 0.043         7.83         99.00
============================================================================

Other analysis techniques were investigated as well, but none was found superior to the one described. Efforts to use weighted 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. Harmonic analysis and fast Fourier transforms have been suggested by various colleagues as promising techniques for analysis of the data, but I have not pursued this avenue. I will attempt to post enough of the raw data at my URL so that other investigators may experiment with their own techniques; perhaps some other method will indeed be more successful than my own in producing a more accurate extrapolation, or a superior error bound.
The error bound formula E_2(x) in (11) is a generalization of that obtained by Brent (1975). As a consequence of a quite different line of reasoning, Brent arrived at the constant 3.5 in place of K_95(x), and believed this to produce an 88 % confidence interval for his estimate for B_2. It now appears that Brent's error estimate was quite conservative. On the other hand, the error bound obtained in (Nicely, 1995) was specifically designed to represent one computed standard deviation at 1e14, and the present estimate for B_2 differs from that value by more than two of those standard deviations. As pointed out above, the present technique, when applied to the portion of the data for x <= 1e14, produces a 95 % confidence interval containing the current best estimate for B_2 (and even containing the entire current best confidence interval); since a 95 % confidence interval corresponds to about +/- 1.96 standard deviations (for a normal distribution), that result implies sigma(1e14) = 0.00000 00028, a more conservative value than the estimate of 0.00000 00021 arrived at (using a different approach) in (Nicely, 1995).
As the upper bound x_2 of computation for pi_2(x) and S_2(x) is extended, corresponding error estimates can be obtained by analyzing the new totality of data to determine K_95(x_2), according to part (B) of the error analysis algorithm, then substituting into (11). Note that there is no need to recompute K_95(x) for any x other than the value for which a new error bound is desired; computation of K_95(x) over the entire sample set was carried out only to explain and validate the algorithm. Indeed, based on the variation exhibited by K_95(x) in Table 2, a rough error estimate could be obtained by simply using K_95(x_2) = K_95(x_0) = 1.380; or if a quite conservative value is desired, use K_95(x_2) = 2.
Finally, it must be emphasized that both the value of B_2 and the associated error estimate obtained in this paper are entirely dependent on the validity of the Hardy-Littlewood approximation (2). All the numerical evidence to date strongly supports this conjecture, but one must maintain some informed skepticism; after all, the numerical evidence to the current level of computations also supports the famous conjecture that Li(x) > pi(x), eventually disproved by Littlewood (1914) himself. Absent a major theoretical breakthrough, it will be difficult indeed to improve significantly on either the estimate or error bound herein presented for Brun's constant. As Shanks and Wrench (1974, p. 299) noted, the calculation of B_2 to eight or nine decimals is (was) extremely difficult---or at least computationally intensive---and twenty decimals of precision remains as remote now as it was then. Equation (11) indicates that computations may have to be extended to 1e17 just to settle the tenth decimal place, and twenty decimals would require calculations out to perhaps 1e36---a figure far exceeding the total number of machine cycles available in the cumulative projected lifetimes of all the CPUs currently on our planet.

Literature Cited

Bohman, Jan. 1973. Some computational results regarding the prime numbers below 2,000,000,000. Nordisk Tidskr. Informationsbehandling (BIT) 13:242-244. Errata (1974) ibid. 14:127. MR 48 #217.
Brent, Richard P. 1975. Irregularities in the distribution of primes and twin primes. Math. Comp. 29(129):43-56, MR 51 #5522. Corrigendum (1976) ibid. 30(133):198, MR 53 #302. Addendum (1976) reviewed ibid. 30:379.
Brun, Viggo. 1919. 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. Bulletin des sciences mathématiques 43:100-104, 124-128.
Fröberg, Carl-Erik. On the sum of inverses of primes and twin primes. Nordisk Tidskr. Informationsbehandling (BIT) 1:15-20.
Glaisher, J. W. L. 1878. An enumeration of prime-pairs. Messenger of Mathematics 8:28-33.
Gruenberger, F. J., and G. Armerding. 1965. Statistics on the first six million primes. UMT 73. Math. Comp. 19:503-505.
Hardy, G. H., and J. E. Littlewood. 1922/1923. Some problems of 'Partitio Numerorum'; III: On the expression of a number as a sum of primes. Acta Math. 44:1-70. See also Collected papers of G. H. Hardy (1966), Clarendon Press, Oxford I:561-630.
Kutrib, Martin, and Jörg Richstein. 1996. Primzahlen: Zwillinge aus dem Parallelrechner. Spektrum der Wissenschaft (February, 1996), 26-31.
Lehmer, Derrick Henry. 1957. Tables concerning the distribution of primes up to 37 millions. UMT 3. Reviewed in MTAC (1959) 13:56-57.
Littlewood, J. E. 1914. Sur la distribution des nombres premiers. Comptes Rendus Acad. Sci. Paris 158:1869-1872.
Nicely, Thomas R. 1995. Enumeration to 1e14 of the twin primes and Brun's constant. Virginia Journal of Science 46(3):195-204 (Fall, 1995). MR1401560 (97e:11014). Electronic reprint available at http://www.trnicely.net/twins/twins.html.
Nicely, Thomas R. 1999. New maximal prime gaps and first occurrences. Math. Comp. 68(227):1311-1315 (July, 1999). MR1627813 (99i:11004). Electronic reprint available at http://www.trnicely.net/gaps/gaps.html.
Richstein, Jörg. 1997. 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.
Riesel, Hans. 1994. Prime numbers and computer methods for factorization. 2nd edition. Birkhäuser Boston. 464 pp. + xvi. ISBN 0-8176-3743-5. MR 95h:11142.
Selmer, Ernst S. 1942. A special summation method in the theory of prime numbers and its application to 'Brun's Sum.' (Norwegian) Nordisk Mat. Tidskr. 24:74-81. MR 8,316g.
Sexton, C. R. Counts of twin primes less than 100000. MTAC 8:47-49.
Shanks, Daniel and Wrench, John W., Jr. 1974 (January). Brun's constant. Math. Comp. 28(125):293-299; corrigenda (1974) ibid. 28:1183. MR 50 #4510.
Sutton, Charles S. 1937. An investigation of the average distribution of twin prime numbers. Journal of Mathematics and Physics 16:1-42.
Weintraub, Sol. 1973. Four tables concerning the distribution of primes. UMT 38. Math. Comp. 27:676-677.
Wrench, John W., Jr. 1961. Evaluation of Artin's constant and the twin-prime constant. Math. Comp. 15:396-398. MR 23 #A1619.

Unpublished Addendum

NOTE: This addendum is not part of the published journal article.

As of 2009, the author believes the error analysis in the paper "Enumeration to 1.6e15 of the twin primes and Brun's constant" (1999) to be more accurate than the one presented in this paper. It obtains a considerably more conservative result, namely, that the error bound at any specific value of x is given by
```
    |B*(x) - B| < E_2(x) = 4.14083/(sqrt(x)*ln(x))
```
See the referenced paper for details. By this formula, the result of the calculations up to 3e15 becomes
```
    B = 1.90216 05823 +/- 0.00000 00022
```
Furthermore, although the cited paper represents this error bound as having a 99 % confidence level, the author is of the opinion (2009) that this error bound is best characterized as representing only "at least one standard deviation," i.e., a confidence level of at least 68.27 %.
R. F. Arenstorf of Vanderbilt University published a proposed proof (26 May 2004) of the twin-primes conjecture. Arenstorf used analytic continuation, the complex Tauberian theorem of Wiener and Ikehara, and other advanced techniques of classical analytic number theory in an attempt to prove a related conjecture of Hardy and Littlewood which has the twin-primes conjecture as a corollary. Following the discovery by Gerald Tenenbaum of a significant error in Lemma 8, the paper was withdrawn (8 June 2004) for revision. My thanks to Richard P. Brent of the Oxford University Computing Laboratory and the Australian National University for directing me to Arenstorf's paper.
See also my paper "New evidence for the infinitude of some prime constellations" (20 July 2004) for some additional implications of Arenstorf's work.
The latest updated count of the twin primes, and the corresponding estimates for Brun's constant, may be found at http://www.trnicely.net/counts.html.
An extended version of Table 1, for values of x to 1e16, is available at http://www.trnicely.net/twins/t2_0000.htm. Also available is a table of values for x between 1e16 and 2e16.
Corresponding tables of the values of pi(x) and its associated functions, determined by direct count as part of this project, are available for the interval 0 to 1e16 and for the interval 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).
A similar project has been initiated by Pascal Sebah, with results available at http://numbers.computation.free.fr. Sebah has carried his computations to 1e16, with our results agreeing within rounding error.

Another project similar to this one was initiated by Patrick Fry, Jeffrey Nesheiwat, and Boleslaw K. Szymanski, titled "Rensselaer's twin primes computing effort" (ca. 1997-2004). Fragmentary results of this work have been found at the following locations (thanks to Guy H. Bearman for these pointers):
- " Experiences with distributed computation of twin primes distribution," P. Fry, J. Nesheiwat, and B. Szymanski, Parallel and Distributed Computing Practices, Volume 2, Number 3, pp. 299-316 (November, 1999).
- " Computing twin primes and Brun's constant: a distributed approach," P. Fry, J. Nesheiwat, and B. Szymanski, Proc. 7th IEEE Int. Symposium on High Performance Distributed Computing, Chicago IL (July, 1998), IEEE Computer Society Press, Los Alamitos CA (1998), pp. 42-49.
- See also http://web.archive.org/web/ 19991004124128/http://www.cs.rpi.edu/research/twinp/.
- See also http://web.archive.org/web/*/ http://www.cs.rpi.edu/research/twinp/.
Note that the documents at these sites appear to be directed almost exclusively at the computer science aspects of the computations. Mathematical results are generally sparse, incomplete, ambiguous, unreliable, and occasionally erroneous. Many of the associated links produce only error messages (404s or equivalents). Caveat emptor !
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.
The count of machine error instances has now reached 58, the latest consisting of numerous errors observed on a Dell desktop (Dimension 3100, 3GHz Pentium 4) circa 16 December 2008. The errors appear to have been caused by a defective memory module; after it was replaced the machine repeated the interval correctly, and has exhibited no further errors since.
The twin-primes constant c_2 in (3) (as well as a number of other constants of interest in number theory) has been recalculated (1999-2000) by Gerhard Niklasch and Pieter Moree to 1002 decimal places. The relevant Web document is no longer accessible.
I have had several inquiries regarding the search for sequences of consecutive twin-prime pairs (including no other intervening unpaired primes), sometimes referred to as "twin-prime constellations" or "twin-prime clusters." This has not thus far been a subject of my research. Individuals who have communicated such results include the following:
- Nigel B. Backhouse (Liverpool, UK).
- Denis DeVries.
- Tony Forbes.
- Lévai Gábor.
- Robert Hein. His results were received as unpublished work, via telefax (3 March 1995).
- Jim Morton.
- Randall L. Rathbun.
The first occurrence of a sequence of eleven consecutive twin-prime pairs has been discovered by Lévai Gábor (01 October 2011), in the interval [789795449254776509, 789795449254776871]. No instance of twelve consecutive twin-prime pairs is presently known.
In addition, Professor Ken Hicks (Ohio University) has discovered (ca. March 2000) an instance of eight twin-prime pairs in arithmetic progression, consisting of the set {(a + nd, a + nd + 2), n=0..7} with a=11543661551 and d=6469693230. Professor Hicks believes this to be the longest known sequence of twin-prime pairs in arithmetic progression.