Category:PrimTrial
unit primTrial
Maybe NativeUint must be typed in older versions to LongWord aka cardinal
unit PrimTrial;
// NativeUInt: LongWord 32-Bit-OS/ Uint64 64-Bit-OS
{$IFDEF FPC}
{$MODE DELPHI}
{$Smartlink ON}
{$OPTIMIZATION ON,ALL}
{$CODEALIGN proc=32}
{$ENDIF}
interface
type
ptPrimeList = array of NativeUint;
procedure InitPrime;
function actPrime :NativeUint;
function isPrime(pr: NativeUint):boolean;
function isAlmostPrime(n: NativeUint;cnt: NativeUint): boolean;
function isSemiprime(n: NativeUint): boolean;
function SmallFactor(pr: NativeUint):NativeUint;
//next prime
function NextPrime: NativeUint;
//next possible prime of number wheel
function NextPosPrim: NativeUint;
//next prime greater equal limit
function PrimeGELimit(Limit:NativeUint):NativeUint;
function PrimeRange(LowLmt,UpLmt:NativeUint): ptPrimeList;
implementation
uses
sysutils;
const
cntsmallPrimes = 6;
smallPrimes : array[0..cntsmallPrimes-1] of NativeUint = (2,3,5,7,11,13);
wheelSize = (2-1)*(3-1)*(5-1)*(7-1)*(11-1)*(13-1);
wheelCircumfence = 2*3*5*7*11*13;
var
deltaWheel : array[0..wheelSize-1] of byte;
WheelIdx : nativeUint;
p,pw : nativeUint;
procedure InitPrime;
//initialies wheel and prime to startposition
Begin
p := 2;
pw := 1;
WheelIdx := 0;
end;
function actPrime :NativeUint;
Begin
result := p;
end;
procedure InitWheel;
//search for numbers that are no multiples of smallprimes
//saving only the distance, to keep size small
var
p0,p1,i,d,res : NativeUint;
Begin
p0 := 1;d := 0;p1 := p0;
repeat
Repeat
p1 := p1+2;// only odd
i := 1;
repeat
res := p1 mod smallPrimes[i];
inc(i)
until (res =0)OR(i >= cntSmallPrimes);
if res <> 0 then
Begin
deltaWheel[d] := p1-p0;
inc(d);
break;
end;
until false;
p0 := p1;
until d >= wheelSize;
end;
function biggerFactor(p: NativeUint):NativeUint;
//trial division by wheel numbers
//reduces count of divisions from 1/2 = 0.5( only odd numbers )
//to 5760/30030 = 0.1918
var
sp : NativeUint;
d : NativeUint;
r : NativeUint;
Begin
sp := 1;d := 0;
repeat
sp := sp+deltaWheel[d];
r := p mod sp;
d := d+1;
//IF d = WheelSize then d := 0;
d := d AND NativeUint(-ord(d<>WheelSize));
IF r = 0 then
BREAK;
until p < sp*sp;
IF r = 0 then
result := sp
else
result := p;
end;
function SmallFactor(pr: NativeUint):NativeUint;
//checking numbers omitted by biggerFactor
const
// Bit Mask for small primes
Check2_61 :UInt64 = 1 shl 2+1 shl 3+1 shl 5+1 shl 7+1 shl 11+1 shl 13+1 shl 17+
1 shl 19+1 shl 23+1 shl 29+1 shl 31+1 shl 37+1 shl 41+1 shl 47+
1 shl 53+1 shl 59+1 shl 61;
var
k : NativeUint;
Begin
result := pr;
if pr < 64 then
if (Uint64(1) shl pr) AND Check2_61 <> 0 then
EXIT;
IF NOT(ODD(pr))then
Begin
result := 2;
EXIT;
end;
For k := 1 to cntSmallPrimes-1 do
Begin
IF pr Mod smallPrimes[k] = 0 then
Begin
result := smallPrimes[k];
EXIT
end;
end;
k := smallPrimes[cntsmallPrimes-1];
IF pr>k*k then
result := biggerFactor(pr);
end;
function isPrime(pr: NativeUint):boolean;
Begin
IF pr > 1 then
isPrime := smallFactor(pr) = pr
else
isPrime := false;
end;
function isAlmostPrime(n: NativeUint;cnt: NativeUint): boolean;
var
fac1,c : NativeUint;
begin
c := 0;
repeat
fac1 := SmallFactor(n);
n := n div fac1;
inc(c);
until (n = 1) OR (c > cnt);
isAlmostPrime := (n = 1) AND (c = cnt);
end;
function isSemiprime(n: NativeUint): boolean;
begin
result := isAlmostPrime(n,2);
end;
function NextPosPrim: NativeUint;
var
WI : NativeUint;
Begin
result := pw+deltaWheel[WheelIdx];
WI := (WheelIdx+1);
WheelIdx := WI AND NativeUint(-ORD(WI<>WheelSize));
pw := result;
end;
function NextPrime: NativeUint;
Begin
IF p >= smallPrimes[High(smallPrimes)]then
Begin
repeat
until isPrime(NextPosPrim);
result := pw;
p := result;
end
else
Begin
result := 0;
while p >= smallPrimes[result] do
inc(result);
result := smallPrimes[result];
p:= result;
end;
end;
function PrimeGELimit(Limit:NativeUint):NativeUint;
//prime greater or equal limit
Begin
IF Limit > wheelCircumfence then
Begin
WheelIdx:= wheelSize-1;
result := (Limit DIV wheelCircumfence)*wheelCircumfence-1;
pw := result;
//the easy way, no prime test
while pw <= Limit do
NextPosPrim;
result := pw;
p := result;
if Not(isPrime(result)) then
result := NextPrime;
end
else
Begin
InitPrime;
repeat
until (NextPosPrim >= limit) AND isPrime(pw);
result := pw;
p := result;
end;
end;
function PrimeRange(LowLmt,UpLmt:NativeUint): ptPrimeList;
var
i,newP : NativeUint;
Begin
IF LowLmt>UpLmt then
Begin
setlength(result,0);
EXIT;
end;
i := 0;
setlength(result,100);
newP := PrimeGELimit(LowLmt);
while newP<= UpLmt do
Begin
result[i]:= newP;
inc(i);
IF i>High(result) then
setlength(result,i*2);
newP := NextPrime;
end;
setlength(result,i);
end;
//initialization
Begin
InitWheel;
InitPrime;
end.
Pages in category "PrimTrial"
The following 11 pages are in this category, out of 11 total.