Meissel–Mertens constant: Difference between revisions
m
→{{header|Phix}}: added analytical/gmp version
mNo edit summary |
m (→{{header|Phix}}: added analytical/gmp version) |
||
Line 152:
24,304,615 0.261497212899 (accurate to 10 d.p.)
(actual value 0.26149721284764278375542683860869)
</pre>
===analytical/gmp===
{{trans|PARI/GP}}
<!--<syntaxhighlight lang="phix">-->
<span style="color: #008080;">without</span> <span style="color: #008080;">js</span> <span style="color: #000080;font-style:italic;">-- no mpfr_zeta[_ui]() in mpfr.js, as yet, or mpfr_log() or mpfr_const_euler(), for that matter</span>
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">moebius</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">prime_factors</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #004600;">true</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">odd</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">))?-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">:+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">Meissel_Mertens</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_set_default_precision</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (d decimal places)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">5</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">p</span> <span style="color: #008080;">in</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">7</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">do</span>
<span style="color: #000080;font-style:italic;">-- z += log(1-1/p)+1/p</span>
<span style="color: #7060A8;">mpfr_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_div_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_si_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">d</span> <span style="color: #008080;">do</span> <span style="color: #000080;font-style:italic;">-- (see note)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">moebius</span><span style="color: #0000FF;">(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">m</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">mpfr_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">p</span> <span style="color: #008080;">in</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">7</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">do</span>
<span style="color: #000080;font-style:italic;">-- q *= 1-power(p,-k)</span>
<span style="color: #7060A8;">mpfr_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_pow_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_si_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000080;font-style:italic;">-- y += moebius(k)*log(zeta(k)*q)/k</span>
<span style="color: #7060A8;">mpfr_zeta_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_div_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #000000;">rp</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000080;font-style:italic;">-- res := EULER+z+y</span>
<span style="color: #000000;">mpfr_const_euler</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">Meissel_Mertens</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1001</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">?</span><span style="color: #7060A8;">shorten</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mpfr_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1001</span><span style="color: #0000FF;">))</span>
<!--</syntaxhighlight>-->
{{out}}
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.
<pre>
"0.261497212847642783...70842383659092665508 (1,003 digits)"
</pre>
| |||