Meissel–Mertens constant: Difference between revisions
(Added Algol 68) |
(→{{header|Wren}}: Added Analytic method.) |
||
Line 351: | Line 351: | ||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
===Summation method=== |
|||
{{libheader|Wren-math}} |
{{libheader|Wren-math}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
Line 386: | Line 387: | ||
100,000,000 0.261497212859 |
100,000,000 0.261497212859 |
||
105,097,565 0.261497212858 |
105,097,565 0.261497212858 |
||
</pre> |
|||
===Analytic method=== |
|||
{{trans|PARI/GP}} |
|||
{{libheader|Wren-gmp}} |
|||
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). |
|||
<syntaxhighlight lang="ecmascript">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))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
0.261497212847642783...70842383659092665508 |
|||
</pre> |
</pre> |
Revision as of 22:34, 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
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
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
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