[Question]
  [
[Brun's constant](https://en.wikipedia.org/wiki/Brun%27s_theorem) is the value to which the sum of the reciprocals of [twin prime](https://en.wikipedia.org/wiki/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
```
  
[Answer]
# [Python 3](https://docs.python.org/3/), ~~78~~ ~~77~~ ~~75~~ ~~70~~ ~~68~~ 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!](https://tio.run/nexus/python3#FYxBCsIwFETXeorZhDbpr00TyaIY7xIpgTbmK0XPH7@LYQbe8FqOz1QfawJTiZ5qnGmPdik3RuJ1rKqY3bipDLmXx@CommLMlQRMRbf8OsDYGI4QCLOV@H/LCM4u59P72PjTd8p/Md6hLnPIHRTEBlFqrdsP "Python 3 – TIO Nexus")
### Background
Recall that [Wilson's theorem](https://en.wikipedia.org/wiki/Wilson%27s_theorem "Wilson's theorem - Wikipedia") states that for all integers **k > 1**,
![](https://i.stack.imgur.com/bmCcO.png)
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](https://arxiv.org/abs/1302.3676 "[1302.3676] Wilson Theorems for Double-, Hyper-, Sub- and Super-factorials"), the authors prove generalizations for double factorials, on which this answer builds. The [double factorial](https://en.wikipedia.org/wiki/Double_factorial "Double factorial - Wikipedia") of an integer **k ≥ 0** is defined by
![](https://i.stack.imgur.com/iCGD4.png)
Theorem 4 of the aforementioned paper states the following.
![](https://i.stack.imgur.com/rhz9x.png)
Elevating both sides of the congruences to the fourth power, we deduce that
![](https://i.stack.imgur.com/NbNii.png)
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
![](https://i.stack.imgur.com/JPDKa.png)
Since
![](https://i.stack.imgur.com/Zt2g7.png)
it follows that
![](https://i.stack.imgur.com/MkN6e.png)
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
![](https://i.stack.imgur.com/7nyHV.png)
where **|** indicates divisibility.
Summing up, for all odd integers **k > 1**
![](https://i.stack.imgur.com/V4Wzw.png)
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.
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~15~~ 14 bytes
```
’ÆRµ_2fµ+2;µİS
```
[Try it online!](https://tio.run/nexus/jelly#@/@oYebhtqBDW@ON0g5t1TayPrT1yIbg/0f3HG5/1LTG/f9/Ix0FMx0FQwMgNgbRQIaZkQEA "Jelly – TIO Nexus")
### 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.
```
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~16~~ 14 bytes (with a little help from @Dennis)
```
’ÆRṡ2_/2+$$ÐḟFİS
```
[Try it online!](https://tio.run/nexus/jelly#@/@oYebhtqCHOxcaeR7pPrT08IQ0I7cjG4L/RxvpmOkYGugYGgNJAx0zI4PYw@2Pmtb8BwA "Jelly – TIO Nexus")
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.
```
[Answer]
# [Brachylog](https://github.com/JCumin/Brachylog), 17 [bytes](https://github.com/JCumin/Brachylog/wiki/Code-page)
```
{>I-₂:I{ṗ/₁}ᵐ}ᶠc+
```
[Try it online!](https://tio.run/nexus/brachylog2#ASMA3P//ez5JLeKCgjpJe@G5ly/igoF94bWQfeG2oGMr//82MjD/Wg "Brachylog – TIO Nexus")
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
```
[Answer]
# [MATL](https://github.com/lmendo/MATL), 16 bytes
```
liqZqtd2=)t2+h/s
```
[Try it online!](https://tio.run/nexus/matl#@5@TGVVYkmJkq1lipJ2hX/z/v5mRAQA "MATL – TIO Nexus")
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
```
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), ~~19~~ 14 bytes (-5 @Emigna)
```
<LDpÏÍDpÏDÌ«zO
```
[Try it online!](https://tio.run/nexus/05ab1e#@2/j41JwuN/l0FKjwOJDK4qBTCPtQ6ur/P//NzQGAA "05AB1E – TIO Nexus")
[Answer]
# 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}&
```
[Answer]
# 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!](https://tio.run/#blfF2)
[Answer]
## 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](http://meta.codegolf.stackexchange.com/a/9067/12012).
### 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)
```
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 19 bytes
```
’ÆRḊµ_Æp=2Tịµ_2;µİS
```
[Try it online!](https://tio.run/nexus/jelly#@/@oYebhtqCHO7oObY0/3FZgaxTycHc3kG1kfWjrkQ3B/6ONdMx0DA10DI2BpIGOmZFB7OH2R01r/gMA "Jelly – TIO Nexus")
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.
[Answer]
# Pyth - ~~22~~ ~~21~~ 17 bytes
I think that refactoring will help.
```
scL1sf.APM_MTm,tt
```
[Test Suite](http://pyth.herokuapp.com/?code=scL1sf.APM_MTm%2Ctt&test_suite=1&test_suite_input=2%0A6%0A10%0A13%0A100%0A620&debug=0).
[Answer]
# [Julia 0.4](http://julialang.org/), ~~48~~ 46 bytes
```
!n=sum(1./[(p=n-1|>primes)∩(p-2);p∩(p+2)])
```
[Try it online!](https://tio.run/nexus/julia#@6@YZ1tcmqthqKcfrVFgm6drWGNXUJSZm1qs@ahjpUaBrpGmdQGYpW2kGav5Py2/SCFPITNPQcNIx0zH0EDH0BhIGuiYGRlocnE6ALXmlaRpKKkalyro2imo6hmapcXkKenk6SjmaXKl5qX8BwA "Julia 0.4 – TIO Nexus")
[Answer]
# [Perl 6](http://perl6.org/), ~~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`.
[Answer]
# 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)))))
```
[Answer]
# 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
```
[Answer]
# [Husk](https://github.com/barbuz/Husk), 14 bytes
```
ṁ\Σfo=2F-Ẋefṗŀ
```
[Try it online!](https://tio.run/##ASIA3f9odXNr///huYFczqNmbz0yRi3huoplZuG5l8WA////MTAw "Husk – 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
```
[Answer]
# [Pari/GP](http://pari.math.u-bordeaux.fr/), 50 bytes
```
n->forprime(p=s=0.,n,isprime(q=p-2)&&s+=1/p+1/q);s
```
[Try it online!](https://tio.run/##JYrBCsMwDEN/xfRQYpq0SQa9FPdHxg49LCOweW6yy74@czuB0JOQbCW7h7QE1Nit6V2k5NfdCFXyo2Wb63/YSVzEvq8DhUmGMO241LaJPL@Gwa2gN/4odkfpIBlGtHCNFmYLwasvRyrM0Z@kumH7AQ "Pari/GP – Try It Online")
[Answer]
# [Vyxal](https://github.com/Vyxal/Vyxal) `sḋ`, ~~13~~ 12 bytes
```
~æ2l'÷ε2=;fĖ
```
[Try it online](https://vyxal.pythonanywhere.com/#WyJz4biLIiwiIiwifsOmMmwnw7fOtTI9O2bEliIsIiIsIjEwMCJd) or [verify all test cases](https://vyxal.pythonanywhere.com/#WyJz4biLQSIsIiIsIn7DpjJsJ8O3zrUyPTtmxJYiLCIiLCIyXG42XG4xMFxuMTNcbjEwMFxuNjIwIl0=).
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
```
[Answer]
# [Bash](https://www.gnu.org/software/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!](https://tio.run/nexus/bash#@5@WX6ShkW1rYp1to2Jona2trZNlm61rpKlZrRCtkJCWmFySX6Sgkq2gklVTnqygW56gYKtgohCrplZhq1Khbaivkg0isqxrrZOA0jk2NjYGKhX///83AwA "Bash – TIO Nexus")
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.*
]