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:

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

{
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
};