Meissel–Mertens constant: Difference between revisions
(→{{header|Go}}: Removed preamble as Wren example now fixed.) |
(Meissel–Mertens constant in FreeBASIC) |
||
Line 38:
:* Details in the Wikipedia article: [https://en.wikipedia.org/wiki/Meissel%E2%80%93Mertens_constant Meissel–Mertens constant]
<br /><br />
=={{header|FreeBASIC}}==
On my CPU i5 it takes about 3.5 minutes
<syntaxhighlight lang="freebasic">'#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</syntaxhighlight>
{{out}}
<pre>MM = 0.2614972131104154</pre>
=={{header|Go}}==
|
Revision as of 10:19, 5 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 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
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
Euler=: 0.57721566490153286
MM=: {{ Euler + +/ (+ ^.@-.)@% p:i. y }}
0j12 ": MM 1e7
0.261497212987
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:
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)
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)"
Wren
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