Meissel–Mertens constant

Meissel–Mertens constant
You are encouraged to solve this task according to the task description, using any language you may know.

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

Go

Translation of: Wren

Library: Go-rcu

Unlike the Wren example, slowly converging towards the correct answer.

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.

Phix

Converges very slowly, I started running out of memory (from generating so many primes) before it showed any signs of getting stuck.

with javascript_semantics -- (but perhaps a bit too slow)
constant mmc = 0.2614972128476427837554268386086958590516,
       smmc = "0.2614972128476427837554268386086958590516"
atom t = 0.57721566490153286, dpa = 1, p10=1
integer pn = 1, adp = 0
string st = "0", fmt = "%.0f"
printf(1,"Primes added         M\n")
printf(1,"------------  --------------\n")
while adp<=10 do
    atom rp = 1/get_prime(pn)
    t += log(1-rp) + rp
    if (t-mmc)<dpa then
        string ft = sprintf(fmt,trunc(t*p10)/p10) -- (as below)
        if ft=st then
            --
            -- We have to artifically calculate a "truncated t", aka tt, 
            -- to prevent say 0.2..299[>5..] being automatically rounded 
            -- by printf() to 0.2..300, otherwise it just "looks wrong".
            --
            atom tt = trunc(t*1e12)/1e12
            printf(1,"%,11d   %0.12f (accurate to %d d.p.)\n", {pn, tt, adp})
            adp += 1
            fmt = sprintf("%%.%df",adp)
            st = smmc[1..2+adp]
            dpa /= 10
            p10 *= 10
        end if
    end if
    pn += 1
end while
printf(1,"(actual value 0.26149721284764278375542683860869)\n")

Output:

(Couldn't be bothered with making it an inner loop to upgrade the "accurate to 4dp" to 5dp)

Primes added         M
------------  --------------
          1   0.384068484341 (accurate to 0 d.p.)
          3   0.288793158252 (accurate to 1 d.p.)
          6   0.269978901636 (accurate to 2 d.p.)
         38   0.261990332075 (accurate to 3 d.p.)
      1,940   0.261499998883 (accurate to 4 d.p.)
      1,941   0.261499997116 (accurate to 5 d.p.)
      5,471   0.261497999951 (accurate to 6 d.p.)
     34,891   0.261497299999 (accurate to 7 d.p.)
    303,447   0.261497219999 (accurate to 8 d.p.)
  9,246,426   0.261497212999 (accurate to 9 d.p.)
 24,304,615   0.261497212899 (accurate to 10 d.p.)
(actual value 0.26149721284764278375542683860869)

analytical/gmp

Translation of: PARI/GP

without js -- no mpfr_zeta[_ui]() in mpfr.js, as yet, or mpfr_log() or mpfr_const_euler(), for that matter
include mpfr.e
function moebius(integer n)
    if n=1 then return 1 end if
    sequence f = prime_factors(n,true)
    for i=2 to length(f) do
        if f[i] = f[i-1] then return 0 end if
    end for
    return iff(odd(length(f))?-1:+1)
end function

function Meissel_Mertens(integer d)
    mpfr_set_default_precision(-d) -- (d decimal places)
    mpfr {res,rp,z,y,q} = mpfr_inits(5)
    for p in {2,3,5,7} do
        -- z += log(1-1/p)+1/p
        mpfr_set_si(rp,1)
        mpfr_div_si(rp,rp,p)
        mpfr_add(z,z,rp)
        mpfr_si_sub(rp,1,rp)
        mpfr_log(rp,rp)
        mpfr_add(z,z,rp)
    end for
    for k=2 to d do -- (see note)
        integer m = moebius(k)
        if m then
            mpfr_set_si(q,1)
            for p in {2,3,5,7} do
                -- q *= 1-power(p,-k)
                mpfr_set_si(rp,p)
                mpfr_pow_si(rp,rp,-k)
                mpfr_si_sub(rp,1,rp)
                mpfr_mul(q,q,rp)
            end for
            -- y += moebius(k)*log(zeta(k)*q)/k
            mpfr_zeta_ui(rp,k)
            mpfr_mul(rp,rp,q)
            mpfr_log(rp,rp)
            mpfr_div_si(rp,rp,k)
            mpfr_mul_si(rp,rp,m)
            mpfr_add(y,y,rp)
        end if
    end for
    -- res := EULER+z+y
    mpfr_const_euler(res)
    mpfr_add(z,z,y)
    mpfr_add(res,res,z)
    return res
end function

mpfr m = Meissel_Mertens(1001)
?shorten(mpfr_get_str(m,10,1001))

Output:

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