Meissel–Mertens constant

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

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

11l

Translation of: Python
```F primes_up_to_limit(Int limit)
[Int] r
I limit >= 2
r.append(2)

V isprime = [1B] * ((limit - 1) I/ 2)
V sieveend = Int(sqrt(limit))
L(i) 0 .< isprime.len
I isprime[i]
Int p = i * 2 + 3
r.append(p)
I i <= sieveend
L(j) ((p * p - 3) >> 1 .< isprime.len).step(p)
isprime[j] = 0B
R r

V euler = 0.57721566490153286
V m = 0.0
L(x) primes_up_to_limit(10'000'000)
m += log(1 - (1 / x)) + (1 / x)

print(‘MM = #.16’.format(euler + m))```
Output:
```MM = 0.2614972157776471
```

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
```

BASIC

BASIC256

Translation of: FreeBASIC
```Euler = 0.5772156649
m = 0
for x = 2 to 1e6    # more prime numbers do not add more precision
if isPrime(x) then m += log(1-(1/x)) + (1/x)
next x
print "MM = "; Euler + m
print Euler
end

function isPrime(v)
if v < 2 then return False
if v mod 2 = 0 then return v = 2
if v mod 3 = 0 then return v = 3
d = 5
while d * d <= v
if v mod d = 0 then return False else d += 2
end while
return True
end function
```
Output:
`MM = 0.26149724673`

Run BASIC

Translation of: FreeBASIC
```function isPrime(n)
if n < 2       then isPrime = 0 : goto [exit]
if n = 2       then isPrime = 1 : goto [exit]
if n mod 2 = 0 then isPrime = 0 : goto [exit]
isPrime = 1
for i = 3 to int(n^.5) step 2
if n mod i = 0 then isPrime = 0 : goto [exit]
next i
[exit]
end function

e = 0.5772156

for x = 2 to 100000    ' more prime numbers do not add more precision
if isPrime(x) then m = m + log(1-(1/x)) + (1/x)
next x
print "MM = "; e + m```
Output:
`MM = 0.261274`

PureBasic

Translation of: FreeBASIC
```Procedure isPrime(v.i)
If     v <= 1    : ProcedureReturn #False
ElseIf v < 4     : ProcedureReturn #True
ElseIf v % 2 = 0 : ProcedureReturn #False
ElseIf v < 9     : ProcedureReturn #True
ElseIf v % 3 = 0 : ProcedureReturn #False
Else
Protected r = Round(Sqr(v), #PB_Round_Down)
Protected f = 5
While f <= r
If v % f = 0 Or v % (f + 2) = 0
ProcedureReturn #False
EndIf
f + 6
Wend
EndIf
ProcedureReturn #True
EndProcedure

OpenConsole()
Euler.d = 0.5772156649 ;0153286

For x.i = 2 To 1e8
If isPrime(x)
m.d = m + Log(1-(1/x)) + (1/x)
EndIf
Next x
PrintN("MM = " + StrD(Euler + m))
PrintN(#CRLF\$ + "--- terminado, pulsa RETURN---"): Input()
CloseConsole()
```
Output:
`MM = 0.2614972129`

True BASIC

Translation of: FreeBASIC
```FUNCTION isPrime (n)
IF n = 2 THEN
LET isPrime = 1
ELSEIF n <= 1 OR REMAINDER(n, 2) = 0 THEN
LET isPrime = 0
ELSE
LET isPrime = 1
FOR i = 3 TO INT(SQR(n)) STEP 2
IF REMAINDER(n, i) = 0 THEN
LET isPrime = 0
EXIT FUNCTION
END IF
NEXT i
END IF
END FUNCTION

LET e = .5772156649

FOR x = 2 to 1e6                  !more prime numbers do not add more precision
IF isPrime(x) = 1 THEN
LET m = m + LOG(1-(1/x)) + (1/x)
END IF
NEXT x
PRINT "MM ="; e + m
END
```
Output:
`MM = 0.26149725`

Yabasic

Translation of: FreeBASIC
```e = 0.5772156

for x = 2 to 1e6    // more prime numbers do not add more precision
if isPrime(x)  m = m + log(1-(1/x)) + (1/x)
next x
print "MM = ", e + m
end

sub isPrime(v)
if v < 2  return False
if mod(v, 2) = 0  return v = 2
if mod(v, 3) = 0  return v = 3
d = 5
while d * d <= v
if mod(v, d) = 0 then return False else d = d + 2 : fi
wend
return True
end sub
```
Output:
`MM = 0.261497`

Dart

Translation of: FreeBASIC
```import 'dart:math';

bool isPrime(var n) {
if (n <= 1) return false;
if (n == 2) return true;
for (var i = 2; i <= sqrt(n); i++) if (n % i == 0) return false;
return true;
}

void main() {
const double euler = 0.57721566490153286;

double m = 0.0;
for (var x = 2; x <= 1e8; x++)
if (isPrime(x)) m += log(1 - (1 / x)) + (1 / x);

print('MM = \${euler + m}');
}
```
Output:
`MM = 0.2614972131057144`

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("------------  --------------")
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

Implementation

```Euler=: 0.57721566490153286
MM=: {{ Euler + +/ (+ ^.@-.)@% p: i. y }}
```

```   0j13 ": MM 1e8
0.2614972128591
```

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:

10 valid digits, prime summation up to 10^9.

```? \p 12
realprecision = 19 significant digits (12 digits displayed)
? MM(1e9)
?
%1 = 0.261497212874
?
? ##
***   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. For some reason the last digit is in some rare cases (rounded) wrong. Whatever, the last digit is per definition the one with the lowest weight.

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

Perl

Library: ntheory
```use v5.36;
use ntheory 'is_prime';

my \$s;
is_prime \$_ and \$s += log(1 - 1/\$_)+1/\$_ for 2 .. 10**9;
say my \$result = \$s + .57721566490153286;
```
Output:
`0.261497212871049`

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,"------------  --------------\n")
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})
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_si_sub(rp,1,rp)
mpfr_log(rp,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)
end if
end for
-- res := EULER+z+y
mpfr_const_euler(res)
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)"
```

Python

Translation of: FreeBASIC
```#!/usr/bin/python

from math import log

def isPrime(n):
for i in range(2, int(n**0.5) + 1):
if n % i == 0:
return False
return True

if __name__ == '__main__':
Euler = 0.57721566490153286
m = 0
for x in range(2, 10_000_000):
if isPrime(x):
m += log(1-(1/x)) + (1/x)

print("MM =", Euler + m)
```
Output:
`MM = 0.26149721577764207`

Raku

Translation of: FreeBASIC
```# 20221011 Raku programming solution

my \$s;
.is-prime and \$s += log(1-1/\$_)+1/\$_ for 2 .. 10**8;
say \$s + .57721566490153286
```
Output:
`0.26149721310571383`

Wren

Summation method

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("------------  --------------")
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
```

Analytic method

Translation of: PARI/GP
Library: Wren-gmp

This agrees with the Phix entry that the 1,000th digit after the decimal point is '8' after rounding (according to OEIS it's actually '7' but the next one is '7' also).

```import "./gmp" for Mpf
import "./math" for Int
import "./fmt" for Fmt

var isSquareFree = Fn.new { |n|
var i = 2
while (i * i <= n) {
if (n%(i*i) == 0) return false
i = (i > 2) ? i + 2 : i + 1
}
return true
}

var mu = Fn.new { |n|
if (n < 1) Fiber.abort("Argument must be a positive integer")
if (n == 1) return 1
var sqFree = isSquareFree.call(n)
var factors = Int.primeFactors(n)
if (sqFree && factors.count % 2 == 0) return 1
if (sqFree) return -1
return 0
}

var meisselMertens = Fn.new { |d|
Mpf.defaultPrec = d
var z = Mpf.zero
var y = Mpf.zero
var r = Mpf.new()
var q = Mpf.new()
var t = Mpf.new()
var m = Mpf.new()
for (p in [2, 3, 5, 7]) {
r.setUi(p).inv
}
for (k in 2..d) {
q.setUi(1)
for (p in [2, 3, 5, 7]) {
r.setUi(p).inv
t.uiSub(1, r.pow(k))
q.mul(t)
}
m.setSi(mu.call(k))
t.zetaUi(k).mul(q).log.mul(m).div(k)
}
}

Fmt.print("\$20a", meisselMertens.call(3300).toString(1001))
```
Output:
```0.261497212847642783...70842383659092665508
```

Mathematica / Wolfram Language

```PrimeNumbers = Select[Range[100000000], PrimeQ[#] &]; (*all primes in the first 100 000 000 numbers, this takes a toll on my computer's CPU and RAM*)
MM = N[Total[Log[1 - 1/PrimeNumbers] + 1/PrimeNumbers] + EulerGamma, 10] (*Calculating it up to a precision of 10, this is correct up to 8 digits*)
AnalyticMMto305 = N[EulerGamma + Sum[MoebiusMu[n]/n Log[Zeta[n]], {n, 2, 1000}], 1000] (*Precise up to 305 digits*)
AnalyticMM = N[EulerGamma + Sum[MoebiusMu[n]/n Log[Zeta[n]], {n, 2, 10000}], 1001] (*Precise up to at least 1000 digits*)
```
Output:
```0.2614972131
0.2614972128476427837554268386086958590515666482611992061920642139249245108973682097141426314342466510516177288764860219977833903242700444245434874019723864066619495570939258171277477421198525880726627206414446423259002354310517723217392566322998031476383162375814905929038228475826597236342201597145878545286546864785868143588282690675811212349253390275898516648963722862532834592270193787663814788872688899729598789511243611212105643183719483028350694889140412496558830385604700978671612478906659270558671778566494523172201279234383779159082832574101983245266507469788929642499636163109342230499137582782542115136325204939772977424363327903550519767079810438569344262666067369673777630199473914034464437524475856770642386781282897163202273071362276688321001599406964213181471899136557333922267566346325514106895289284293459552062339943989806700523705145085201774716630458161117596603698964857884739475619061322196278674275162487770072794776455644351959504021182511495417580620912301758280597008866059
0.26149721284764278375542683860869585905156664826119920619206421392492451089736820971414263143424665105161772887648602199778339032427004442454348740197238640666194955709392581712774774211985258807266272064144464232590023543105177232173925663229980314763831623758149059290382284758265972363422015971458785446941586825460538918007031787714156680620570605257601785334398970354507934530971953511716888598019955346947142883673537117910619342522616975101911159537244599605203558051780574237201332999961769676911386909654186249097435916294862238555389898241954857937738258646582212506260380084370067541379219020626760709633535981989783010762417792511961619355361391684002933280522289185167238258837930443067100391254985761418536020400457460311825670423438456551983202200477824746954606715454777572171338072595463648319687279859427306787306509669454587505942593547068846408425666008833035029366514525328713339609172639368543886291288200447611698748441593459920236225093315001729474600911978170842383659092665508

```