Category:PrimTrial

UNIT for primTrial

Works with: Free Pascal
Works with: Delphi

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=16}
{\$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.
```

UNIT for primsieve

segmented sieve of erathostenes.Not that bad.

```unit primsieve;
{\$IFDEF FPC}
{\$IFEND}
{segmented sieve of Erathostenes using only odd numbers}
{using presieved sieve of small primes, to reduce the most time consuming}
interface
procedure InitPrime;
procedure NextSieve;
function LastPrimeInBlock:Uint64;
function SieveStart:Uint64;
function SieveSize :LongInt;
function Nextprime: Uint64;
function StartCount :Uint64;
function TotalCount :Uint64;
function PosOfPrime: Uint64;

implementation
uses
sysutils;
const
smlPrimes :array [0..10] of Byte = (2,3,5,7,11,13,17,19,23,29,31);
cSieveSize = 16384 * 4; //<= High(Word)+1 // Level I Data Cache
type
tSievePrim = record
svdeltaPrime:word;//diff between actual and new prime
svSivOfs:word;    //Offset in sieve
svSivNum:LongWord;//1 shl (1+16+32) = 5.6e14
end;
tpSievePrim = ^tSievePrim;

var
//sieved with primes 3..maxPreSievePrime.here about 255255 Byte
{\$ALIGN 32}
preSieve :array[0..3*5*7*11*13*17-1] of Byte;//must be > cSieveSize
{\$ALIGN 32}
Sieve :array[0..cSieveSize-1] of Byte;
{\$ALIGN 32}
//prime = FoundPrimesOffset + 2*FoundPrimes[0..FoundPrimesCnt]
FoundPrimes : array[0..cSieveSize] of word;
{\$ALIGN 32}
sievePrimes : array[0..78498] of tSievePrim;// 1e6^2 ->1e12
//  sievePrimes : array[0..664579] of tSievePrim;// maximum 1e14
FoundPrimesOffset : Uint64;
LastPrimeInSieve,
FoundPrimesCnt,
FoundPrimesIdx,
FoundPrimesTotal,
SieveNum,
SieveMaxIdx,
preSieveOffset,
LastInsertedSievePrime :NativeUInt;

procedure CopyPreSieveInSieve; forward;
procedure CollectPrimes; forward;
procedure sieveOneSieve; forward;
procedure Init0Sieve; forward;
procedure SieveOneBlock; forward;

//****************************************
procedure preSieveInit;
var
i,pr,j,umf : NativeInt;
Begin
fillchar(preSieve[0],SizeOf(preSieve),#1);
i := 1;
pr := 3;// starts with pr = 3
umf := 1;
repeat
IF preSieve[i] =1 then
Begin
pr := 2*i+1;
j := i;
repeat
preSieve[j] := 0;
inc(j,pr);
until j> High(preSieve);
umf := umf*pr;
end;
inc(i);
until (pr = maxPreSievePrime)OR(umf>High(preSieve)) ;
preSieveOffset := 0;
end;

function InsertSievePrimes(PrimPos:NativeInt):NativeInt;
var
delta    :NativeInt;
i,pr,loLmt : NativeUInt;
begin
i := 0;
//ignore first primes already sieved with
if SieveNum = 0 then
pr :=0;
loLmt := Uint64(SieveNum)*(2*cSieveSize);
delta := loLmt-LastInsertedSievePrime;
with sievePrimes[PrimPos] do
Begin
pr := FoundPrimes[i]*2+1;
svdeltaPrime := pr+delta;
delta := pr;
end;

inc(PrimPos);
for i := i+1 to FoundPrimesCnt-1 do
Begin
IF PrimPos > High(sievePrimes) then
BREAK;
with sievePrimes[PrimPos] do
Begin
pr := FoundPrimes[i]*2+1;
svdeltaPrime := (pr-delta);
delta := pr;
end;
inc(PrimPos);
end;
LastInsertedSievePrime := loLmt+pr;
result := PrimPos;
end;

procedure CalcSievePrimOfs(lmt:NativeUint);
//lmt High(sievePrimes)
var
i,pr : NativeUInt;
sq : Uint64;
begin
pr := 0;
i := 0;
repeat
with sievePrimes[i] do
Begin
pr := pr+svdeltaPrime;
IF sqr(pr)  < (cSieveSize*2) then
Begin
svSivNum := 0;
svSivOfs := (pr*pr-1) DIV 2;
end
else
Begin
SieveMaxIdx := i;
pr := pr-svdeltaPrime;
BREAK;
end;
end;
inc(i);
until i > lmt;

for i := i to lmt do
begin
with sievePrimes[i] do
Begin
pr := pr+svdeltaPrime;
sq := sqr(pr);
svSivNum := sq DIV (2*cSieveSize);
svSivOfs := ( (sq - Uint64(svSivNum)*(2*cSieveSize))-1)DIV 2;
end;
end;
end;

procedure sievePrimesInit;
var
i,j,pr,PrimPos:NativeInt;
Begin
LastInsertedSievePrime := 0;
preSieveOffset := 0;
SieveNum :=0;
CopyPreSieveInSieve;
//normal sieving of first block sieve
repeat
while Sieve[i] = 0 do
inc(i);
pr := 2*i+1;
inc(i);
j := ((pr*pr)-1) DIV 2;
if j > High(Sieve) then
BREAK;
repeat
Sieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
until false;

CollectPrimes;
PrimPos := InsertSievePrimes(0);
//correct for SieveNum = 0
CalcSievePrimOfs(PrimPos);
Init0Sieve;
sieveOneBlock;
//now start collect with SieveNum = 1
IF PrimPos < High(sievePrimes) then
repeat
sieveOneBlock;
CollectPrimes;
dec(SieveNum);
PrimPos := InsertSievePrimes(PrimPos);
inc(SieveNum);
until PrimPos > High(sievePrimes);
Init0Sieve;
end;

procedure Init0Sieve;
begin
FoundPrimesTotal :=0;
preSieveOffset := 0;
SieveNum :=0;
CalcSievePrimOfs(High(sievePrimes));
end;

procedure CopyPreSieveInSieve;
var
lmt : NativeInt;
Begin
lmt := preSieveOffset+cSieveSize;
lmt := lmt-(High(preSieve)+1);
IF lmt<= 0 then
begin
Move(preSieve[preSieveOffset],Sieve[0],cSieveSize);
if lmt <> 0 then
inc(preSieveOffset,cSieveSize)
else
preSieveOffset := 0;
end
else
begin
Move(preSieve[preSieveOffset],Sieve[0],cSieveSize-lmt);
Move(preSieve[0],Sieve[cSieveSize-lmt],lmt);
preSieveOffset := lmt
end;
end;

procedure sieveOneSieve;
var
sp:tpSievePrim;
pSieve :pByte;
i,j,pr,sn,dSievNum :NativeUint;
Begin
pr := 0;
sn := sieveNum;
sp := @sievePrimes[0];
pSieve := @Sieve[0];
For i := SieveMaxIdx downto 0 do
with sp^ do
begin
pr := pr+svdeltaPrime;
IF svSivNum = sn then
Begin
j := svSivOfs;
repeat
pSieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
dSievNum := j DIV cSieveSize;
svSivOfs := j-dSievNum*cSieveSize;
svSivNum := sn+dSievNum;
//        svSivNum := svSivNum+dSievNum;
end;
inc(sp);
end;
i := SieveMaxIdx+1;
repeat
if i > High(SievePrimes) then
BREAK;
with sp^ do
begin
if svSivNum > sn then
Begin
SieveMaxIdx := I-1;
Break;
end;
pr := pr+svdeltaPrime;
j := svSivOfs;
repeat
Sieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
dSievNum := j DIV cSieveSize;
svSivOfs := j-dSievNum*cSieveSize;
svSivNum := sn+dSievNum;
end;
inc(i);
inc(sp);
until false;
end;

procedure CollectPrimes;
//extract primes to FoundPrimes
var
pSieve : pbyte;
pFound : pWord;
i,idx : NativeUint;
Begin
FoundPrimesOffset := SieveNum*(2*cSieveSize);
FoundPrimesIdx := 0;
pFound :=@FoundPrimes[0];
i := 0;
idx := 0;
IF SieveNum = 0 then
//include small primes used to pre-sieve
Begin
repeat
pFound[idx]:= (smlPrimes[idx]-1) DIV 2;
inc(idx);
until smlPrimes[idx]>maxPreSievePrime;
i := (smlPrimes[idx] -1) DIV 2;
end;
//grabbing the primes without if then -> reduces time extremly
//primes are born to let branch-prediction fail.
pSieve:= @Sieve[Low(Sieve)];
repeat
//store every value until a prime aka 1 is found
pFound[idx]:= i;
inc(idx,pSieve[i]);
inc(i);
until i>High(Sieve);
FoundPrimesCnt:= idx;
inc(FoundPrimesTotal,Idx);
LastPrimeInSieve := FoundPrimesOffset+pSieve[idx]*2+1;
end;

procedure SieveOneBlock;inline;
begin
CopyPreSieveInSieve;
sieveOneSieve;
CollectPrimes;
inc(SieveNum);
end;

procedure NextSieve;inline;
Begin
SieveOneBlock;
end;
function LastPrimeInBlock:Uint64;inline;
begin
result := LastPrimeInSieve;
end;
function Nextprime:Uint64;
Begin
result := FoundPrimes[FoundPrimesIdx]*2+1+FoundPrimesOffset;
if (FoundPrimesIdx=0) AND (sievenum = 1) then
inc(result);
inc(FoundPrimesIdx);
If FoundPrimesIdx>= FoundPrimesCnt then
SieveOneBlock;
end;

function PosOfPrime: Uint64;inline;
Begin
result := FoundPrimesTotal-FoundPrimesCnt+FoundPrimesIdx;
end;

function StartCount : Uint64 ;inline;
begin
result := FoundPrimesTotal-FoundPrimesCnt;
end;

function TotalCount :Uint64;inline;
begin
result := FoundPrimesTotal;
end;

function SieveSize :LongInt;inline;
Begin
result := 2*cSieveSize;
end;

function SieveStart:Uint64;inline;
Begin
result := (SieveNum-1)*2*cSieveSize;
end;

procedure InitPrime;inline;
Begin
Init0Sieve;
SieveOneBlock;
end;

begin
preSieveInit;
sievePrimesInit;
InitPrime;
end.

the test program

program test;
{\$IFDEF FPC}
{\$MODE objFPC}{\$Optimization ON,ALL}
{\$IFEND}

uses
primsieve;
var
cnt,p,lmt : Uint64;
Begin
lmt := 1000*1000*1000;
p := 0;
while TotalCount < lmt do
Begin
NextSieve;
inc(p);
If p AND (4096-1) = 0 then
write(p:8,TotalCount:15,#13);
end;
cnt := StartCount;
repeat
p := NextPrime;
inc(cnt);
until cnt >= lmt;
writeln(cnt:14,p:14);
end.
```

UNIT for prime decomposition

```unit
primeDecomp;

// gets factors of consecutive integers fast
// limited to 1.2e11
{\$IFDEF FPC}
{\$ENDIF}
{\$IFDEF WINDOWS}
{\$APPTYPE CONSOLE}
{\$ENDIF}
//prime decomposition
interface
const
//HCN(86) > 1.2E11 = 128,501,493,120     count of divs = 4096   7 3 1 1 1 1 1 1 1
HCN_DivCnt  = 4096;
type
tItem     = Uint64;
tDivisors = array [0..HCN_DivCnt] of tItem;
tpDivisor = pUint64;
const
//used odd size for test only
SizePrDeFe = 576*64;//*SizeOf(tprimeFac) <= 32kb level I or 2 Mb ~ level 2 cache
type
tdigits = array [0..31] of Uint32;
//the first number with 11 different prime factors =
//2*3*5*7*11*13*17*19*23*29*31 = 2E11
//56 byte
tprimeFac = packed record
pfSumOfDivs,
pfRemain  : Uint64;
pfDivCnt  : Uint32;
pfMaxIdx  : Uint32;
pfpotPrimIdx : array[0..9] of word;
pfpotMax  : array[0..11] of byte;
end;
tpPrimeFac = ^tprimeFac;

tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac;
tPrimes = array[0..65535] of Uint32;

procedure InitSmallPrimes;
function  Init_Sieve(n:NativeUint):boolean;
function  GetNextPrimeDecomp:tpPrimeFac;
function  smplPrimeDecomp(n:Uint64):tprimeFac;
procedure GetDivisors(pD:tpPrimeFac;var Divs:tDivisors);
procedure AllFacsOut(var Divs:tdivisors;proper:boolean=true);
function  OutPots(pD:tpPrimeFac;n:NativeInt):Ansistring;

implementation
var
{\$ALIGN 8}
SmallPrimes: tPrimes;
{\$ALIGN 32}
PrimeDecompField :tPrimeDecompField;
pdfIDX,pdfOfs: NativeInt;
SmallPrimesDone : boolean = false;

procedure InitSmallPrimes;
//get primes. #0..65535.Sieving only odd numbers
const
MAXLIMIT = (821641-1) shr 1;
var
pr : array[0..MAXLIMIT] of byte;
p,j,d,flipflop :NativeUInt;
Begin
if SmallPrimesDone then
EXIT;
SmallPrimes[0] := 2;
fillchar(pr[0],SizeOf(pr),#0);
p := 0;
repeat
repeat
p +=1
until pr[p]= 0;
j := (p+1)*p*2;
if j>MAXLIMIT then
BREAK;
d := 2*p+1;
repeat
pr[j] := 1;
j += d;
until j>MAXLIMIT;
until false;

SmallPrimes[1] := 3;
SmallPrimes[2] := 5;
j := 3;
d := 7;
flipflop := (2+1)-1;//7+2*2,11+2*1,13,17,19,23
p := 3;
repeat
if pr[p] = 0 then
begin
SmallPrimes[j] := d;
inc(j);
end;
d += 2*flipflop;
p+=flipflop;
flipflop := 3-flipflop;
until (p > MAXLIMIT) OR (j>High(SmallPrimes));
SmallPrimesDone := true;
end;

function OutPots(pD:tpPrimeFac;n:NativeInt):Ansistring;
var
s: String[31];
chk,p,i: NativeInt;
Begin
str(n,s);
result := s+' :';
with pd^ do
begin
str(pfDivCnt:3,s);
result += s+' : ';
chk := 1;
For n := 0 to pfMaxIdx-1 do
Begin
if n>0 then
result += '*';
p := SmallPrimes[pfpotPrimIdx[n]];
chk *= p;
str(p,s);
result += s;
i := pfpotMax[n];
if i >1 then
Begin
str(pfpotMax[n],s);
result += '^'+s;
repeat
chk *= p;
dec(i);
until i <= 1;
end;

end;
p := pfRemain;
If p >1 then
Begin
str(p,s);
chk *= p;
result += '*'+s;
end;
str(chk,s);
result += '_chk_'+s+'<';
str(pfSumOfDivs,s);
result += '_SoD_'+s+'<';
end;
end;

function smplPrimeDecomp(n:Uint64):tprimeFac;
var
pr,i,pot,fac,q :NativeUInt;
Begin
with result do
Begin
pfDivCnt := 1;
pfSumOfDivs := 1;
pfRemain := n;
pfMaxIdx := 0;
pfpotPrimIdx[0] := 1;
pfpotMax[0] := 0;

i := 0;
while i < High(SmallPrimes) do
begin
pr := SmallPrimes[i];
q := n DIV pr;
//if n < pr*pr
if pr > q then
BREAK;
if n = pr*q then
Begin
pfpotPrimIdx[pfMaxIdx] := i;
pot := 0;
fac := pr;
repeat
n := q;
q := n div pr;
pot+=1;
fac *= pr;
until n <> pr*q;
pfpotMax[pfMaxIdx] := pot;
pfDivCnt *= pot+1;
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfMaxIdx);
end;
inc(i);
end;
pfRemain := n;
if n > 1 then
Begin
pfDivCnt *= 2;
pfSumOfDivs *= n+1
end;
end;
end;

function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt;
//n must be multiple of base aka n mod base must be 0
var
q,r: Uint64;
i : NativeInt;
Begin
fillchar(dgt,SizeOf(dgt),#0);
i := 0;
n := n div base;
result := 0;
repeat
r := n;
q := n div base;
r  -= q*base;
n := q;
dgt[i] := r;
inc(i);
until (q = 0);
//searching lowest pot in base
result := 0;
while (result<i) AND (dgt[result] = 0) do
inc(result);
inc(result);
end;

function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt;
var
q :NativeInt;
Begin
result := 0;
q := dgt[result]+1;
if q = base then
repeat
dgt[result] := 0;
inc(result);
q := dgt[result]+1;
until q <> base;
dgt[result] := q;
result +=1;
end;

function SieveOneSieve(var pdf:tPrimeDecompField):boolean;
var
dgt:tDigits;
i,j,k,pr,fac,n,MaxP : Uint64;
begin
n := pdfOfs;
if n+SizePrDeFe >= sqr(SmallPrimes[High(SmallPrimes)]) then
EXIT(FALSE);
//init
for i := 0 to SizePrDeFe-1 do
begin
with pdf[i] do
Begin
pfDivCnt := 1;
pfSumOfDivs := 1;
pfRemain := n+i;
pfMaxIdx := 0;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := 0;
end;
end;
//first factor 2. Make n+i even
i := (pdfIdx+n) AND 1;
IF (n = 0) AND (pdfIdx<2)  then
i := 2;

repeat
with pdf[i] do
begin
j := BsfQWord(n+i);
pfMaxIdx := 1;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := j;
pfRemain := (n+i) shr j;
pfSumOfDivs := (Uint64(1) shl (j+1))-1;
pfDivCnt := j+1;
end;
i += 2;
until i >=SizePrDeFe;
//i now index in SmallPrimes
i := 0;
maxP := trunc(sqrt(n+SizePrDeFe))+1;
repeat
//search next prime that is in bounds of sieve
if n = 0 then
begin
repeat
inc(i);
pr := SmallPrimes[i];
k := pr-n MOD pr;
if k < SizePrDeFe then
break;
until pr > MaxP;
end
else
begin
repeat
inc(i);
pr := SmallPrimes[i];
k := pr-n MOD pr;
if (k = pr) AND (n>0) then
k:= 0;
if k < SizePrDeFe then
break;
until pr > MaxP;
end;

//no need to use higher primes
if pr*pr > n+SizePrDeFe then
BREAK;

//j is power of prime
j := CnvtoBASE(dgt,n+k,pr);
repeat
with pdf[k] do
Begin
pfpotPrimIdx[pfMaxIdx] := i;
pfpotMax[pfMaxIdx] := j;
pfDivCnt *= j+1;
fac := pr;
repeat
pfRemain := pfRemain DIV pr;
dec(j);
fac *= pr;
until j<= 0;
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfMaxIdx);
k += pr;
j := IncByBaseInBase(dgt,pr);
end;
until k >= SizePrDeFe;
until false;

//correct sum of & count of divisors
for i := 0 to High(pdf) do
Begin
with pdf[i] do
begin
j := pfRemain;
if j <> 1 then
begin
pfSumOFDivs *= (j+1);
pfDivCnt *=2;
end;
end;
end;
result := true;
end;

function NextSieve:boolean;
begin
dec(pdfIDX,SizePrDeFe);
inc(pdfOfs,SizePrDeFe);
result := SieveOneSieve(PrimeDecompField);
end;

function GetNextPrimeDecomp:tpPrimeFac;
begin
if pdfIDX >= SizePrDeFe then
if Not(NextSieve) then
EXIT(NIL);
result := @PrimeDecompField[pdfIDX];
inc(pdfIDX);
end;

function Init_Sieve(n:NativeUint):boolean;
//Init Sieve pdfIdx,pdfOfs are global
begin
pdfIdx := n MOD SizePrDeFe;
pdfOfs := n-pdfIdx;
result := SieveOneSieve(PrimeDecompField);
end;

procedure InsertSort(pDiv:tpDivisor; Left, Right : NativeInt );
var
I, J: NativeInt;
Pivot : tItem;
begin
for i:= 1 + Left to Right do
begin
Pivot:= pDiv[i];
j:= i - 1;
while (j >= Left) and (pDiv[j] > Pivot) do
begin
pDiv[j+1]:=pDiv[j];
Dec(j);
end;
pDiv[j+1]:= pivot;
end;
end;

procedure GetDivisors(pD:tpPrimeFac;var Divs:tDivisors);
var
pDivs : tpDivisor;
pPot : UInt64;
i,len,j,l,p,k: Int32;
Begin
pDivs := @Divs[0];
pDivs[0] := 1;
len := 1;
l   := 1;
with pD^ do
Begin
For i := 0 to pfMaxIdx-1 do
begin
//Multiply every divisor before with the new primefactors
//and append them to the list
k := pfpotMax[i];
p := SmallPrimes[pfpotPrimIdx[i]];
pPot :=1;
repeat
pPot *= p;
For j := 0 to len-1 do
Begin
pDivs[l]:= pPot*pDivs[j];
inc(l);
end;
dec(k);
until k<=0;
len := l;
end;
p := pfRemain;
If p >1 then
begin
For j := 0 to len-1 do
Begin
pDivs[l]:= p*pDivs[j];
inc(l);
end;
len := l;
end;
end;
//Sort. Insertsort much faster than QuickSort in this special case
InsertSort(pDivs,0,len-1);
//end marker
pDivs[len] :=0;
end;

procedure AllFacsOut(var Divs:tdivisors;proper:boolean=true);
var
k,j: Int32;
Begin
k := 0;
j := 1;
if Proper then
j:= 2;
repeat
IF Divs[j] = 0 then
BREAK;
write(Divs[k],',');
inc(j);
inc(k);
until false;
writeln(Divs[k]);
end;

BEGIN
InitSmallPrimes;
END.
```

Pages in category "PrimTrial"

The following 11 pages are in this category, out of 11 total.