Meissel–Mertens constant: Difference between revisions
Line 50:
};</syntaxhighlight>
{{out}}
Running 10^9 summations to get 9 valid digits
<pre>? \p10
realprecision = 19 significant digits (10 digits displayed)
Line 79:
};</syntaxhighlight>
{{out}}
1000 valid digits
<pre>? Meissel_Mertens(1001)
%1 = 0.26149721284764278375542683860869585905156664826119920619206421392492451089736820971414263143424665105161772887648602199778339032427004442454348740197238640666194955709392581712774774211985258807266272064144464232590023543105177232173925663229980314763831623758149059290382284758265972363422015971458785446941586825460538918007031787714156680620570605257601785334398970354507934530971953511716888598019955346947142883673537117910619342522616975101911159537244599605203558051780574237201332999961769676911386909654186249097435916294862238555389898241954857937738258646582212506260380084370067541379219020626760709633535981989783010762417792511961619355361391684002933280522289185167238258837930443067100391254985761418536020400457460311825670423438456551983202200477824746954606715454777572171338072595463648319687279859427306787306509669454587505942593547068846408425666008833035029366514525328713339609172639368543886291288200447611698748441593459920236225093315001729474600911978170842383659092665509
|
Revision as of 07:36, 4 October 2022
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 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
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)
Wren
Wren's only native number type is a 64 bit float (15 or 16 digits accuracy) and it appears from the following results that, no matter how many primes we use, we're not going to be able to improve on 8 valid digits for M with this particular number type.
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 = sum + (1-rp).log + rp
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.261497213008 20,000,000 0.261497213008 30,000,000 0.261497213009 40,000,000 0.261497213008 50,000,000 0.261497213009 60,000,000 0.261497213009 70,000,000 0.261497213009 80,000,000 0.261497213009 90,000,000 0.261497213009 100,000,000 0.261497213009 105,097,565 0.261497213009