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
11l
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
ALGOL W
begin % compute an approximation to the Meissel-Mertens constant %
integer MAX_PRIME;
MAX_PRIME := 10000000;
begin
logical array primes ( 0 :: MAX_PRIME );
begin % construct a sieve of odd primes %
integer i, ip;
for i := 1 until MAX_PRIME do primes( i ) := true;
ip := 3;
i := 1;
while i + ip <= MAX_PRIME do begin
if primes( i ) then begin
for s := i + ip step ip until MAX_PRIME do primes( s ) := false
end if_primes__i ;
ip := ip + 2;
i := i + 1
end while_i_plus_ip_le_MAX_PRIME
end;
begin % sum the reciprocals of the primes %
integer pCount, lastP, p, p10;
long real sum, eulersConstant;
procedure showProgress ;
write( s_w := 0, i_w := 8
, r_format := "A", r_w := 14, r_d := 12
, "after ", pCount, " primes, the approximation is: "
, ( sum + eulersConstant ), ", last prime considered: "
, i_w := 1
, lastP
);
pCount := 1;
lastP := 0;
sum := longLn( 0.5 ) + 0.5;
p := 1;
p10 := 10;
% Euler's constant from the wikipedia, truncated for long real %
eulersConstant := 0.5772156649015328606 % 0651209008240243104215933593992 %;
for i := 1 until MAX_PRIME do begin
p := p + 2;
if primes( i ) then begin
long real rp;
rp := 1 / long p;
sum := sum + longLn( 1 - rp ) + rp;
pCount := pCount + 1;
lastP := p;
if pCount = p10 then begin
showProgress;
p10 := if p10 < 1000000 THEN p10 * 10 else p10 + 1000000
end if_pCount_eq_p10
end if_primes__i
end for_i ;
showProgress
end
end
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
Arturo
meisselMertens: function [depth][
Euler: 0.57721566490153286
m: (1//2) + ln 1-1//2
loop range.step:2 3 depth 'x ->
if prime? x ->
m: m + (1//x) + ln 1-1//x
return m + Euler
]
print meisselMertens 10000000
- Output:
0.2614972157776471
BASIC
BASIC256
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
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
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
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
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
C++
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
std::vector<double> list_prime_reciprocals(const int32_t& limit) {
const int32_t half_limit = ( limit % 2 == 0 ) ? limit / 2 : 1 + limit / 2;
std::vector<bool> composite(half_limit);
for ( int32_t i = 1, p = 3; i < half_limit; p += 2, ++i ) {
if ( ! composite[i] ) {
for ( int32_t a = i + p; a < half_limit; a = a + p ) {
composite[a] = true;
}
}
}
std::vector<double> result(composite.size());
result[0] = 0.5;
for ( int32_t i = 1, p = 3; i < half_limit; p += 2, ++i ) {
if ( ! composite[i] ) {
result.emplace_back(1.0 / p);
}
}
return result;
}
int main() {
std::vector<double> prime_reciprocals = list_prime_reciprocals(100000000);
const double euler = 0.577215664901532861;
double sum = 0.0;
for ( double reciprocal : prime_reciprocals ) {
sum += reciprocal + log(1.0 - reciprocal);
}
const double meissel_mertens = euler + sum;
std::cout << "The Meissel-Mertens constant = " << std::setprecision(8) << meissel_mertens << std::endl;
}
- Output:
The Meissel-Mertens constant = 0.26149721
Delphi
This program doesn't test the limits of numbers in Delphi. The limit is the processing time. The time to process 10^7 is 3 minutes. The time to process 10^8 is nearly two hours. As a result, process beyond 10^8 is pretty much impossible and a different strategy would be required
function IsPrime(N: int64): boolean;
{Fast, optimised prime test}
var I,Stop: int64;
begin
if (N = 2) or (N=3) then Result:=true
else if (n <= 1) or ((n mod 2) = 0) or ((n mod 3) = 0) then Result:= false
else
begin
I:=5;
Stop:=Trunc(sqrt(N+0.0));
Result:=False;
while I<=Stop do
begin
if ((N mod I) = 0) or ((N mod (I + 2)) = 0) then exit;
Inc(I,6);
end;
Result:=True;
end;
end;
function GetNextPrime(var Start: integer): integer;
{Get the next prime number after Start}
{Start is passed by "reference," so the
{original variable is incremented}
begin
repeat Inc(Start)
until IsPrime(Start);
Result:=Start;
end;
function MeisselMertens(Depth: integer; Prog: TProgress): extended;
{Calculate MM value up a certain Depth}
var I,P: integer;
const Euler = 0.57721566490153286;
begin
Result:=0;
P:=1;
for I:=1 to Depth do
begin
P:=GetNextPrime(P);
Result:=Result+Ln(1-(1/P)) + (1/P);
if Assigned(Prog) and ((I mod 10000)=0) then
HandleProgress(MulDiv(100,I,Depth));
end;
Result:=Result+Euler;
end;
procedure ShowMeisselMertens(Memo: TMemo; Prog: TProgress);
var I,IT,Digits: integer;
var M,Last: extended;
begin
Memo.Lines.Add('Primes Digits M');
Memo.Lines.Add('-----------------------------------------');
Last:=0;
IT:=1;
{Calculate MM to specified Power of 10}
for I:=1 to 7 do
begin
IT:=IT*10;
M:=MeisselMertens(IT,Prog);
{Calculated Digits of accuracy}
Digits:=Trunc(abs(Log(abs(M-Last))));
Memo.Lines.Add(Format('10^%2d %7d %25.18f',[I,Digits,M]));
Last:=M
end;
end;
- Output:
Primes Digits M ----------------------------------------- 10^ 1 0 0.265160104017311294 10^ 2 2 0.261624626172936147 10^ 3 3 0.261503563616850571 10^ 4 5 0.261497594498432488 10^ 5 6 0.261497238454297260 10^ 6 7 0.261497214692305277 10^ 7 8 0.261497212987251183 10^ 8 9 0.261497212858582276
Dart
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
EasyLang
fastfunc isprim num .
if num mod 2 = 0
if num = 2
return 1
.
return 0
.
i = 3
while i <= sqrt num
if num mod i = 0
return 0
.
i += 2
.
return 1
.
func log x .
return log10 x / log10 2.7182818284590452354
.
euler = 0.5772156649
for x = 2 to 1e6
if isprim x = 1
m += log (1 - (1 / x)) + (1 / x)
.
.
numfmt 11 0
print "mm = " & euler + m
- Output:
mm = 0.26149724673
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
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
Implementation
Euler=: 0.57721566490153286
MM=: {{ Euler + +/ (+ ^.@-.)@% p: i. y }}
Task example
0j13 ": MM 1e8
0.2614972128591
That said, a more literal implementation, like
mm=: (% 10x(^#)10#.inv]) 26149721284764278375542683860869585905156664826119920619206421392x
would be more precise:
0j65":mm
0.26149721284764278375542683860869585905156664826119920619206421392
Java
import java.util.ArrayList;
import java.util.BitSet;
import java.util.List;
public final class MeisselMertensConstant {
public static void main(String[] aArgs) {
List<Double> primeReciprocals = listPrimeReciprocals(1_000_000_000);
final double euler = 0.577_215_664_901_532_861;
double sum = 0.0;
for ( double reciprocal : primeReciprocals ) {
sum += reciprocal + Math.log(1.0 - reciprocal);
}
final double meisselMertens = euler + sum;
System.out.println(String.format("%s%.9f", "The Meissel-Mertens constant = ", meisselMertens));
}
private static List<Double> listPrimeReciprocals(int aLimit) {
BitSet sieve = new BitSet(aLimit + 1);
sieve.set(2, aLimit + 1);
for ( int i = 2; i <= Math.sqrt(aLimit); i = sieve.nextSetBit(i + 1) ) {
for ( int j = i * i; j <= aLimit; j += i ) {
sieve.clear(j);
}
}
List<Double> result = new ArrayList<Double>(sieve.cardinality());
for ( int i = 2; i >= 0; i = sieve.nextSetBit(i + 1) ) {
result.add(1.0 / i);
}
return result;
}
}
- Output:
The Meissel-Mertens constant = 0.261497213
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
Maxima
meissel_mertens:%gamma+lsum(log(1-(1/i))+(1/i),i,primes(2,10000000)),numer;
- Output:
0.2614972157776471
Nim
import std/[math, strformat, strutils]
proc initPrimes(N: static int): seq[int] =
## Initialize the list of primes.
const M = 2 * N - 1
var composite = newSeq[bool](N)
composite[0] = true # 1 is not prime.
# Conversions from index to value and value to index.
template index(n: int): int = (n - 1) shr 1
template value(idx: int): int = idx shl 1 + 1
# Fill the sieve.
var n = 3
while n * n <= M:
if not composite[n.index]:
for k in countup(n * n, M, 2 * n):
composite[k.index] = true
inc n, 2
# Build list of primes.
result = @[2]
for idx in 0..composite.high:
if not composite[idx]:
result.add idx.value
const N = 2^30
let primes = initPrimes(N)
echo "Primes added M"
echo "──────────── ──────────────"
const γ = 0.57721566490153286 # Euler–Mascheroni constant.
let primeCount = primes.len
var sum = 0.0
var count = 0
for p in primes:
let rp = 1 / p
sum += ln(1 - rp) + rp
inc count
if count mod 10_000_000 == 0 or count == primeCount:
echo &"{insertSep($count):>11} {sum+γ:.12}"
- 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
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
use v5.36;
use ntheory qw(forprimes);
my $s;
forprimes { $s += log(1 - 1/$_)+1/$_ } 1e9;
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,"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
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)"
Python
#!/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
# 20221011 Raku programming solution
my $s;
.is-prime and $s += log(1-1/$_)+1/$_ for 2 .. 10**8;
say $s + .57721566490153286
- Output:
0.26149721310571383
REXX
Libraries: How to use
Library: Numbers
Library: Functions
Library: Constants
Library: Settings
Library: Abend
Library: Sequences
Below program shows 3 methods to calculate the Meissel-Mertens constant.
- BruteForce: Check all odd integers 3...999999 on primality (Miller-Rabin) and use defining sum.
- UsingSieve: Get all primes below 1000000 (Sieve) and use defining sum.
- Analytic: Fast converging sum formula given by Wolfram MathWorld.
include Settings
say version; say 'Meissel-Mertens constant'; say
arg n
if n = '' then
n = 16
numeric digits n
fact. = 0
call Time('r'); a = BruteForce(); e = Format(Time('e'),,3)
say 'BruteForce' a '('e 'seconds)'
call Time('r'); a = UsingSieve(); e = Format(Time('e'),,3)
say 'UsingSieve' a '('e 'seconds)'
call Time('r'); a = Analytic(); e = Format(Time('e'),,3)
say 'Analytic ' a '('e 'seconds)'
call Time('r'); a = TrueValue(); e = Format(Time('e'),,3)
say 'True value' a '('e 'seconds)'
exit
BruteForce:
procedure expose glob.
numeric digits Digits()+2
y = 0.5-Ln(2)
do n = 3 by 2 to 1000000
if Prime(n) then do
q = 1/n; t = Ln(1-q)+q; y = y+t
end
end
y = Euler()+y
numeric digits Digits()-2
return y+0
UsingSieve:
procedure expose glob. prim.
numeric digits Digits()+2
n = Primes(1000000); y = 0
do i = 1 to n
q = 1/prim.prime.i; t = Ln(1-q)+q; y = y+t
end
y = Euler()+y
numeric digits Digits()-2
return y+0
Analytic:
procedure expose glob. fact.
numeric digits Digits()+2
y = 0; v = 0
do n = 2 to 1000
t = Moebius(n) * Ln(Zeta(n)) / n
if t <> 0 then do
y = y+t
if y = v then
leave
v = y
end
end
y = Euler()+y
numeric digits Digits()-2
return y+0
TrueValue:
procedure
return 0.261497212847642783755426838608695859051566648261199206192064213924924510897368209714142631434246651051617+0
include Constants
include Functions
include Numbers
- Output:
Runs with 16 digits: REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022 Meissel-Mertens constant BruteForce 0.2614972467355791 (19.263 seconds) UsingSieve 0.2614972467355791 (1.510 seconds) Analytic 0.2614972128476428 (0.040 seconds) True value 0.2614972128476428 (0.000 seconds) Runs with 111 digits: REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022 111 digits Meissel-Mertens constant BruteForce 0.261497246735500132677260976437201892720087109609945002636762560229390985475526828746177041776354307483184476679 (57.470 seconds) UsingSieve 0.261497246735500132677260976437201892720087109609945002636762560229390985475526828746177041776354307483184476679 (1.830 seconds) Analytic 0.261497212847642783755426838608695859051566648261199206192064213924924510897368209714142631434246651051617728876 (48.760 seconds) True value 0.261497212847642783755426838608695859051566648261199206192064213924924510897368209714142631434246651051617 (0.000 seconds)
Using primes only 7 correct decimals are achieved. UsingSieve is faster, but requires extra memory to store all primes.
In spite of the expensive functions Moebius and Zeta, the analytic solution outperforms. The Analytic results suggest that the OEIS entry is truncated, not rounded.
RPL
BPRIM?
is defined at Primality by trial division#RPL
≪ → n ≪ 0.5 3 n FOR j IF j R→B BPRIM? THEN j INV + END 2 STEP n LN LN - ≫ ≫ 'MMCON' STO
10000 MMCON
- Output:
1: 0.262733140866
Sidef
var sum = 0
1e7.each_prime {|p|
with (1f/p) {|t|
sum += (log(1 - t) + t)
}
}
say sum+Num.EulerGamma
- Output:
0.26149721577767111119422410228297206467931376306
Wren
Summation method
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
Analytic method
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
t.uiSub(1, r).log.add(r)
z.add(t, z)
}
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)
y.add(t, y)
}
return Mpf.euler.add(z).add(y)
}
Fmt.print("$20a", meisselMertens.call(3300).toString(1001))
- Output:
0.261497212847642783...70842383659092665508
XPL0
Runs in 98 seconds on Pi4.
func IsPrime(N); \Return 'true' if N is prime
int N, D;
[if N < 2 then return false;
if (N&1) = 0 then return N = 2;
if rem(N/3) = 0 then return N = 3;
D:= 5;
while D*D <= N do
[if rem(N/D) = 0 then return false;
D:= D+2;
if rem(N/D) = 0 then return false;
D:= D+4;
];
return true;
];
def Euler = 0.57721566490153286;
real M; int P;
[M:= 0.;
for P:= 2 to 100_000_000 do
if IsPrime(P) then
M:= M + Ln(1. - 1./float(P)) + 1./float(P);
Format(1, 16);
Text(0, "MM = "); RlOut(0, Euler + M); CrLf(0);
]
- Output:
MM = 0.2614972131104150
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.26149721310.2614972128476427837554268386086958590515666482611992061920642139249245108973682097141426314342466510516177288764860219977833903242700444245434874019723864066619495570939258171277477421198525880726627206414446423259002354310517723217392566322998031476383162375814905929038228475826597236342201597145878545286546864785868143588282690675811212349253390275898516648963722862532834592270193787663814788872688899729598789511243611212105643183719483028350694889140412496558830385604700978671612478906659270558671778566494523172201279234383779159082832574101983245266507469788929642499636163109342230499137582782542115136325204939772977424363327903550519767079810438569344262666067369673777630199473914034464437524475856770642386781282897163202273071362276688321001599406964213181471899136557333922267566346325514106895289284293459552062339943989806700523705145085201774716630458161117596603698964857884739475619061322196278674275162487770072794776455644351959504021182511495417580620912301758280597008866059 0.26149721284764278375542683860869585905156664826119920619206421392492451089736820971414263143424665105161772887648602199778339032427004442454348740197238640666194955709392581712774774211985258807266272064144464232590023543105177232173925663229980314763831623758149059290382284758265972363422015971458785446941586825460538918007031787714156680620570605257601785334398970354507934530971953511716888598019955346947142883673537117910619342522616975101911159537244599605203558051780574237201332999961769676911386909654186249097435916294862238555389898241954857937738258646582212506260380084370067541379219020626760709633535981989783010762417792511961619355361391684002933280522289185167238258837930443067100391254985761418536020400457460311825670423438456551983202200477824746954606715454777572171338072595463648319687279859427306787306509669454587505942593547068846408425666008833035029366514525328713339609172639368543886291288200447611698748441593459920236225093315001729474600911978170842383659092665508