Brun's constant (WP 34S)

[Gerson W. Barbosa]

In 1919 the Norwegian mathematician Viggo Brun proved that the sum of the reciprocals of the twin primes (pair of prime numbers which differ by 2) converges to a finite value now known as Brun's constant (B₂). Unlike some constants related to divergent sums of reciprocals of integers, like the Euler-Mascheroni constant (harmonic series) and Mertens constant (prime numbers), which are respectively known to millions and thousands of digits, the Brun's constant is known to 9 or 10 digits only. The sum converges very slowly, so slowly that the sum will not reach the value 1.9 until all the reciprocals of the twin prime pairs up to 10⁵³⁰ are summed up. Thus, an indirect method is used, assuming the Twin Prime Conjecture is true:

     B*2 = B2(p) + 4*C2/log(p)

where

     B2(p) = sum of the reciprocals of the twin prime pairs up to p

     C2 = Twin primes constant (0.66016181584686957..)

     B*2 = Approximation to Brun's constant

B^*₂ has been calculated for p up to 10¹⁵ and 10¹⁶ by Thomas R. Nicely (1999) and Pascal Sebah (2012), respectively:

http://www.trnicely.net/twins/twins2.html

http://numbers.computation.free.fr/Constants/Primes/twin.pdf

By the way, the famous Pentium bug was discovered by Dr. Nicely in 1994 when calculating B₂ for p up to 10¹⁴.

Let's now compute a the constant to a few digits on the WP-34S:

001 LBL A		018 x<>y   		

002 STO 03		019 RCL L		

003 2			020 x>? 03			

004 +/-			021 SKIP 007

005 STO 02		022 ||

006 0			023 RCL L

007 STO 01		024 x<>y		

008 +			025 1/x	

009 RCL L		026 STO+ 01		

010 DEC X		027 x<>y

011 NEXTP		028 BACK 018		

012 ENTER^		029 RCL 00 		

013 ENTER^		030 RCL 03		

014 NEXTP		031 LN

015 -			032 /

016 x<>? 02		033 RCL+ 01

017 BACK 008            034 END

    2.64064726339 STO 00       ; 4*C2 

   10 A  -->  2.02300901133    ; B21 

 RCL 01  -->  0.87619047619    ; (1/3 + 1/5) + (1/5 + 1/7)

  100 A  -->  1.90439963329    ; B22  

 RCL 01  -->  1.33099036572    ; (1/3 + 1/5) + (1/5 + 1/7) + (1/11 + 1/13) + ... + (1/83 + 1/89)

EEX 3 A  -->  1.90030530861    ; B23 

 RCL 01  -->  1.51803246356    ; (1/3 + 1/5) + (1/5 + 1/7) + (1/11 + 1/13) + ... + (1/991 + 1/997) 

EEX 4 A  -->  1.90359819122    ; B24  

EEX 5 A  -->  1.90216329186    ; B25 

EEX 6 A  -->  1.90191335333    ; B26 

EEX 7 A  -->  1.90218826322    ; B27  (about 15 minutes on the WP-34S emulator)

The most accurate value of B₂ to date is

1.902160583104

from which the last two or three digits are uncertain.

Now let's find a couple of suitable approximations:

On squaring

1.902160583104

we get

3.61821488391

Notice the three digits after the decimal point resemble those of the golden-ratio. Let's add 2 to the built-in Phi constant:

1.61803398875

2 +   -->

1.61803398875

Let's divide the value obtained earlier by this one:

/   --->

1.00004999819

Not bad! But there's more. Let's take the square of B₂ again

3.61821488391

and divide it by the golden-ratio:

1.61803398875

/   -->

2.23617977686

Notice this is close to the square root of 5. So, let's square it:

2.23617977686

ENTER *   -->

5.00049999444

One 9 better! These allow for the following nice approximations:

B2 ~ Sqrt(1.00005*(Phi + 2)) = 1.90216058482

and

B2 ~ Sqrt(Phi*Sqrt(5.0005)) = 1.90216058363

The latter is a great mnemonic aid: Phi, a constant related to the square root of 5, and a 5-digit number, beginning and ending with 5 and zeroes in the middle.

P.S.: No bugs were found on the WP-34S when running the above program :-)

----------------------------------

P.S.: Some optimization:

001 LBL A		019 x!=? 02   		

002 STO 04		020 BACK 007		

003 # 048		021 x<> L			

004 SDR 002		022 x>? 04

005 °->G		023 SKIP 007

006 STO 01		024 ||

007 2   		025 y<> L		

008 +/-			026 1/x	

009 STO 02		027 STO+ 01		

010 5   		028 x<>y      

011 +   		029 RCL+ 03		

012 STO 03		030 BACK 015 		

013 RCL L		031 RCL 00		

014 DEC X		032 RCL 04

015 NEXTP		033 LN

016 FILL 	      	034 /

017 NEXTP		035 RCL+ 01  

018 -                   036 END       

2.64064726339 STO 00 

Emulador @ 1.86 GHz:

  10^5:      1.90216329186       (   6.1 s)

  10^6:      1.90191335333       (  63.1 s)

  10^7:      1.90218826322       ( 759.6 s)

2*10^7:      1.90217962170       (1692.3 s)

_{Edited to add P.S.}

Edited: 1 Jan 2013, 7:19 p.m. after one or more responses were posted

[Paul Dale]

Not all that relevant to the discussion but this was one of the constants that didn't make it into the 34S constant table :-)

There were others.

Still a nice short piece of methematics to get a very nice approximation.

- Pauli

[Gerson W. Barbosa]

What about an efficient NEXTTP function on the next WP-43S? :-)

Gerson.

[Paul Dale]

Just do it as user code. NEXTP is a keystroke program on the 34S after all :)

- Pauli