Legendre prime counting function: Difference between revisions
Content added Content deleted
GordonBGood (talk | contribs) (→{{header|Erlang}}: add Elm contribution above...) |
(Added solution for Pascal) |
||
Line 2,916: | Line 2,916: | ||
As counting range is increased above these trivial ranges, a better algorithm such as that of LMO will likely perform much better than this and use much less memory (to proportional to the cube root of the range), although at the cost of greater code complexity and more lines of code (about 500 LOC). The trade-off between algorithms of the Legendre type and those of the Meissel type is that the latter decreases the time complexity but at the cost of more time spent sieving; modern follow-on work to LMO produces algorithms which are able to tune the balance between Legendre-type and Meissel-type algorithm in trading off the execution time costs between the different parts of the given algorithm for the best in performance, with the latest in Kim Walisch's tuning of Xavier Gourdon's work being about five to ten times faster than LMO (plus includes multi-threading abilities). |
As counting range is increased above these trivial ranges, a better algorithm such as that of LMO will likely perform much better than this and use much less memory (to proportional to the cube root of the range), although at the cost of greater code complexity and more lines of code (about 500 LOC). The trade-off between algorithms of the Legendre type and those of the Meissel type is that the latter decreases the time complexity but at the cost of more time spent sieving; modern follow-on work to LMO produces algorithms which are able to tune the balance between Legendre-type and Meissel-type algorithm in trading off the execution time costs between the different parts of the given algorithm for the best in performance, with the latest in Kim Walisch's tuning of Xavier Gourdon's work being about five to ten times faster than LMO (plus includes multi-threading abilities). |
||
=={{header|Pascal}}== |
|||
<syntaxhighlight lang="pascal"> |
|||
// Rosetta Code task "Legendre prime counting function". |
|||
// Solution for Free Pascal (Lazarus) or Delphi. |
|||
program LegendrePrimeCount; |
|||
{$IFDEF FPC} // Free Pascal |
|||
{$MODE Delphi} |
|||
{$ELSE} // Delphi |
|||
{$APPTYPE CONSOLE} |
|||
{$ENDIF} |
|||
// Optimization needs to be on if the program is to finish in a reasonable |
|||
// length of time. See "Comments on Task" on the Rosetta Code website. |
|||
{$DEFINE SIMPLE_OPTIMIZATION} |
|||
uses SysUtils, Types; |
|||
{------------------------------------------------------------------- |
|||
Function to return an array of primes up to the passed-in limit. |
|||
Uses a straightforward Eratosthenes sieve. |
|||
TIntegerDynArray must be 0-based. To be compatible with Rosetta Code, |
|||
the first prime 2 is at result[1], and result[0] is not used. |
|||
} |
|||
function FindPrimes( limit : integer) : Types.TIntegerDynArray; |
|||
var |
|||
deleted : array of boolean; |
|||
j, k, p, resultSize : integer; |
|||
begin |
|||
if (limit < 2) then begin |
|||
SetLength( result, 1); |
|||
exit; |
|||
end; |
|||
SetLength( deleted, limit + 1); // 0..limit |
|||
deleted[0] := true; |
|||
for j := 1 to limit do deleted[j] := false; |
|||
p := 2; |
|||
while (p*p <= limit) do begin |
|||
j := 2*p; |
|||
while (j <= limit) do begin |
|||
deleted[j] := true; |
|||
inc( j, p); |
|||
end; |
|||
repeat inc(p) |
|||
until (p > limit) or (not deleted[p]); |
|||
end; |
|||
resultSize := 0; |
|||
for j := 0 to limit do |
|||
if not deleted[j] then inc( resultSize); |
|||
SetLength( result, resultSize); |
|||
k := 0; |
|||
for j := 0 to limit do begin |
|||
if not deleted[j] then begin |
|||
result[k] := j; |
|||
inc(k); |
|||
end; |
|||
end; |
|||
end; |
|||
{----------------------------------------------------------------------------- |
|||
Function to count primes up to the passed-in limit, by Legendre's method. |
|||
Iterative, using a stack. Each item in the stack is a term phi(x,a) along |
|||
with a sign. If the top item on the stack can be evaluated easily, it is |
|||
popped off and its value is added to the result. Else the top item is |
|||
replaced by two items according to the formual in the task description. |
|||
} |
|||
function CountPrimes( n : integer) : integer; |
|||
type |
|||
TPhiTerm = record |
|||
IsNeg : boolean; |
|||
x : integer; |
|||
a : integer; |
|||
end; |
|||
const |
|||
STACK_SIZE = 100; // 10 is enough for n = 10^9 |
|||
var |
|||
primes : Types.TIntegerDynArray; |
|||
nrPrimes : integer; |
|||
stack : array [0..STACK_SIZE - 1] of TPhiTerm; |
|||
sp : integer; // stack pointer, points to first free entry |
|||
tos : TPhiTerm; // top of stack |
|||
begin |
|||
primes := FindPrimes( Trunc( Sqrt( n + 0.5))); |
|||
nrPrimes := Length( primes) - 1; // primes[0] is not used |
|||
result := nrPrimes - 1; // initialize total |
|||
// Push initial entry onto stack |
|||
with stack[0] do begin |
|||
IsNeg := false; |
|||
x := n; |
|||
a := nrPrimes; |
|||
end; |
|||
sp := 1; |
|||
while sp > 0 do begin |
|||
tos := stack[sp - 1]; |
|||
{$IFDEF SIMPLE_OPTIMIZATION} |
|||
// Using optimization described in "Comments on Task" |
|||
if tos.x = 0 then begin // top of stack = 0 |
|||
dec(sp); // pop top of stack, no change to result |
|||
end |
|||
else if (tos.a > 0) and (tos.x < primes[tos.a]) then begin // top of stack = 1 |
|||
dec( sp); // pop top of stack, update result |
|||
if tos.IsNeg then dec( result) |
|||
else inc( result); |
|||
end |
|||
else if tos.a = 0 then begin |
|||
{$ELSE} |
|||
// Using only the task description, i.e. only phi(x,0) = x |
|||
if tos.a = 0 then begin |
|||
{$ENDIF} |
|||
dec( sp); // pop top of stack, update result |
|||
if tos.IsNeg then dec( result, tos.x) |
|||
else inc( result, tos.x); |
|||
end |
|||
else begin |
|||
// Replace top of stack by two items as in the task description, |
|||
// namely phi(x, a - 1) and -phi(x div primes[a], a - 1) |
|||
if (sp >= STACK_SIZE) then |
|||
raise SysUtils.Exception.Create( 'Legendre phi stack overflow'); |
|||
with stack[sp - 1] do begin |
|||
IsNeg := tos.isNeg; |
|||
x := tos.x; |
|||
a := tos.a - 1; |
|||
end; |
|||
with stack[sp] do begin |
|||
IsNeg := not tos.IsNeg; |
|||
x := tos.x div primes[tos.a]; |
|||
a := tos.a - 1; |
|||
end; |
|||
inc(sp); |
|||
end; |
|||
end; |
|||
end; |
|||
{----------------------------------------------------------- |
|||
Main routine |
|||
} |
|||
var |
|||
power, limit, count : integer; |
|||
begin |
|||
WriteLn( 'Limit Count'); |
|||
limit := 1; |
|||
for power := 0 to 9 do begin |
|||
if power > 0 then limit := 10*limit; |
|||
count := CountPrimes( limit); |
|||
WriteLn( SysUtils.Format( '10^%d %10d', [power, count])) |
|||
end; |
|||
end. |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Limit Count |
|||
10^0 0 |
|||
10^1 4 |
|||
10^2 25 |
|||
10^3 168 |
|||
10^4 1229 |
|||
10^5 9592 |
|||
10^6 78498 |
|||
10^7 664579 |
|||
10^8 5761455 |
|||
10^9 50847534 |
|||
</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |