unit Unit1;
{$mode objfpc}{$H+}
interface
uses
Classes, math, SysUtils, FileUtil, Forms, Controls, Graphics,
StdCtrls;
type
{ TForm1 }
TForm1 = class(TForm)
btnCalculaSieve: TButton;
Edit1: TEdit;
Memo1: TMemo;
Memo2: TMemo;
procedure btnCalculaSieveClick(Sender: TObject);
private
procedure segmented_sieve(const limit: int64);
public
end;
{ std_vector }
generic std_vector<T> = class
type TarrT = array of T;
strict private
FItems: TarrT;
FCount: T;
function GetItem(aIndex: T): T;
procedure SetItem(const aIndex: T;const AValue: T);
public
procedure Add(const Value: T);
property Items[aIndex: T]: T read GetItem write SetItem; default;
property Count: T read FCount;
end;
var
Form1: TForm1;
implementation
{$R *.frm}
const
/// Set your CPU's L1 data cache size (in bytes) here
L1D_CACHE_SIZE : integer = 32768;
init_count: array [boolean] of byte = (1, 0); // for emulating (limit < 2) ? 0 : 1;
{ std_vector implementation }
function std_vector.GetItem(aIndex: T): T;
begin
result := FItems[aIndex];
end;
procedure std_vector.SetItem(const aIndex: T;const AValue: T);
begin
if FItems[aIndex] <> aValue then
FItems[aIndex] := aValue;
end;
procedure std_vector.Add(const Value: T);
begin
inc(FCount);
SetLength(FItems, FCount);
FItems[FCount-1] := Value;
end;
{ TForm1 }
procedure TForm1.segmented_sieve(const limit: int64);
var
sqrt,
segment_size,
i, j, k,
S, N, _high:int64;
count: int64 = 0;
_low: int64 = 0;
is_prime, /// std::vector<char> is_prime
sieve: array of byte; /// std::vector<char> sieve
primes, // can be TFPList or use TList instead
_next: specialize std_vector<QWord>;
begin
sqrt := Trunc(System.sqrt(limit));
segment_size := Max(sqrt, L1D_CACHE_SIZE);
inc(count, init_count[limit < 2]);
S := 3;
N := 3;
// generates small primes <= sqrt
SetLength(is_prime, sqrt);
//FillByte(is_prime[0], length(is_prime), 1); // Little improvement 3 next lines
FillWord(is_prime[0], length(is_prime) div 2, $100);
FillWord(is_prime[0], 2, $101);
i := 3;
while (i*i <= sqrt) do begin
if boolean(is_prime[i]) then begin
j := i*i;
while (j <= sqrt) do begin
is_prime[j] := 0;
inc(j,i);
end;
end;
inc(i);
end;
// vector used to sieve
SetLength(sieve, segment_size);
primes := specialize std_vector<longWord>.Create;
_next := specialize std_vector<longWord>.Create;
while (_low <= limit) do begin
FillByte(sieve[0], length(sieve), 1);
// actual segment = interval [_low, _high]
_high := math.min(_low + segment_size - 1, limit);
// adds new primes in sieve <= sqrt(_high)
while (S * S <= _high) do begin
if boolean(is_prime[S]) then begin
primes.Add(S);
_next.Add(integer(S * S - _low));
end;
inc(S, 2);
end;
// sieves the actual segment
i := 0;
while (i < primes.Count) do begin
j := _next[i];
k := primes[i] * 2;
while (j < segment_size) do begin
sieve[j] := 0;
inc(j,k);//j+=k
end;
_next[i] := j - segment_size;
inc(i);
end;
while (N <= _high) do begin
if boolean(sieve[N - _low]) then
inc(count);
inc(N, 2);
end;
inc(_low, segment_size);//_low += segment_size;
end;
Form1.Caption := intToStr(count);
primes.Free;
_next.Free;
end;
procedure TForm1.btnCalculaSieveClick(Sender: TObject);
var
lim: Int64;
DT: QWord;
begin
DT := GetTickCount64;
Caption := 'Calculating...';
lim := StrToInt(Edit1.Text);
segmented_sieve(lim);
DT := GetTickCount64 - DT;
Caption := Caption + ' : [' + IntToStr(DT) + '] millisecs';
end;
end.