Meissel–Mertens constant: Difference between revisions
Added Easylang
(Added Easylang) |
|||
| (8 intermediate revisions by 4 users not shown) | |||
Line 286:
<pre>MM = 0.261497</pre>
=={{header|C++}}==
<syntaxhighlight lang="c++">
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
std::vector<double> list_prime_reciprocals(const int32_t& limit) {
const int32_t half_limit = ( limit % 2 == 0 ) ? limit / 2 : 1 + limit / 2;
std::vector<bool> composite(half_limit);
for ( int32_t i = 1, p = 3; i < half_limit; p += 2, ++i ) {
if ( ! composite[i] ) {
for ( int32_t a = i + p; a < half_limit; a = a + p ) {
composite[a] = true;
}
}
}
std::vector<double> result(composite.size());
result[0] = 0.5;
for ( int32_t i = 1, p = 3; i < half_limit; p += 2, ++i ) {
if ( ! composite[i] ) {
result.emplace_back(1.0 / p);
}
}
return result;
}
int main() {
std::vector<double> prime_reciprocals = list_prime_reciprocals(100'000'000);
const double euler = 0.577'215'664'901'532'861;
double sum = 0.0;
for ( double reciprocal : prime_reciprocals ) {
sum += reciprocal + log(1.0 - reciprocal);
}
const double meissel_mertens = euler + sum;
std::cout << "The Meissel-Mertens constant = " << std::setprecision(8) << meissel_mertens << std::endl;
}
</syntaxhighlight>
{{ out }}
<pre>
The Meissel-Mertens constant = 0.26149721
</pre>
=={{header|Delphi}}==
Line 402 ⟶ 447:
{{out}}
<pre>MM = 0.2614972131057144</pre>
=={{header|EasyLang}}==
{{trans|BASIC256}}
<syntaxhighlight>
fastfunc isprim num .
if num mod 2 = 0
if num = 2
return 1
.
return 0
.
i = 3
while i <= sqrt num
if num mod i = 0
return 0
.
i += 2
.
return 1
.
func log x .
return log10 x / log10 2.7182818284590452354
.
euler = 0.5772156649
for x = 2 to 1e6
if isprim x = 1
m += log (1 - (1 / x)) + (1 / x)
.
.
numfmt 11 0
print "mm = " & euler + m
</syntaxhighlight>
{{out}}
<pre>
mm = 0.26149724673
</pre>
=={{header|FreeBASIC}}==
Line 490 ⟶ 572:
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.ArrayList;
import java.util.BitSet;
import java.util.List;
public final class MeisselMertensConstant {
Line 507 ⟶ 584:
double sum = 0.0;
for ( double reciprocal : primeReciprocals ) {
sum +=
}
final double
System.out.println(String.format("%s%.9f", "The Meissel-Mertens constant = ",
}
Line 518 ⟶ 595:
sieve.set(2, aLimit + 1);
for ( int j = i * i; j <= aLimit; j += i ) {
sieve.clear(j);
Line 537 ⟶ 613:
{{ out }}
<pre>
The Meissel-Mertens constant = 0.
</pre>
Line 553 ⟶ 629:
@show meissel_mertens(100_000_000) # meissel_mertens(100000000) = 0.2614972128591237
</syntaxhighlight>
=={{header|Maxima}}==
<syntaxhighlight lang="maxima">
meissel_mertens:%gamma+lsum(log(1-(1/i))+(1/i),i,primes(2,10000000)),numer;
</syntaxhighlight>
{{out}}
<pre>
0.2614972157776471
</pre>
=={{header|Nim}}==
{{trans|Wren}}
<syntaxhighlight lang="Nim">import std/[math, strformat, strutils]
proc initPrimes(N: static int): seq[int] =
## Initialize the list of primes.
const M = 2 * N - 1
var composite = newSeq[bool](N)
composite[0] = true # 1 is not prime.
# Conversions from index to value and value to index.
template index(n: int): int = (n - 1) shr 1
template value(idx: int): int = idx shl 1 + 1
# Fill the sieve.
var n = 3
while n * n <= M:
if not composite[n.index]:
for k in countup(n * n, M, 2 * n):
composite[k.index] = true
inc n, 2
# Build list of primes.
result = @[2]
for idx in 0..composite.high:
if not composite[idx]:
result.add idx.value
const N = 2^30
let primes = initPrimes(N)
echo "Primes added M"
echo "──────────── ──────────────"
const γ = 0.57721566490153286 # Euler–Mascheroni constant.
let primeCount = primes.len
var sum = 0.0
var count = 0
for p in primes:
let rp = 1 / p
sum += ln(1 - rp) + rp
inc count
if count mod 10_000_000 == 0 or count == primeCount:
echo &"{insertSep($count):>11} {sum+γ:.12}"
</syntaxhighlight>
{{out}}
<pre>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
</pre>
=={{header|PARI/GP}}==
Line 603 ⟶ 751:
{{libheader|ntheory}}
<syntaxhighlight lang="perl" line>use v5.36;
use ntheory
my $s;
say my $result = $s + .57721566490153286;</syntaxhighlight>
{{out}}
Line 773 ⟶ 921:
<pre>
1: 0.262733140866
</pre>
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">var sum = 0
1e7.each_prime {|p|
with (1f/p) {|t|
sum += (log(1 - t) + t)
}
}
say sum+Num.EulerGamma</syntaxhighlight>
{{out}}
<pre>
0.26149721577767111119422410228297206467931376306
</pre>
Line 780 ⟶ 941:
{{libheader|Wren-fmt}}
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.
<syntaxhighlight lang="
import "./fmt" for Fmt
Line 817 ⟶ 978:
{{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="
import "./math" for Int
import "./fmt" for Fmt
| |||