# 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? Jan 22, 2017 at 0:58
• @LegionMammal978 Yes, I'll clarify.
– user45941
Jan 22, 2017 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). Jan 22, 2017 at 1:05
• @GregMartin While that's true, it's also the best approximation we currently have.
– user45941
Jan 22, 2017 at 1:06
• 5 is the only prime that appears in two prime pairs Jan 22, 2017 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 = \frac{\mod(\left(\left(k - 4\right)!!\right)^4 , k - 2)}{k - 2}\$$, as proven in the first paragraph, $$\j = \frac{1}{k - 2}\$$ if $$\k - 2\$$ is prime and $$\j = 0\$$ if not. Likewise, since $$\\mod(m,k) = \left(\left(k - 2\right)!!\right)^4\$$ equals $$\1\$$ if $$\k\$$ is prime and $$\0\$$ if not, $$\\mod(-m,k) = k - 1\$$ if $$\k\$$ is prime and $$\\mod(-m , k) = 0\$$ if not. Therefore, -m%k*j*2/k evaluates to $$\\frac{2\left(k - 1\right)}{k\left(k - 2\right)} = \frac{\left(k - 2\right) + k}{k(k - 2)} = \frac{1}{k} + \frac{1}{k - 2}\$$ if the pair $$\\left(k - 2, k\right)\$$ 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 $$\k^4\$$ since $$\mk^4 = \left(\left(k - 2\right)!!\right)^4k^4 = \left(k!!\right)^4\$$, and pass the current value of $$\\frac{\mod(m,k)}{k}\$$ – which equals $$\\frac{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 $$\\frac{1}{k} + \frac{1}{k - 2}\$$ such $$\\left(k - 2, k\right)\$$ 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, 2017 at 3:05
• Yes, that works. Thanks! Jan 24, 2017 at 3:12
• Another save
– xnor
Jan 24, 2017 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, 2017 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. Jan 24, 2017 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. Jan 22, 2017 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]
%   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


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

<LDpÏÍDpÏDÌ«zO


Try it online!

• <LDpÏÍDpÏDÌ«zO to save 5 bytes. Jan 24, 2017 at 8:31

# 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}& 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. ⚐ Jan 22, 2017 at 1:00
• I wondered if "exact fraction" meant Mathematica :) Jan 22, 2017 at 1:01
• -1 byte: If[PrimeQ/@(i&&(g=i-2)),1/i+1/g,0]~Sum~{i,#-1}& Jan 22, 2017 at 1:37
• One can get an arbitrary-precision number by adding two bytes, N@, in front of the code. Jan 22, 2017 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. Jan 22, 2017 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. Jan 22, 2017 at 10:53
• @Arnauld Thanks. For some reason I didn't know that +++ doesn't always throw an error... Jan 23, 2017 at 21:02

# 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


# Julia 0.4, 48 46 bytes

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


Try it online!

# 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)))))


# 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


# Husk, 14 bytes

ṁ\Σfo=2F-Ẋefṗŀ


Try it online!

### Explanation

ṁ\Σfo=2F-Ẋefṗŀ
ŀ   lowered range, 0..n-1
fṗ    filter primes
Ẋe      get each consecutive pair
f   F-        filter where difference
o=2          equals 2
Σ              concat, list of pairs to flattened list
ṁ\               map 1/n and then sum


# Pari/GP, 50 bytes

n->forprime(p=s=0.,n,isprime(q=p-2)&&s+=1/p+1/q);s


Try it online!

# Vyxalsḋ, 13 12 bytes

~æ2l'÷ε2=;fĖ


Explanation:

~             # Filter the range 0 to n-1 by:
æ            #   Is prime
2l          # Get the overlapping groups of 2
'    ;    # Filter those by:
÷ε       #   Absolute difference of the list
2=     #   Equals 2
f   # Flatten to make one list
Ė  # Take the reciprocal of each
# Sum with the s flag
# Convert to decimal format with the ḋ flag


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