unit sieves;
{$mode objfpc}{$H+}
{$modeswitch advancedrecords}
{$inline on}
interface
uses
SysUtils;
const
MAX_LIMIT = High(DWord) - 131072; //to avoid overflow
type
TDWordArray = array of DWord;
{ uses bit-packed array }
function BpSieve(aLimit: DWord): TDWordArray;
{ uses bit-packed array and skips even numbers}
function OddSieve(aLimit: DWord): TDWordArray;
{ segmented sieve version }
function SegmentSieve(aLimit: DWord): TDWordArray;
implementation
type
TBitArray = record
private
FBits: array of Byte;
function GetBit(aIndex: DWord): Boolean; inline;
procedure SetBit(aIndex: DWord; aValue: Boolean); inline;
public
constructor Create(aSize: DWord);
property Bits[aIndex: DWord]: Boolean read GetBit write SetBit; default;
end;
function TBitArray.GetBit(aIndex: DWord): Boolean;
begin
Result := Boolean(FBits[aIndex shr 3] and (1 shl (aIndex and 7)));
end;
procedure TBitArray.SetBit(aIndex: DWord; aValue: Boolean);
begin
if aValue then
FBits[aIndex shr 3] := FBits[aIndex shr 3] or (1 shl (aIndex and 7))
else
FBits[aIndex shr 3] := FBits[aIndex shr 3] and not(1 shl (aIndex and 7));
end;
constructor TBitArray.Create(aSize: DWord);
begin
SetLength(FBits, aSize shr 3 + Ord(aSize and 7 <> 0));
FillChar(Pointer(FBits)^, Length(FBits), $ff);
end;
function EstimatePrimeCount(aLimit: DWord): DWord;
begin
if aLimit < 2 then
exit(0);
if aLimit <= 200 then
Result := Trunc((1.6 * aLimit)/Ln(aLimit)) + 1
else
Result := Trunc(aLimit/(Ln(aLimit) - 2)) + 1;
end;
function BpSieve(aLimit: DWord): TDWordArray;
var
IsPrime: TBitArray;
I, J, SqrtBound: DWord;
Count: Integer = 0;
begin
if aLimit > MAX_LIMIT then
raise Exception.Create('Prime limit exceeded');
if aLimit < 2 then
exit(nil);
IsPrime := TBitArray.Create(Succ(aLimit));
SqrtBound := Trunc(Sqrt(aLimit));
SetLength(Result, EstimatePrimeCount(aLimit));
for I := 2 to aLimit do
if IsPrime[I] then
begin
Result[Count] := I;
Inc(Count);
if I <= SqrtBound then
begin
J := I * I;
repeat
IsPrime[J] := False;
Inc(J, I);
until J > aLimit;
end;
end;
SetLength(Result, Count);
end;
function OddSieve(aLimit: DWord): TDWordArray;
var
IsPrime: TBitArray;
I: DWord = 3;
J, SqrtBound: DWord;
Count: Integer = 1;
begin
if aLimit > MAX_LIMIT then
raise Exception.Create('Prime limit exceeded');
if aLimit < 2 then
exit(nil);
IsPrime := TBitArray.Create(aLimit div 2);
SqrtBound := Trunc(Sqrt(aLimit));
SetLength(Result, EstimatePrimeCount(aLimit));
Result[0] := 2;
while I <= aLimit do
begin
if IsPrime[(I - 3) shr 1] then
begin
Result[Count] := I;
Inc(Count);
if I <= SqrtBound then
begin
J := I * I;
repeat
IsPrime[(J - 3) shr 1] := False;
Inc(J, I shl 1);
until J > aLimit;
end;
end;
Inc(I, 2);
end;
SetLength(Result, Count);
end;
function SegmentSieve(aLimit: DWord): TDWordArray;
const
SEG_SIZE = $8000;
var
Segment: array[0..Pred(SEG_SIZE)] of Boolean;
FirstPrimes, Primes: TDWordArray;
I, J, Prime: DWord;
Count: Integer = 0;
begin
if aLimit > MAX_LIMIT then
raise Exception.Create('Prime limit exceeded');
if aLimit <= SEG_SIZE * 2 then
exit(OddSieve(aLimit));
I := Trunc(Sqrt(aLimit)) + 1;
FirstPrimes := OddSieve(I);
SetLength(Primes, EstimatePrimeCount(aLimit) - Length(FirstPrimes));
Dec(I);
while I < aLimit do
begin
FillChar(Segment, SEG_SIZE, Byte(True));
for Prime in FirstPrimes do
begin
J := I mod Prime;
if J <> 0 then
J := Prime - J;
while J < SEG_SIZE do
begin
Segment[J] := False;
Inc(J, Prime);
end;
end;
for J := 0 to Pred(SEG_SIZE) do
if (J + I <= aLimit) and Segment[J] then
begin
Primes[Count] := J + I;
Inc(Count);
end;
Inc(I, SEG_SIZE);
end;
SetLength(FirstPrimes, Length(FirstPrimes) + Count);
Move(Primes[0], FirstPrimes[Length(FirstPrimes) - Count], Count * SizeOf(DWord));
Result := FirstPrimes;
end;
end.