ubigintsv3.pas

Delphi for fun library v12, latest. This is the library for manuplating list, combination-permutati

{************* IsZero *************}
function TInteger.IsZero: Boolean;
begin
  result := Sign=0;
end;

{************* AbsCompare *************}
function TInteger.AbsCompare(i2: Tinteger): integer;
  {compare absolute values ingoring signs - to Tinteger}
var
  i: integer;
begin

  Result := 0;
  if (self.Sign = 0) and (i2.Sign = 0) then
    Result := 0
  else if length(fDigits) > length(i2.fDigits) then
    Result := +1
  else if length(fDigits) < length(i2.fDigits) then
    Result := -1
  else {equal length}
    for i := high(fDigits) downto 0 do
    begin
      if fDigits[i] > i2.fDigits[i] then
      begin
        Result := +1;
        break;
      end
      else if fDigits[i] < i2.fDigits[i] then
      begin
        Result := -1;
        break;
      end;
    end;
end;

{************* AbsCompare *************}
function TInteger.AbsCompare(I2: int64): integer;
  {compare absolute values ingoring signs - to Tinteger}
var
  i3: Tinteger;
begin
  //i3 := tinteger.Create;
  //i3.Assign(i2);
  i3:=GetNextScratchPad(i2);
  Result := AbsCompare(i3);
  releaseScratchPad(i3);
  //i3.Free;
end;


{*********** Factorial *******}
procedure TInteger.Factorial;
{Compute factorial - number must be less than max integer value}
var
  n: int64;
  i: integer;
begin
  n := 0;
  if (Compare(high(integer)) >= 0) or (Sign < 0) then
    exit;
  if IsZero then
  begin  {0! =1 by definition}
    AssignOne;
    exit;
  end;
  for i := high(fDigits) downto 0 do
  begin
    n := n * Base + fDigits[i];
  end;
  Dec(n);
  while n > 1 do
  begin
    Mult(n);
    Dec(n);
    {provide a chance to cancel long running ops}
    if (n and $4f) = $4f then
      application.ProcessMessages;
  end;
end;

{************** ConvertToDecimalStirng ********}
function TInteger.ConvertToDecimalString(commas: boolean): string;
var
  i, j, NumCommas, CurPos, StopPos, b, DigCount, NewPos, Top: integer;
  n, nn, last: int64;
  c: byte;

begin
  if length(fDigits) = 0 then
    AssignZero;
  if IsZero then
  begin
    Result := '0';
    exit;
  end;
  Result := '';
  b      := GetBasePower;
  Top    := high(self.Digits);
  Last   := fDigits[Top];
  DigCount := Top * b + 1 + trunc(Math.log10(Last));
  Dec(Top);
  if Sign > 0 then
  begin
    CurPos := DigCount;
    SetLength(Result, CurPos);
    StopPos := 0;
  end
  else
  begin
    CurPos := DigCount + 1;
    SetLength(Result, CurPos);
    Result[1] := '-';
    StopPos   := 1;
  end;
  for i := 0 to Top do
  begin
    n := fDigits[i];
    for j := 1 to b do
    begin
      nn := n div 10;
      c  := n - nn * 10;
      n  := nn;
      Result[CurPos] := char($30 + c);
      Dec(CurPos);
    end;
  end;
  repeat
    nn   := Last div 10;
    c    := Last - nn * 10;
    Last := nn;
    Result[CurPos] := char($30 + c);
    Dec(CurPos);
  until CurPos <= StopPos;

  if Commas = True then
  begin
    CurPos    := Length(Result);
    NumCommas := (DigCount - 1) div 3;
    if NumCommas > 0 then
    begin
      NewPos := CurPos + NumCommas;
      SetLength(Result, NewPos);
      for i := 1 to NumCommas do
      begin
        Result[NewPos]     := Result[CurPos];
        Result[NewPos - 1] := Result[CurPos - 1];
        Result[NewPos - 2] := Result[CurPos - 2];
        Result[NewPos - 3] := ThousandSeparator;
        Dec(NewPos, 4);
        Dec(CurPos, 3);
      end;
    end;
  end;
end;


{********* ConvertToInt64 **********}
function TInteger.ConvertToInt64(var n: int64): boolean;
var
  i: integer;
  savesign: integer;
begin
  Result   := False;
  savesign := Sign;
  Sign     := +1;
  if (self.Compare(high(n) {9223372036854775806}) < 1) then
  begin
    n := 0;
    for i := high(fDigits) downto 0 do
      n := Base * n + fDigits[i];
    if savesign < 0 then n := -n;
    Result := True;
  end;
  Sign := savesign;
end;




{********* Pow **************}
procedure Tinteger.Pow(const exponent: int64);
{raise self to the "exponent" power}
var
  i2: tinteger;
  n:  int64;
  s:  integer;
begin
  n := exponent;
  if (n <= 0) then
  begin
    if n = 0 then
      AssignOne
    else
      AssignZero;
    exit;
  end;
  s := Sign;
  if ((s < 0) and not (odd(n))) then
    s := 1;
  Sign := 1;

  //i2 := TInteger.Create;
  //i2.AssignOne;
  i2:=GetNextScratchPad(1);
  if (n >= 1) then
    if n = 1 then
      i2.Assign(self)
    else
    begin
      repeat
        if (n and $1) = 1 then i2.Mult(self);
        n := N shr 1;
        application .processmessages;
        Square;
      until (n = 1);
      i2.Mult(self);
    end;
  Assign(i2);
  Sign := s;
  //i2.Free;
  releaseScratchPad(i2);
end;



//partially rewritten by hk Oct 2005
procedure Tinteger.ModPow(const i2, m: Tinteger);
{self^I2 modulo m}
var
  ans, e, one: Tinteger;
  hulp: integer;
begin
  {if (length(i2.fdigits) = 0) or (length(m.fdigits) = 0) then
    exit;}
  hulp := i2.Getsign;
  if hulp <= 0 then
    if hulp = 0 then
    begin
      AssignOne;
      exit;
    end
    else
      exit;
  if m.IsZero then
    exit;
  //one := tinteger.Create;
  //one.AssignOne;
  One:=GetNextScratchPad(1);
  //ans := tinteger.Create;
  //ans.AssignOne;
  Ans:=GetNextScratchpad(1);
  //e   := tinteger.Create;
  //e.Assign(i2);
  E:=GetNextScratchPad(i2);
  hulp := e.Compare(one);
  if hulp >= 0 then
    if hulp = 0 then
    begin
      Modulo(m);
      ans.Assign(self);
    end
    else
    begin
      repeat
        if e.IsOdd then
        begin
          ans.Mult(self);
          ans.Modulo(m);
        end;
        e.divide2;
        Square;
        Modulo(m);
      until (e.Compare(one) = 0);
      ans.Mult(self);
      ans.Modulo(m);
    end;
  Assign(ans);
  Releasescratchpad(e);
  Releasescratchpad(ans);
  Releasescratchpad(One);

  //ans.Free;
  //e.Free;
  //one.Free;
end;


{square root}
procedure Tinteger.Sqroot;
begin
  NRoot(2);
end;

{****************** GCD ***************}
procedure Tinteger.Gcd(const I2: tinteger);
 {greatest common divisor}
 {revised by Hans Aug 2005 to handle 0 "I2" value }
 {revised by hk Oct 2005: swapping by hashing}
var
  gcdarr:  array[0..1] of tinteger;
  a, b, h: integer;
begin
  gcdarr[0] := GetNextScratchPad(self); //tinteger.Create;
  gcdarr[1] := GetNextScratchPad(i2); //tinteger.Create;
  if AbsCompare(i2) = 1 then
  begin
    //gcdarr[0].Assign(self);  {already assigned above}
    //gcdarr[1].Assign(i2);
  end
  else
  begin
    gcdarr[1].Assign(self);
    gcdarr[0].Assign(i2);
  end;
  if gcdarr[1].IsZero then
  begin
    AssignZero;
    //gcdarr[0].Free;
    //gcdarr[1].Free;
    releasescratchpad(gcdarr[1]);
    releasescratchpad(gcdarr[0]);
    exit;
  end;
  a := 0;
  b := 1;
  repeat
    gcdarr[a].Modulo(gcdarr[b]);
    if not gcdarr[a].IsZero then
    begin
      h := a;
      a := b;
      b := h;
    end;
  until gcdarr[a].IsZero;
  Assign(gcdarr[b]);
  self.AbsoluteValue;
  //gcdarr[a].Free;
  //gcdarr[b].Free;
  releasescratchpad(gcdarr[1]);
  releasescratchpad(gcdarr[0]);
end;


{*************** GCD (Int64) ***********}
procedure Tinteger.Gcd(const I2: int64);
var
  h: tinteger;
begin
  //h := tinteger.Create;
  //h.Assign(i2);
  h:=GetNextScratchPad(I2);
  Gcd(h);
  //h.Free;
  releasescratchpad(h);
end;

{*********** IsOdd *********}
function Tinteger.IsOdd: boolean;
  {Return true if self is an odd integer}
begin
  Result := (fDigits[0] and $1) = 1;
end;


 (*rewritten by hk Oct 2005*)
 {************** IsProbablyPrime ***********}
function Tinteger.IsProbablyPrime: boolean;
  //miller rabin probabilistic primetest with 10 random bases;
var
  n_1, one, i3, worki2: tinteger;
  Base, lastdigit, j, t, under10: integer;
  probableprime: boolean;
  p, work, s, diff, factorlimit: int64;

  function witness(Base: integer; e, n: tinteger): boolean;
  var
    it, h: tinteger;
    i:     integer;
  begin
    it := Getnextscratchpad(Base); //tinteger.Create;
    h  := GetNextScratchPad; //tinteger.Create;

    it.Assign(Base);
    it.ModPow(e, n);
    Result := True;
    for i := 1 to t do
    begin
      h.Assign(it);
      h.Square;
      h.Modulo(n);
      if h.Compare(one) = 0 then
        if it.Compare(one) <> 0 then
          if it.Compare(n_1) <> 0 then
            Result := False;
      it.Assign(h);
    end;
    if it.Compare(one) <> 0 then
      Result := False;
    //it.Free;
    //h.Free;
    releaseScratchpad(h);
    releasescratchpad(it);
  end;

  function is_smallprime: boolean;
  var
    j, w, diff: integer;
    prime: boolean;
    i: int64;
  begin
    prime := True;
    ConvertToInt64(i);
    w    := trunc(sqrt(0.0 + i));
    j    := 11;
    diff := 2;
    while prime and (j <= w) do
    begin
      if (i mod j) = 0 then
        prime := False;
      Inc(j, diff);
      if diff = 2 then
        diff := 4
      else
        diff := 2;
    end;
    Result := prime;
  end;

begin
  //one := tinteger.Create;
  //n_1 := tinteger.Create;
  one:=GetnextScratchPad(1);
  n_1:= GetNextScratchPad;
  I3:=GetNextScratchPad;
  Worki2:=GetNextScratchPad;

  probableprime := True;
  factorlimit := self.Base;
  if factorlimit > 200 then
    factorlimit := 200;
  if factorlimit = 10 then
    factorlimit := 32;
  {lastdigit := worki2.shiftright mod 10; }
  lastdigit := fDigits[0] mod 10;
  under10 := Compare(10);
  if under10 = -1 then
    if (lastdigit <> 2) and (lastdigit <> 3) and (lastdigit <> 5) and
      (lastdigit <> 7) then
    begin
      probableprime := False;
    end;
  if under10 >= 0 then
  begin
    if (lastdigit <> 1) and (lastdigit <> 3) and (lastdigit <> 7) and
      (lastdigit <> 9) then
      probableprime := False
    else
    begin
      p := 3;

      while probableprime and (p < 8) do
      begin
        worki2.Assign(self);
        worki2.divmodsmall(p, work);
        if work = 0 then
          probableprime := False;
        p := p + 4;
      end;
    end;
  end;
  if probableprime and (under10 > 0) and (Compare(120) = 1) then
    if Compare(1000000) = -1 then
      probableprime := is_smallprime
    else
    begin
      p    := 11;
      diff := 2;
      while probableprime and (p < factorlimit) do
      begin
        worki2.Assign(self);
        worki2.divmodsmall(p, work);
        if work = 0 then
          probableprime := False;
        p := p + diff;
        if diff = 2 then
          diff := 4
        else
          diff := 2;
      end;
      if probableprime then
      begin
        one.AssignOne;
        n_1.Assign(self);
        n_1.Subtract(1);
        t := 0;
        i3.Assign(n_1);
        repeat
          if not (i3.IsOdd) then
          begin
            Inc(t);
            i3.divide2;
          end;
        until i3.IsOdd;
        worki2.Assign(self);
        j := 10;
        if Compare(1000000) > 0 then
          s := 1000000
        else
          worki2.ConvertToInt64(s);
        s := s - 1;
        while (j > 0) and probableprime do
        begin
          repeat
            Base := random(s);
          until Base > 1;
          probableprime := witness(Base, i3, worki2);
          Dec(j);
        end;
        //n_1.Free;
        //one.Free;
        releaseScratchPad(worki2);
        releaseScratchPad(i3);
        releasescratchpad(n_1);
        releasescratchpad(one);

      end;
    end;
  Result := probableprime;
end;



{************ InvMod **********}
procedure Tinteger.InvMod(I2: Tinteger);
{calculates the number n such that n*self=1 mod I2, result in self}
var
  r, q, fn, fv, h, n, m: Tinteger; s1,s2,s3:integer;
begin
//Extended Euclidean Algoritm.
(* (it's not necessary to calculate the gcd beforehand,
   the extended euclidean algoritm does that itself.)
  r.Assign(self);
  r.Gcd(I2);
  {  1st operand must be non-zero and operands must be relatively prime}
  //  if (compare(0) = 0) or (r.compare(1) <> 0) then
  if (CompareZero = 0) or (r.Compare(1) <> 0) then

  begin
    AssignZero;
    // start add by Charles Doumar
    r.Free;
    // end
    exit;
  end;

 *)
 {both operands must be non-zero}
 if (IsZero) or (i2.IsZero) then
 begin
   assignzero;
   exit
 end;
  r  := Getnextscratchpad(1); //tinteger.Create;
  q  := Getnextscratchpad(0); //tinteger.Create;
  fn := Getnextscratchpad(1); //tinteger.Create;
  fv := Getnextscratchpad(0); //tinteger.Create;
  h  := Getnextscratchpad(0); //tinteger.Create;
  n  := Getnextscratchpad(self); //tinteger.Create;
  m  := Getnextscratchpad(I2); //tinteger.Create;
  //n.Assign(self);
  //m.Assign(I2);
  m.sign:=1;
  s1:=self.sign;
  //first make all operands positive;
  if s1<0 then    {change restclass}
  begin
    n.changesign;
    n.Add(m);
  end;
  n.modulo(m);
  //start euclidean algorithm
  repeat
    q.Assign(m);
    q.DivideRem(n,r);
    if not r.IsZero then
    begin
      m.Assign(n);
      n.Assign(r);
    end;
    h.Assign(fn);
    fn.Mult(q);
    fv.Subtract(fn);
    fn.Assign(fv);
    fv.Assign(h);
  until r.IsZero;
  {n contains the gcd, gcd must be 1; h contains invmod}
  if n.Compare(1)<>0 then assign(0)
  else
  begin