RosettaCodeData/Task/Smith-numbers/Pascal/smith-numbers.pas

program SmithNum;
{$IFDEF FPC}
  {$MODE objFPC} //result and  useful for x64
  {$CODEALIGN PROC=64}
{$ENDIF}
uses
  sysutils;
type
  tdigit  = byte;
  tSum    = LongInt;
const
  base = 10;
  //maxDigitCnt *(base-1) <= High(tSum)
  //maxDigitCnt <= High(tSum) DIV (base-1);
  maxDigitCnt = 16;

  StartPrimNo = 6;
  csegsieveSIze = 2*3*5*7*11*13;//prime 0..5
type
  tDgtSum = record
              dgtNum : LongInt;
              dgtSum : tSum;
              dgts   : array[0..maxDigitCnt-1] of tdigit;
            end;
  tNumFactype = word;
  tnumFactor = record
                 numfacCnt: tNumFactype;
                 numfacts : array[1..15] of tNumFactype;
               end;
  tpnumFactor= ^tnumFactor;

  tsieveprim = record
                 spPrim   : Word;
                 spDgtsum : Word;
                 spOffset : LongWord;
               end;
  tpsieveprim = ^tsieveprim;

  tsievePrimarr  = array[0..6542-1] of tsieveprim;
  tsegmSieve     = array[1..csegsieveSIze] of tnumFactor;

var
  Primarr:tsievePrimarr;
  copySieve,
  actSieve : tsegmSieve;
  PrimDgtSum :tDgtSum;
  PrimCnt : NativeInt;

function IncDgtSum(var ds:tDgtSum):boolean;
//add 1 to dgts and corrects sum of Digits
//return if overflow happens
var
  i : NativeInt;
Begin
  i := High(ds.dgts);
  inc(ds.dgtNum);
  repeat
    IF ds.dgts[i] < Base-1 then
    //add one and done
    Begin
      inc(ds.dgts[i]);
      inc(ds.dgtSum);
      BREAK;
    end
    else
    Begin
      ds.dgts[i] := 0;
      dec(ds.dgtSum,Base-1);
    end;
    dec(i);
 until i < Low(ds.dgts);
 result := i < Low(ds.dgts)
end;

procedure OutDgtSum(const ds:tDgtSum);
var
  i : NativeInt;
Begin
  i := Low(ds.dgts);
  repeat
    write(ds.dgts[i]:3);
    inc(i);
  until i > High(ds.dgts);
  writeln(' sum of digits :  ',ds.dgtSum:3);
end;

procedure OutSieve(var s:tsegmSieve);
var
  i,j : NativeInt;
Begin
  For i := Low(s) to High(s) do
    with s[i] do
    Begin
      write(i:6,numfacCnt:4);
      For j := 1 to numfacCnt do
        write(numFacts[j]:5);
      writeln;
    end;
end;

procedure SieveForPrimes;
// sieve for all primes < High(Word)
var
  sieve : array of byte;
  pS : pByte;
  p,i   : NativeInt;
Begin
  setlength(sieve,High(Word));
  Fillchar(sieve[Low(sieve)],length(sieve),#0);
  pS:= @sieve[0]; //zero based
  dec(pS);// make it one based
  //sieve
  p := 2;
  repeat
    i := p*p;
    IF i> High(Word) then
      BREAK;
    repeat pS[i] := 1; inc(i,p); until i > High(Word);
    repeat inc(p) until pS[p] = 0;
  until false;
  //now fill array of primes
  fillchar(PrimDgtSum,SizeOf(PrimDgtSum),#0);
  IncDgtSum(PrimDgtSum);//1
  i := 0;
  For p := 2 to High(Word) do
  Begin
    IncDgtSum(PrimDgtSum);
    if pS[p] = 0 then
    Begin
      with PrimArr[i] do
      Begin
        spOffset := 2*p;//start at 2*prime
        spPrim   := p;
        spDgtsum := PrimDgtSum.dgtSum;
      end;
      inc(i);
    end;
  end;
  PrimCnt := i-1;
end;

procedure MarkWithPrime(SpIdx:NativeInt;var sf:tsegmSieve);
var
  i : NativeInt;
  pSf :^tnumFactor;
  MarkPrime : NativeInt;
Begin
  with Primarr[SpIdx] do
  Begin
    MarkPrime := spPrim;
    i :=  spOffSet;
    IF i <= csegsieveSize then
    Begin
      pSf := @sf[i];
      repeat
        pSf^.numFacts[pSf^.numfacCnt+1] := SpIdx;
        inc(pSf^.numfacCnt);
        inc(pSf,MarkPrime);
        inc(i,MarkPrime);
      until i > csegsieveSize;
    end;
    spOffset := i-csegsieveSize;
  end;
end;

procedure InitcopySieve(var cs:tsegmSieve);
var
  pr: NativeInt;
Begin
  fillchar(cs[Low(cs)],sizeOf(cs),#0);
  For Pr := 0 to 5 do
  Begin
    with Primarr[pr] do
     spOffset := spPrim;//mark the prime too
    MarkWithPrime(pr,cs);
  end;
end;

procedure MarkNextSieve(var s:tsegmSieve);
var
  idx: NativeInt;
Begin
  s:= copySieve;
  For idx := StartPrimNo to PrimCnt do
    MarkWithPrime(idx,s);
end;

function DgtSumInt(n: NativeUInt):NativeUInt;
var
  r : NativeUInt;
Begin
  result := 0;
  repeat
    r := n div base;
    inc(result,n-base*r);
    n := r
  until r = 0;
end;

{function DgtSumOfFac(pN: tpnumFactor;dgtNo:tDgtSum):boolean;}
function TestSmithNum(pN: tpnumFactor;dgtNo:tDgtSum):boolean;
var
  i,k,r,dgtSumI,dgtSumTarget : NativeUInt;
  pSp:tpsieveprim;
  pNumFact : ^tNumFactype;
Begin
  i := dgtNo.dgtNum;
  dgtSumTarget :=dgtNo.dgtSum;

  dgtSumI := 0;
  with pN^ do
  Begin
    k := numfacCnt;
    pNumFact := @numfacts[k];
  end;

  For k := k-1 downto 0 do
  Begin
    pSp := @PrimArr[pNumFact^];
    r := i DIV pSp^.spPrim;
    repeat
      i := r;
      r := r DIV pSp^.spPrim;
      inc(dgtSumI,pSp^.spDgtsum);
    until (i - r* pSp^.spPrim) <> 0;
    IF dgtSumI > dgtSumTarget then
    Begin
      result := false;
      EXIT;
    end;
    dec(pNumFact);
  end;
  If i <> 1 then
    inc(dgtSumI,DgtSumInt(i));
  result := dgtSumI = dgtSumTarget
end;

function CheckSmithNo(var s:tsegmSieve;var dgtNo:tDgtSum;Lmt:NativeInt=csegsieveSIze):NativeUInt;
var
  pNumFac : tpNumFactor;
  i : NativeInt;
Begin
  result := 0;
  i := low(s);
  pNumFac := @s[i];
  For i := i to lmt do
  Begin
    incDgtSum(dgtNo);
    IF pNumFac^.numfacCnt<> 0 then
      IF TestSmithNum(pNumFac,dgtNo) then
      Begin
        inc(result);
        //Mark as smith number
        inc(pNumFac^.numfacCnt,1 shl 15);
      end;
    inc(pNumFac);
  end;
end;

const
  limit = 100*1000*1000;
var
  actualNo :tDgtSum;
  i,s : NativeInt;
Begin
  SieveForPrimes;
  InitcopySieve(copySieve);
  i := 1;
  s:= -6;//- 2,3,5,7,11,13

  fillchar(actualNo,SizeOf(actualNo),#0);
  while i < Limit-csegsieveSize do
  Begin
    MarkNextSieve(actSieve);
    inc(s,CheckSmithNo(actSieve,actualNo));
    inc(i, csegsieveSize);
  end;
  //check the rest
  MarkNextSieve(actSieve);
  inc(s,CheckSmithNo(actSieve,actualNo,Limit-i+1));
  write(s:8,' smith-numbers up to ',actualNo.dgtnum:10);
end.