Meissel–Mertens constant: Difference between revisions
Added Easylang
SqrtNegInf (talk | contribs) (Added Perl) |
(Added Easylang) |
||
| (21 intermediate revisions by 11 users not shown) | |||
Line 38:
:* Details in the Wikipedia article: [https://en.wikipedia.org/wiki/Meissel%E2%80%93Mertens_constant Meissel–Mertens constant]
<br /><br />
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">
F primes_up_to_limit(Int limit)
[Int] r
I limit >= 2
r.append(2)
V isprime = [1B] * ((limit - 1) I/ 2)
V sieveend = Int(sqrt(limit))
L(i) 0 .< isprime.len
I isprime[i]
Int p = i * 2 + 3
r.append(p)
I i <= sieveend
L(j) ((p * p - 3) >> 1 .< isprime.len).step(p)
isprime[j] = 0B
R r
V euler = 0.57721566490153286
V m = 0.0
L(x) primes_up_to_limit(10'000'000)
m += log(1 - (1 / x)) + (1 / x)
print(‘MM = #.16’.format(euler + m))
</syntaxhighlight>
{{out}}
<pre>
MM = 0.2614972157776471
</pre>
=={{header|ALGOL 68}}==
Line 100 ⟶ 133:
</pre>
=={{header|Arturo}}==
<syntaxhighlight lang="arturo">meisselMertens: function [depth][
Euler: 0.57721566490153286
m: (1//2) + ln 1-1//2
loop range.step:2 3 depth 'x ->
if prime? x ->
m: m + (1//x) + ln 1-1//x
return m + Euler
]
print meisselMertens 10000000</syntaxhighlight>
{{out}}
<pre>0.2614972157776471</pre>
=={{header|BASIC}}==
Line 236 ⟶ 285:
{{out}}
<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}}==
{{works with|Delphi|6.0}}
{{libheader|SysUtils,StdCtrls}}
This program doesn't test the limits of numbers in Delphi. The limit is the processing time. The time to process 10^7 is 3 minutes. The time to process 10^8 is nearly two hours. As a result, process beyond 10^8 is pretty much impossible and a different strategy would be required
<syntaxhighlight lang="Delphi">
function IsPrime(N: int64): boolean;
{Fast, optimised prime test}
var I,Stop: int64;
begin
if (N = 2) or (N=3) then Result:=true
else if (n <= 1) or ((n mod 2) = 0) or ((n mod 3) = 0) then Result:= false
else
begin
I:=5;
Stop:=Trunc(sqrt(N+0.0));
Result:=False;
while I<=Stop do
begin
if ((N mod I) = 0) or ((N mod (I + 2)) = 0) then exit;
Inc(I,6);
end;
Result:=True;
end;
end;
function GetNextPrime(var Start: integer): integer;
{Get the next prime number after Start}
{Start is passed by "reference," so the
{original variable is incremented}
begin
repeat Inc(Start)
until IsPrime(Start);
Result:=Start;
end;
function MeisselMertens(Depth: integer; Prog: TProgress): extended;
{Calculate MM value up a certain Depth}
var I,P: integer;
const Euler = 0.57721566490153286;
begin
Result:=0;
P:=1;
for I:=1 to Depth do
begin
P:=GetNextPrime(P);
Result:=Result+Ln(1-(1/P)) + (1/P);
if Assigned(Prog) and ((I mod 10000)=0) then
HandleProgress(MulDiv(100,I,Depth));
end;
Result:=Result+Euler;
end;
procedure ShowMeisselMertens(Memo: TMemo; Prog: TProgress);
var I,IT,Digits: integer;
var M,Last: extended;
begin
Memo.Lines.Add('Primes Digits M');
Memo.Lines.Add('-----------------------------------------');
Last:=0;
IT:=1;
{Calculate MM to specified Power of 10}
for I:=1 to 7 do
begin
IT:=IT*10;
M:=MeisselMertens(IT,Prog);
{Calculated Digits of accuracy}
Digits:=Trunc(abs(Log(abs(M-Last))));
Memo.Lines.Add(Format('10^%2d %7d %25.18f',[I,Digits,M]));
Last:=M
end;
end;
</syntaxhighlight>
{{out}}
<pre>
Primes Digits M
-----------------------------------------
10^ 1 0 0.265160104017311294
10^ 2 2 0.261624626172936147
10^ 3 3 0.261503563616850571
10^ 4 5 0.261497594498432488
10^ 5 6 0.261497238454297260
10^ 6 7 0.261497214692305277
10^ 7 8 0.261497212987251183
10^ 8 9 0.261497212858582276
</pre>
=={{header|Dart}}==
Line 259 ⟶ 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 336 ⟶ 561:
<syntaxhighlight lang="j"> 0j13 ": MM 1e8
0.2614972128591</syntaxhighlight>
That said, a more literal implementation, like
<syntaxhighlight lang=J>mm=: (% 10x(^#)10#.inv]) 26149721284764278375542683860869585905156664826119920619206421392x</syntaxhighlight>
would be more precise:
<syntaxhighlight lang=J> 0j65":mm
0.26149721284764278375542683860869585905156664826119920619206421392</syntaxhighlight>
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.ArrayList;
import java.util.BitSet;
import java.util.List;
public final class MeisselMertensConstant {
public static void main(String[] aArgs) {
List<Double> primeReciprocals = listPrimeReciprocals(1_000_000_000);
final double euler = 0.577_215_664_901_532_861;
double sum = 0.0;
for ( double reciprocal : primeReciprocals ) {
sum += reciprocal + Math.log(1.0 - reciprocal);
}
final double meisselMertens = euler + sum;
System.out.println(String.format("%s%.9f", "The Meissel-Mertens constant = ", meisselMertens));
}
private static List<Double> listPrimeReciprocals(int aLimit) {
BitSet sieve = new BitSet(aLimit + 1);
sieve.set(2, aLimit + 1);
for ( int i = 2; i <= Math.sqrt(aLimit); i = sieve.nextSetBit(i + 1) ) {
for ( int j = i * i; j <= aLimit; j += i ) {
sieve.clear(j);
}
}
List<Double> result = new ArrayList<Double>(sieve.cardinality());
for ( int i = 2; i >= 0; i = sieve.nextSetBit(i + 1) ) {
result.add(1.0 / i);
}
return result;
}
}
</syntaxhighlight>
{{ out }}
<pre>
The Meissel-Mertens constant = 0.261497213
</pre>
=={{header|Julia}}==
Line 350 ⟶ 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 400 ⟶ 751:
{{libheader|ntheory}}
<syntaxhighlight lang="perl" line>use v5.36;
use ntheory
my $s;
say my $result = $s + .57721566490153286;</syntaxhighlight>
{{out}}
Line 556 ⟶ 907:
{{out}}
<pre>0.26149721310571383</pre>
=={{header|RPL}}==
<code>'''BPRIM?'''</code> is defined at [[Primality by trial division#RPL]]
{{works with|Halcyon Calc|4.2.8}}
≪ → n
≪ 0.5 3 n '''FOR''' j
'''IF''' j R→B '''BPRIM? THEN''' j INV + '''END''' 2 '''STEP'''
n LN LN -
≫ ≫ ''''MMCON'''' STO
10000 '''MMCON'''
{{out}}
<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>
=={{header|Wren}}==
Line 562 ⟶ 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 599 ⟶ 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
Line 654 ⟶ 1,033:
<pre>
0.261497212847642783...70842383659092665508
</pre>
=={{header|XPL0}}==
Runs in 98 seconds on Pi4.
<syntaxhighlight lang "XPL0">func IsPrime(N); \Return 'true' if N is prime
int N, D;
[if N < 2 then return false;
if (N&1) = 0 then return N = 2;
if rem(N/3) = 0 then return N = 3;
D:= 5;
while D*D <= N do
[if rem(N/D) = 0 then return false;
D:= D+2;
if rem(N/D) = 0 then return false;
D:= D+4;
];
return true;
];
def Euler = 0.57721566490153286;
real M; int P;
[M:= 0.;
for P:= 2 to 100_000_000 do
if IsPrime(P) then
M:= M + Ln(1. - 1./float(P)) + 1./float(P);
Format(1, 16);
Text(0, "MM = "); RlOut(0, Euler + M); CrLf(0);
]</syntaxhighlight>
{{out}}
<pre>
MM = 0.2614972131104150
</pre>
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">PrimeNumbers = Select[Range[100000000], PrimeQ[#] &]; (*all primes in the first 100 000 000 numbers, this takes a toll on my computer's CPU and RAM*)
MM = N[Total[Log[1 - 1/PrimeNumbers] + 1/PrimeNumbers] + EulerGamma, 10] (*Calculating it up to a precision of 10, this is correct up to 8 digits*)
AnalyticMMto305 = N[EulerGamma + Sum[MoebiusMu[n]/n Log[Zeta[n]], {n, 2, 1000}], 1000] (*Precise up to 305 digits*)
AnalyticMM = N[EulerGamma + Sum[MoebiusMu[n]/n Log[Zeta[n]], {n, 2, 10000}], 1001] (*Precise up to at least 1000 digits*)</syntaxhighlight>
{{out}}<pre>0.2614972131
0.2614972128476427837554268386086958590515666482611992061920642139249245108973682097141426314342466510516177288764860219977833903242700444245434874019723864066619495570939258171277477421198525880726627206414446423259002354310517723217392566322998031476383162375814905929038228475826597236342201597145878545286546864785868143588282690675811212349253390275898516648963722862532834592270193787663814788872688899729598789511243611212105643183719483028350694889140412496558830385604700978671612478906659270558671778566494523172201279234383779159082832574101983245266507469788929642499636163109342230499137582782542115136325204939772977424363327903550519767079810438569344262666067369673777630199473914034464437524475856770642386781282897163202273071362276688321001599406964213181471899136557333922267566346325514106895289284293459552062339943989806700523705145085201774716630458161117596603698964857884739475619061322196278674275162487770072794776455644351959504021182511495417580620912301758280597008866059
0.26149721284764278375542683860869585905156664826119920619206421392492451089736820971414263143424665105161772887648602199778339032427004442454348740197238640666194955709392581712774774211985258807266272064144464232590023543105177232173925663229980314763831623758149059290382284758265972363422015971458785446941586825460538918007031787714156680620570605257601785334398970354507934530971953511716888598019955346947142883673537117910619342522616975101911159537244599605203558051780574237201332999961769676911386909654186249097435916294862238555389898241954857937738258646582212506260380084370067541379219020626760709633535981989783010762417792511961619355361391684002933280522289185167238258837930443067100391254985761418536020400457460311825670423438456551983202200477824746954606715454777572171338072595463648319687279859427306787306509669454587505942593547068846408425666008833035029366514525328713339609172639368543886291288200447611698748441593459920236225093315001729474600911978170842383659092665508
</pre>
| |||