Category:PrimTrial

From Rosetta Code

Units for Freepascal

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}
  {$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.

UNIT for prime decomposition

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.