Approximate Brun's Constant

Question

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

Answer 1

Python 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!

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)!!)⁴ = 1!!⁴ = 1 = m. In fact, the equality m = ((k - 2)!!)⁴ will hold at all times. j is a float and will always be equal to ((k - 4)!!)⁴ % (k - 2) / (k - 2).

While k < n, the right argument of and will get evaluated. Since j = ((k - 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)!!)⁴ 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 k⁴ since mk⁴ = ((k - 2)!!)⁴k⁴ = (k!!)⁴, 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.}

Answer 2

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.

Answer 3

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.

Answer 4

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

Answer 5

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

Answer 6

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 7

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)

Answer 8

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!

Answer 9

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

<LDpÏÍDpÏDÌ«zO

Try it online!

Answer 10

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.

Answer 11

Pyth - 22 21 17 bytes

I think that refactoring will help.

scL1sf.APM_MTm,tt

Test Suite.

Answer 12

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.

Answer 13

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 14

Julia 0.4, 48 46 bytes

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

Try it online!

Answer 15

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.

Answer 16

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