Meissel–Mertens constant

Task

Calculate Meissel–Mertens constant up to a precision your language can handle.

Motivation

Analogous to Euler's constant, which is important in determining the sum of reciprocal natural numbers, Meissel-Mertens' constant is important in calculating the sum of reciprocal primes.

Example

We consider the finite sum of reciprocal natural numbers:

1 + 1/2 + 1/3 + 1/4 + 1/5 ... 1/n

this sum can be well approximated with:

log(n) + E

where E denotes Euler's constant: 0.57721...

log(n) denotes the natural logarithm of n.

Now consider the finite sum of reciprocal primes:

1/2 + 1/3 + 1/5 + 1/7 + 1/11 ... 1/p

this sum can be well approximated with:

log( log(p) ) + M

where M denotes Meissel-Mertens constant: 0.26149...

See

Details in the Wikipedia article: Meissel–Mertens constant

ALGOL 68

Sieving odd primes only, not using extended precision.
Note that to run this with Algol 68G, a large heap size must be specified, e.g.: by adding -heap 256M to the command line.

BEGIN # compute an approximation to the Meissel-Mertens constant             #
    # construct a sieve of odd primes                                        #
    [ 0 : 10 000 000 ]BOOL primes;
    BEGIN
        FOR i TO UPB primes DO primes[ i ] := TRUE  OD;
        INT ip := 1;
        FOR i WHILE i + ( ip +:= 2 ) <= UPB primes DO
            IF primes[ i ] THEN
                FOR s FROM i + ip BY ip TO UPB primes DO primes[ s ] := FALSE OD
            FI
        OD
    END;
    # sum the reciprocals of the primes                                      #
    INT       p count := 1;
    INT       last p  := 0;
    LONG REAL sum     := long ln( 0.5 ) + 0.5;
    INT       p       := 1;
    INT       p10     := 10;
    # Euler's constant from the wikipedia, truncated for LONG REAL           #
    LONG REAL eulers constant = 0.5772156649015328606 # 0651209008240243104215933593992 #;
    FOR i TO UPB primes DO
        p +:= 2;
        IF primes[ i ] THEN
            LONG REAL rp = 1 / LENG p;
            sum     +:= long ln( 1 - rp ) + rp;
            p count +:= 1;
            last p   := p;
            IF p count = p10 THEN
                print( ( "after ", whole( p count, -8 ), " primes, the approximation is: "
                       , fixed( sum + eulers constant, -14, 12 )
                       , ", last prime considered: ", whole( last p, 0 )
                       , newline
                       )
                     );
                p10 := IF p10 < 1 000 000 THEN p10 * 10 ELSE p10 + 1 000 000 FI
            FI
        FI   
    OD;
    print( ( "after ", whole( p count, -8 ), " primes, the approximation is: "
           , fixed( sum + eulers constant, -14, 12 )
           , ", last prime considered: ", whole( last p, 0 )
           , newline
           )
         )
END

Output:

after       10 primes, the approximation is: 0.265160104017, last prime considered: 29
after      100 primes, the approximation is: 0.261624626173, last prime considered: 541
after     1000 primes, the approximation is: 0.261503563617, last prime considered: 7919
after    10000 primes, the approximation is: 0.261497594498, last prime considered: 104729
after   100000 primes, the approximation is: 0.261497238454, last prime considered: 1299709
after  1000000 primes, the approximation is: 0.261497214692, last prime considered: 15485863
after  1270607 primes, the approximation is: 0.261497214255, last prime considered: 19999999

FreeBASIC

On my CPU i5 it takes about 3.5 minutes

'#include "isprime.bas"

Const As Double Euler = 0.57721566490153286

Dim As Double m = 0
For x As Ulongint = 2 To 1e8
    If isPrime(x) Then m += Log(1-(1/x)) + (1/x)
Next x
Print "MM ="; Euler + m
Sleep

Output:

MM = 0.2614972131104154

Go

Translation of: Wren

Library: Go-rcu

package main

import (
    "fmt"
    "math"
    "rcu"
)

func contains(a []int, f int) bool {
    for _, e := range a {
        if e == f {
            return true
        }
    }
    return false
}

func main() {
    const euler = 0.57721566490153286
    primes := rcu.Primes(1 << 31)
    pc := len(primes)
    sum := 0.0
    fmt.Println("Primes added         M")
    fmt.Println("------------  --------------")
    for i, p := range primes {
        rp := 1.0 / float64(p)
        sum += math.Log(1.0-rp) + rp
        c := i + 1
        if (c%1e7) == 0 || c == pc {
            fmt.Printf("%11s   %0.12f\n", rcu.Commatize(c), sum+euler)
        }
    }
}

Output:

Primes added         M
------------  --------------
 10,000,000   0.261497212987
 20,000,000   0.261497212912
 30,000,000   0.261497212889
 40,000,000   0.261497212878
 50,000,000   0.261497212871
 60,000,000   0.261497212867
 70,000,000   0.261497212864
 80,000,000   0.261497212862
 90,000,000   0.261497212861
100,000,000   0.261497212859
105,097,565   0.261497212858

J

   Euler=: 0.57721566490153286
   MM=: {{ Euler + +/ (+ ^.@-.)@% p:i. y }}
   0j12 ": MM 1e7
0.261497212987

Julia

Off by one in the 11th digit after 10^8 primes.

using Base.MathConstants  # sets constant γ = 0.5772156649015...
using Primes

""" Approximate the Meissel-Mertons constant. """
function meissel_mertens(iterations = 100_000_000)
    return mapreduce(p ->(d = 1/p; log(1 - d) + d), +, primes(prime(iterations))) + γ
end

@show meissel_mertens(100_000_000) # meissel_mertens(100000000) = 0.2614972128591237

PARI/GP

Summation method

{
MM(t)=
  my(s=0);
  forprime(p = 2, t,
    s += log(1.-1./p)+1./p
  );
  Euler+s
};

Output:

Running 10^9 summations to get 9 valid digits.

? \p10
   realprecision = 19 significant digits (10 digits displayed)
? MM(1e9)
?
%1 = 0.2614972129
?
? ##
  ***   last result: cpu time 1min, 18,085 ms, real time 1min, 18,094 ms.
?

Analytic method

The Analytic method requires high precision calculation of Riemann zeta function.

{
Meissel_Mertens(d)=
  default(realprecision, d);
  my(prec = default(realprecision), z = 0, y = 0, q);
  forprime(p = 2 , 7, 
    z += log(1.-1./p)+1./p
  );
  for(k = 2, prec,
    q = 1;
    forprime(p = 2, 7, 
      q *= 1.-p^-k
    );
    y += moebius(k)*log(zeta(k)*q)/k
  );
  Euler+z+y
};

Output:

1000 valid digits.

? Meissel_Mertens(1001)
%1 = 0.26149721284764278375542683860869585905156664826119920619206421392492451089736820971414263143424665105161772887648602199778339032427004442454348740197238640666194955709392581712774774211985258807266272064144464232590023543105177232173925663229980314763831623758149059290382284758265972363422015971458785446941586825460538918007031787714156680620570605257601785334398970354507934530971953511716888598019955346947142883673537117910619342522616975101911159537244599605203558051780574237201332999961769676911386909654186249097435916294862238555389898241954857937738258646582212506260380084370067541379219020626760709633535981989783010762417792511961619355361391684002933280522289185167238258837930443067100391254985761418536020400457460311825670423438456551983202200477824746954606715454777572171338072595463648319687279859427306787306509669454587505942593547068846408425666008833035029366514525328713339609172639368543886291288200447611698748441593459920236225093315001729474600911978170842383659092665509
? ##
  ***   last result: cpu time 283 ms, real time 284 ms.

Agrees with PARI/GP except for the very last digit being 8 instead of 9. Note that playing with precision and/or iterations (which seems to differ quite wildly) gave utterly incorrect results... not that I actually comprehend what the algorithm is doing. The 1,003 includes 2 from the "0.", fairly obviously.

"0.261497212847642783...70842383659092665508 (1,003 digits)"

Wren

Library: Wren-math

Library: Wren-fmt

It will be seen that this is converging to the correct answer though we'd need to add a lot more primes to obtain a valid 11th digit.

import "./math" for Int
import "./fmt" for Fmt

var euler = 0.57721566490153286
var primes = Int.primeSieve(2.pow(31))
var pc = primes.count
var sum = 0
var c = 0
System.print("Primes added         M")
System.print("------------  --------------")
for (p in primes) {
    var rp = 1/p
    sum = (1-rp).log + rp + sum
    c = c + 1
    if ((c % 1e7) == 0 || c == pc) Fmt.print("$,11d   $0.12f", c, sum + euler)
}

Output:

Primes added         M
------------  --------------
 10,000,000   0.261497212987
 20,000,000   0.261497212912
 30,000,000   0.261497212889
 40,000,000   0.261497212878
 50,000,000   0.261497212871
 60,000,000   0.261497212867
 70,000,000   0.261497212864
 80,000,000   0.261497212862
 90,000,000   0.261497212861
100,000,000   0.261497212859
105,097,565   0.261497212858