Meissel–Mertens constant: Difference between revisions
→{{header|Wren}}: Added Analytic method.
(Added Algol 68) |
(→{{header|Wren}}: Added Analytic method.) |
||
Line 351:
=={{header|Wren}}==
===Summation method===
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
Line 386 ⟶ 387:
100,000,000 0.261497212859
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>
| |||