Problem Statement

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

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.

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

Let's start with the solution:

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);
]