Approximate Brun's Constant

Brun's constant is the value to which the sum of the reciprocals of twin prime pairs (1/p and 1/(p+2) where p and p+2 are both prime) converges. It is approximately 1.902160583104.

Given a positive integer N, approximate Brun's constant by summing the reciprocals of the twin prime pairs where both primes in the pair are less than N, and output the approximation.

Rules

• N will be a positive integer within the representable range for your language.
• The output must be as accurate as possible to the true value, within the limits of your language's floating-point implementation, ignoring any potential issues due to floating-point arithmetic inaccuracies. If your language is capable of arbitrary-precision arithmetic, it must be at least as precise as IEEE 754 double-precision arithmetic.
• Alternatively, an exact fraction may be output in any consistent, unambiguous format.
• If a prime appears in multiple twin prime pairs (e.g. 5, part of both (3, 5) and (5, 7)), its reciprocal contributes to the sum each time.

Test Cases

2 -> 0
6 -> 0.5333333333333333
10 -> 0.8761904761904762
13 -> 0.8761904761904762
100 -> 1.3309903657190867
620 -> 1.4999706034568274
100000 -> 1.67279958482774

• Can an exact fraction be output? – LegionMammal978 Jan 22 '17 at 0:58
• @LegionMammal978 Yes, I'll clarify. – Mego Jan 22 '17 at 0:58
• Side note: the value 1.902160583104... for Brun's constant is only conjectured; not even the first significant figure has been rigorously computed (that is, it's not even known whether it's greater than or less than 2). – Greg Martin Jan 22 '17 at 1:05
• @GregMartin While that's true, it's also the best approximation we currently have. – Mego Jan 22 '17 at 1:06
• 5 is the only prime that appears in two prime pairs – Christian Sievers Jan 22 '17 at 3:41

Python 3, 7877757068 62 bytes

f=lambda n,k=3,m=1,j=0:k<n and-m%k*j*2/k+f(n,k+2,m*k**4,m%k/k)


Thanks to @xnor for golfing off 2 4 bytes and paving the way for 4 more!

Try it online!

Background

Recall that Wilson's theorem states that for all integers k > 1,

where a ≡ b (mod d) means that a - b is evenly divisible by d, i.e., a and b have the same residue when divided by d.

In Wilson Theorems for Double-, Hyper-, Sub- and Super-factorials, the authors prove generalizations for double factorials, on which this answer builds. The double factorial of an integer k ≥ 0 is defined by

Theorem 4 of the aforementioned paper states the following.

Elevating both sides of the congruences to the fourth power, we deduce that

for all odd primes p. Since 1!! = 1, the equivalence holds also for p = 2.

Now, doing the same to Wilson's theorem reveals that

Since

it follows that

whenever p is prime.

Now, let k be an odd, positive, composite integer. By definition, there exist integers a, b > 1 such that k = ab.

Since k is odd, so are a and b. Thus, both occur in the sequence 1, 3, …, k - 2 and

where | indicates divisibility.

Summing up, for all odd integers k > 1

where p(k) = 1 if k is prime and p(k) = 0 if k is composite.

How it works

When the function f is called with a single argument, k, m, and j are initialized as 3, 1, and 0.

Note that ((k - 2)!!)4 = 1!!4 = 1 = m. In fact, the equality m = ((k - 2)!!)4 will hold at all times. j is a float and will always be equal to ((k - 4)!!)4 % (k - 2) / (k - 2).

While k < n, the right argument of and will get evaluated. Since j = ((k - 4)!!)4 % (k - 2) / (k - 2), as proven in the first paragraph, j = 1/(k - 2) if k - 2 is prime and j = 0 if not. Likewise, since m % k = ((k - 2)!!)4 equals 1 if k is prime and 0 if not, -m % k = k - 1 if k is prime and -m % k = 0 if not. Therefore, -m%k*j*2/k evaluates to 2(k - 1)/(k(k - 2)) = ((k - 2) + k)/(k(k - 2)) = 1/k + 1/(k - 2) if the pair (k - 2, k) consists of twin primes and to 0 if not.

After computing the above, we add the result to the return value of the recursive call f(n,k+2,m*k**4,m%k/k). k gets incremented by 2 so it only takes odd values‡†, we multiply m by k4 since mk4 = ((k - 2)!!)4k4 = (k!!)4, and pass the current value of m % k / k – which equals 1/k if the "old" k is a prime and 0 if not – as parameter j to the function call.

Finally, once k is equal to or greater than n, f will return False and the recursion stops. The return value of f(n) will be the sum of all 1/k + 1/(k - 2) such (k - 2, k) is a twin prime pair and k < n, as desired.

The results from the Background paragraph hold only for odd integers. Since even integers cannot be twin primes, we can skip them safely.

• I think your expression is the same as m%k*(j/k+j/(k-2)). – xnor Jan 24 '17 at 3:05
• Yes, that works. Thanks! – Dennis Jan 24 '17 at 3:12
• Another save – xnor Jan 24 '17 at 3:25
• Nice observation that ((k-2)!!)^4 = p(k) modulo p for odd p. I haven't worked through your argument, but here's one I came up with (that might be the same in essence). Work modulo p in the set {1,2,..,p-1}, the evens are exactly the negatives of the odds. So, prod(odds) = ± prod(evens). Wilson's Theorem tells us that prod(all) = - p(k). Since prod(all) = prod(odds) * prod(evens) = prod(odds) * ± prod(evens), we have prod(odds)^2 = ±p(k) and so prod(odds)^4 = p(k)^2 = p(k). – xnor Jan 24 '17 at 3:43
• Nice! I tried expressing the sum as a single fraction, but computing part of it in j hadn't occurred to me. Thanks again! Your proof is a lot simpler than the one from the paper. – Dennis Jan 24 '17 at 4:05

Jelly, 15 14 bytes

’ÆRµ_2fµ+2;µİS


Try it online!

How it works

’ÆRµ_2fµ+2;µİS  Main link. Argument: n

’               Decrement; yield n-1.
ÆR             Prime range; yield all primes in [1, ..., n-1].
µ            New chain. Argument: r (prime range)
_2          Subtract 2 from all primes.
f         Filter; keep all p-2 that appear in r.
µ        New chain. Argument: t (filtered range)
+2      Add 2 to all primes in s.
;     Concatenate with s.
µ    New chain. Argument: t (twin primes)
İ   Take the inverses.
S  Sum.


Jelly, 16 14 bytes (with a little help from @Dennis)

’ÆRṡ2_/2+$$ÐḟFİS  Try it online! While trying to improve my previous answer, I thought up a totally different algorithm, and it comes in somewhat shorter. I'm using a different post for it, as is the standard here for an answer that uses a different technique. Dennis suggesting replacing _/2+$$Ðḟ with Iċ¥Ðf2; I'd completely forgotten about the possibility of a dyadic filter. As such, this algorithm now ties with the one that Dennis' answer used.

Explanation

’ÆRṡ2Iċ¥Ðf2FİS
’                  Decrement.
ÆR                Primes from 2 to the argument inclusive
(i.e. 2 to the original input exclusive).
ṡ2              Take overlapping slices of size 2.
Ðf         Keep only elements where the following is true:
¥           {the second parse of, which parses like this}
Iċ   2          the differences (I) contain (ċ) 2
F       Flatten.
İ      Take 1/x {for every list element}.
S     Sum.

• 2_/2+Ðḟ can become Iċ¥Ðf2. – Dennis Jan 22 '17 at 2:29

Brachylog, 17 bytes

{>I-₂:I{ṗ/₁}ᵐ}ᶠc+


Try it online!

This is the brand new version of Brachylog, with a shiny code page!

Explanation

{            }ᶠ        Find all valid outputs of the predicate in brackets
c+      Output is the sum of that list after flattening it

>I                    Input > I
-₂:I                The list [I-2, I]
{   }ᵐ          Map:
ṗ/₁              Must be prime and the output is its inverse


MATL, 16 bytes

liqZqtd2=)t2+h/s


Try it online!

Consider input 13 as an example.

l     % Push 1
%   STACK: 1
i     % Input N
%   STACK: 1, 13
q     % Subtract 1
%   STACK: 1, 12
Zq    % Primes up to that
%   STACK: 1, [2 3 5 7 11]
t     % Duplicate
%   STACK: 1, [2 3 5 7 11], [2 3 5 7 11]
d     % Consecutive differences
%   STACK: 1, [2 3 5 7 11], [1 2 2 4]
2=    % Compare with 2, element-wise
%   STACK: 1, [2 3 5 7 11], [0 1 1 0]
)     % Use as logical index to select elements from array
%   STACK: 1, [3 5]
t     % Duplicate
%   STACK: 1, [3 5], [3 5]
2+    % Add 2, element-wise
%   STACK: 1, [3 5], [5 7]
h     % Concatenate horizontally
%   STACK: 1, [3 5 5 7]
/     % Divide, element-wise
%   STACK: [0.3333 0.2 0.2 0.1429]
s     % Sum of array. Implicitly display
%   STACK: 0.8762


Mathematica, 48 47 bytes

Thanks to JungHwan Min for saving 1 byte!

If[PrimeQ/@(i&&(g=i-2)),1/i+1/g,0]~Sum~{i,#-1}&


Unnamed function taking a positive integer as input and returning an exact fraction; for example, If[PrimeQ/@(i&&(g=i-2)),1/i+1/g,0]~Sum~{i,#-1}&[10] returns 92/105.

If[PrimeQ/@(i&&(g=i-2)),1/i+1/g,0] tests whether both i and i-2 are prime, returning the sum of their reciprocals if so and 0 if not. ~Sum~{i,#-1}& then returns the sum of those contributions for all values of i less than the input.

Previous submission:

If[And@@PrimeQ@{i,g=i-2},1/i+1/g,0]~Sum~{i,#-1}&

• Now that's just spooky. I give up. ⚐ – LegionMammal978 Jan 22 '17 at 1:00
• I wondered if "exact fraction" meant Mathematica :) – Greg Martin Jan 22 '17 at 1:01
• -1 byte: If[PrimeQ/@(i&&(g=i-2)),1/i+1/g,0]~Sum~{i,#-1}& – JungHwan Min Jan 22 '17 at 1:37
• One can get an arbitrary-precision number by adding two bytes, N@, in front of the code. – JungHwan Min Jan 22 '17 at 1:56
• Nice golfing of the condition! It's true that N returns a decimal approximation to a real number; however, it requires extra bytes to display more than 6 sig figs or so, and no matter how many sig figs are displayed, it's still less accurate than the fraction itself. – Greg Martin Jan 22 '17 at 2:40

Octave, 45 bytes

@(n)sum(all(isprime(m=[h=3:n-1;h-2]))*m'.^-1)


Explanation:

m=[h=3:n-1;h-2]             generate an concatenate two ranges 3:n-1 and 1:n-3
rec=m'.^-1                  transpose and reciprocal
idx=all(isprime(m))         create a logical [0 1 ..] array  if both ranges are prime set 1 else set 0
sum1 = idx * rec            matrix multiplication(extrat elements with logical index and sum along the first dimension)
sum(sum1)                   sum along the second dimension


Try It Online!

JavaScript (ES6), 67 66 bytes

Saved 1 byte thanks to @Arnauld

f=n=>--n>1&&((p=x=>n%--x?p(x):x==1)(n)&&p(n-=2)&&1/n+++1/++n)+f(n)


Outputs false for test case 2, which is allowed by default.

Test snippet

f=n=>--n>1&&((p=x=>n%--x?p(x):x==1)(n)&&p(n-=2)&&1/n+++1/++n)+f(n)

g=n=>console.log("Input:",n,"Output:",f(n))

g(2)
g(6)
g(10)
g(13)
g(100)
g(620)

• I think 1/n+++1/++n saves a byte. – Arnauld Jan 22 '17 at 10:53
• @Arnauld Thanks. For some reason I didn't know that +++ doesn't always throw an error... – ETHproductions Jan 23 '17 at 21:02

05AB1E, 19 14 bytes (-5 @Emigna)

<LDpÏÍDpÏDÌ«zO


Try it online!

• <LDpÏÍDpÏDÌ«zO to save 5 bytes. – Emigna Jan 24 '17 at 8:31

Jelly, 19 bytes

’ÆRḊµ_Æp=2Tịµ_2;µİS


Try it online!

I have a feeling that this is improvable, but I can't immediately see how.

Explanation

’ÆRḊµ_Æp=2Tịµ_2;µİS
ÆR                  Generate all primes from 2 to n inclusive
’                    Subtract 1
Ḋ                 Remove first element
’ÆRḊ                 Generate all primes from 3 to n-1 exclusive

_Æp             Subtract the previous prime (i.e. calculate the prime gap)
=2           Compare to 2
Tị         Take elements of the input where the comparison is true
_Æp=2Tị         Filter a list of primes to the latter halves of prime pairs

_2      Subtract 2
;     Append
_2;     Append the list to the list with 2 subtracted from it
İ   Take reciprocals
S  Sum
İS  Take the sum of the reciprocals


The µ connect all these portions together pipeline-style, with each taking the output of the one before as its input.

Pyth - 2221 17 bytes

I think that refactoring will help.

scL1sf.APM_MTm,tt


Perl 6, 59 51 bytes

{sum 1 «/»grep((*-(2&0)).is-prime,^$_).flatmap:{$_-2,$_}} {sum 1 «/»grep(*.all.is-prime,(-2..*Z ^$_)).flat}


-2..* Z ^$_ zips the infinite list -2, -1, 0, 1, ... with the list 0, 1, ...$_-1 ($_ being the argument to the function), producing the list (-2, 0), (-1, 1), (0, 2), ..., ($_-3, $_-1). (Obviously none of these numbers less than 3 can be in a prime pair, but 3..* Z 5..^$_ is a few bytes longer, and none of the extra numbers are prime.)

The grep selects only those pairs where all (that is, both) numbers are prime, and the flat flattens them into a plain list of numbers.

«/» is the division hyperoperator; with the list on the right and 1 on the left, it turns the list of prime pairs into their reciprocals, which is then summed by sum.

Clojure, 147 bytes

(fn[n](let[p #(if(> % 2)(<(.indexOf(for[a(range 2 %)](mod % a))0)0))](reduce +(for[a(range 2 n)](if(and(p a)(p(- a 2)))(+(/ 1 a)(/ 1(- a 2)))0)))))


And Clojure comes dead last, as usual.

Ungolfed:

; Returns the primality of a number.
(defn prime? [n]
(if (> n 2)
(< (.indexOf (for [a (range 2 n)] (mod n a)) 0) 0)))

; Calculates the actual Brun's Constant. ' (Stupid highlighter)
(defn brunsconst [n]
; Adds all of the entries together
(reduce
+
; For a in range(2, n):
(for [a (range 2 n)]
(let [b (- a 2)]
; If both a and a-2 are prime:
(if (and (prime? a) (prime? b))
; Place (1/a + 1/a-2) on the array, else 0
(+ (/ 1 a) (/ 1 b)) 0)))))


Julia 0.4, 48 46 bytes

!n=sum(1./[(p=n-1|>primes)∩(p-2);p∩(p+2)])


Try it online!

Bash + GNU utilities, 86 85 bytes

for((k=4;k<$1;k++,j=k-2)){ [ factor$k $j|wc -w = 4 ]&&x=$x+1/$k+1/$j;};bc -l<<<0$x  Try it online! Constructs a large arithmetic expression and then feeds it to bc -l to evaluate it. Edit: Mistakenly left in a$(...) pair from an old version with nested command substitution; changed to backticks to save a byte.

APL NARS, 216 bytes, 108 chars

  r←z n;h;i;k;v
i←0⋄n-←1⋄h←1+⍳n-1⋄→B
A:k←i⊃h⋄h←k∪(0≠k∣h)/h
B:→A×⍳(⍴h)≥i+←1
r←+/÷(v-2),v←(h=1⌽h+2)/h


this would use the "Crivello di Eratostene" for find the sublist in 1..arg of request primes. Test:

  z¨2 6 10 13 100 620
0 0.5333333333 0.8761904762 0.8761904762 1.330990366 1.499970603
z 100000
1.672799585