Legendre prime counting function: Difference between revisions
Added solution for Pascal
GordonBGood (talk | contribs) (→{{header|Erlang}}: add Elm contribution above...) |
(Added solution for Pascal) |
||
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).
=={{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}}==
| |||