Rosetta Code

Category:PrimTrial: Difference between revisions

Category Discussion

Category:PrimTrial (view source)

Revision as of 08:44, 11 July 2023

18,520 bytes added ,  10 months ago
no edit summary
m (changed check for small primes in SmallFactor via Bit Mask 2..61)
No edit summary
 
(One intermediate revision by the same user not shown)
Line 1:
 
unit primTrial
<b>Units for Freepascal</b>
=={{header|UNIT for primTrial}}==
{{works with|Free Pascal}} {{works with|Delphi}}
Maybe NativeUint must be typed in older versions to LongWord aka cardinal
<syntaxhighlight lang="pascal">unit PrimTrial;
// NativeUInt: LongWord 32-Bit-OS/ Uint64 64-Bit-OS
{$IFDEF FPC}
{$MODE DELPHI} {$Smartlink ON} {$OPTIMIZATION ON,ALL} {$CODEALIGN proc=16}
{$Smartlink ON}
{$OPTIMIZATION ON,ALL}
{$CODEALIGN proc=32}
{$ENDIF}
 
Line 115 ⟶ 114:
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 51+1 shl 5753+1 shl 59+1 shl 61;
var
k : NativeUint;
Line 253 ⟶ 252:
InitPrime;
end.
</syntaxhighlight>
=={{header|UNIT for primsieve}}==
segmented sieve of erathostenes.Not that bad.
<syntaxhighlight lang="pascal">
unit primsieve;
{$IFDEF FPC}
{$MODE objFPC}{$Optimization ON,ALL}{$Smartlink ON}{$CODEALINGN proc=16}
{$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);
maxPreSievePrimeNum = 7;
maxPreSievePrime = 17;//smlPrimes[maxPreSievePrimeNum];
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
i := maxPreSievePrimeNum;
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
i := 1; // start with 3
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.</syntaxhighlight>
=={{header|UNIT for prime decomposition}}==
 
<syntaxhighlight lang="pascal">
unit
primeDecomp;
 
// gets factors of consecutive integers fast
// limited to 1.2e11
{$IFDEF FPC}
{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$Smartlink ON}{$COPERATORS ON}
{$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.
</syntaxhighlight>
Horst
132

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/338624...346203"