ubigintsv3.pas

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

  if s1<0 then h.changesign; {change restclass again if we did before}
  {Now assure that invmod has the same sign as self so
   that invmod*self will be positive}
  s3:=h.sign;
  s2:=i2.sign;
  i2.sign:=1;
  if (s3<0) and (s1>0) then h.add(i2);
  if (s3>0) and (s1<0) then h.Subtract(i2);
  Assign(h);
  i2.sign:=s2;
  end;
  (*
  r.Free;
  q.Free;
  fn.Free;
  fv.Free;
  h.Free;
  n.Free;
  m.Free;
  *)
  ReleaseScratchPad(m);
  ReleaseScratchPad(n);
  ReleaseScratchPad(h);
  ReleaseScratchPad(fv);
  ReleaseScratchPad(fn);
  ReleaseScratchPad(q);
  ReleaseScratchPad(r);

end;

{************ NRoot **********}
procedure TInteger.NRoot(const root: integer);
var
  Parta, Partb, Partc, Guess: tinteger;
  Limit, Counter, nthroot, a: integer;
  num1, num2: int64;
  str1: string;
begin
  if (root < 2) or (self.Sign <> 1) then
  begin
    AssignZero;
    exit;
  end;
  // if number within range of int64 then compute there
  if self.ConvertToInt64(num1) = True then
  begin
    num2 := trunc(power(num1, 1 / root));
    if power(num2,root)>num1 then num2:=num2-1;
    if power(num2 + 1, root) <= num1 then
      self.Assign(num2 + 1)
    else
      self.Assign(num2);
    exit;
  end;

  Limit   := length(self.fDigits) * GetBasePower * 3 + 700;
  Counter := 0;
  nthroot := root - 1;
  //PartA   := tinteger.Create;
  //parta.Assign(2);
  PartA   := GetNextScratchpad(2);
  parta.Pow(root);
  if Compare(parta) < 0 then
  begin
    AssignOne;
    //parta.Free;
    releaseScratchPad(PartA);
    exit;
  end;
  PartB := GetNextScratchPad(0); // tinteger.Create;
  PartC := GetNextScratchPad(0); //tinteger.Create;
  Guess := GetNextScratchPad(0); //Tinteger.Create;

  //initial guess
  a     := high(self.fDigits) * GetBasePower div root - 1;
  if a = 0 then
    str1 := '2'
  else
    str1 := '1';
  str1 := str1 + stringofchar('0', a);
  guess.Assign(str1);
  parta.Assign(guess);
  parta.Pow(root);
  while self.Compare(parta) > 0 do
  begin
    guess.Mult(10);
    parta.Assign(guess);
    parta.Pow(root);
  end;
  parta.Assign(guess);
  PARTA.Divide(10);
  parta.Pow(root);
  while self.Compare(parta) <= 0 do
  begin
    guess.Divide(10);
    parta.Assign(guess);
    PARTA.Divide(10);
    parta.Pow(root);
  end;
  // start of newton
  repeat
    Inc(Counter);
    parta.Assign(guess);
    parta.Mult(nthroot);
    partb.Assign(self);
    Partc.Assign(guess);
    partc.Pow(nthroot);
    partb.Divide(partc);
    parta.Add(partb);
    parta.Divide(root);
    partc.Assign(guess);  //old guess
    partb.Assign(guess);
    partb.Subtract(1);
    guess.Assign(parta); //new guess
  until (Counter > 4) and ((Counter > Limit) or (partc.Compare(guess) = 0) or
      (partb.Compare(guess) = 0));
  if Counter > Limit then
    self.AssignZero
  else
  begin
    if partc.Compare(guess) = 0 then
      self.Assign(guess)
    else
    begin
      // assign smallest value to part c;
      if partc.Compare(guess) > 0 then
        partb.Assign(guess)
      else
        partb.Assign(partc);
      // ensure not too big
      partc.Assign(partb);
      partc.Pow(root);
      while (self.Compare(partc) < 0) and (Counter < (limit + 50)) do
      begin
        partb.Subtract(1);
        partc.Assign(partb);
        partc.Pow(root);
        Inc(Counter);
      end;
      // ensure not too small
      partc.Assign(partb);
      partc.Add(1);
      partc.Pow(root);
      while (self.Compare(partc) >= 0) and (Counter < (limit + 50)) do
      begin
        partb.Add(1);
        partc.Assign(partb);
        partc.Add(1);
        partc.Pow(root);
        Inc(Counter);
      end;
      if Counter < (limit + 50) then
        self.Assign(partb)
      else
        self.AssignZero;
    end;
  end;
  //PartA.Free;
  //PartB.Free;
  //PartC.Free;
  //Guess.Free;
  ReleaseScratchpad(Guess);
  ReleaseScratchpad(PartC);
  ReleaseScratchpad(PartB);
  ReleaseScratchpad(PartA);
end;

function TInteger.GetBase: integer;
begin
  Result := Base;
end;

(*
function TInteger.IsZero: boolean;
begin
  Result := self.Digits[high(self.Digits)] = 0;
end;
*)

procedure TInteger.AssignOne;
begin
  SetLength(fDigits, 1);
  self.Sign := 1;
  self.fDigits[0] := 1;
end;

procedure TInteger.AssignZero;
begin
  SetLength(fDigits, 1);
  self.Sign := 0;
  self.fDigits[0] := 0;
end;

(**** additions by hk Oct 2005***)
procedure TInteger.assignsmall(i2: int64);
begin
  if system.abs(i2) >= Base then
    Assign(i2)
  else
  if i2 = 0 then  AssignZero
  else
  begin
    if i2 < 0 then
    begin
      self.Sign := -1;
      i2 := -i2;
    end
    else self.Sign := +1;
    SetLength(self.fDigits, 1);
    self.fDigits[0] := i2;
  end;
end;

{************** Divide2 *************}
procedure TInteger.divide2;
{Fast divide by 2 }
var
  i: integer;
begin
  for i := high(fDigits) downto 1 do
  begin
    if (fDigits[i] and 1) = 1
    then  Inc(fDigits[i - 1], Base);
    fDigits[i] := fDigits[i] shr 1;
  end;
  fDigits[0] := fDigits[0] shr 1;
  Trim;
end;

{*********** DivModSmall **************}
procedure TInteger.divmodsmall(d: int64; var rem: int64);
{DivMod for d < base }
var
  i, dsign: integer;
  k, n, nn: int64;
  dh, remh: tinteger;
begin
  if (d = 0) or (d = 1) then
  begin
    rem := 0;
    exit;
  end;
  if d = 2 then
  begin
    rem := (fDigits[0] and 1);
    divide2;
    exit;
  end;
  if system.abs(d) >= self.Base then
  begin
    //dh   := tinteger.Create;
    //remh := tinteger.Create;
    dh:=getNextScratchpad(d);
    remh:=getNextScratchPad(rem);
    //dh.Assign(d);
    //remh.Assign(rem);
    DivideRem(dh, remh);
    remh.ConvertToInt64(rem);
    //dh.Free;
    //remh.Free;
    releasescratchpad(remh);
    releaseScratchPad(dh);
  end
  else
  begin
    if d < 0 then
    begin
      dsign := -1;
      d     := -d;
    end
    else dsign := 1;
    for i := high(fDigits) downto 1 do
    begin
      n  := fDigits[i];
      nn := n div d;
      k  := n - nn * d;
      fDigits[i] := nn;
      if k > 0 then
        fDigits[i - 1] := fDigits[i - 1] + k * Base;
    end;
    n   := fDigits[0];
    nn  := n div d;
    rem := n - nn * d;
    fDigits[0] := nn;
    if self.Sign < 0 then  rem := -rem;  {Set remainder sign to dividend sign GDD}
    self.Sign := self.Sign * dsign;
    //if self.Sign < 0 then
   //   rem := -rem;
    Trim;
  end;
end;

procedure TInteger.AbsoluteValue;
begin
  if self.Sign < 0 then
    self.Sign := 1;
end;

{***************** FastMult ***************}
procedure TInteger.FastMult(const I2: TInteger);
(*
    This version is based on "fast Karatsuba multiplication", 21 Jan 1999, Carl Burch, [e-mail]
    Converted to Delphi by Charles Doumar
*)

const
  MaxArraySize = $7fffffff;
  MaxInt64ArrayElements = MaxArraySize div sizeof(int64); type
  TStaticInt64Array = array[0..MaxInt64ArrayElements - 1] of int64;
  PInt64Array = ^TStaticInt64Array;
  Int64Array = array of int64;
var
  HR: THandle;
  PA, PB, PR: PInt64Array;
  Alen, Blen, imin, imax, d: integer;

  procedure DoCarry(PR: PInt64Array; const base, d: integer);
  var
    i: integer;
    Carry: int64;
  begin
    Carry := 0;
    for i := 0 to d - 1 do
    begin
      Pr[i] := PR[i] + Carry;
      if Pr[i] < 0 then
        Carry := (Pr[i] + 1) div base - 1
      else
        Carry := PR[i] div base;
      Pr[i] := Pr[i] - Carry * base;
    end;
    if Carry <> 0 then
    begin
      pr[d - 1] := pr[d - 1] + Carry * base;
      ShowMessage('Error in FastMath DoCarry');
    end;
  end;


  procedure GradeSchoolMult(PA, PB, PR: PInt64Array; d: integer); const
  ConstShift = 48;
  var
    i, j, k,ALow,BLow: integer;
    Product, carry: int64;
  begin
    fillchar(pr[0],sizeof(int64)*2*d,0);
    Alow := d-1;
    BLow := ALow;
    while (Alow > 0) and (PA[ALow] = 0) do
      dec(ALow);
    while (BLow > 0) and (PB[Blow] = 0) do
      dec(BLow);
    k:=0;
    for i := 0 to ALow do
    begin
      Carry := 0;
      for j := 0 to BLow do
      begin
        k := i + j;
        Product := Pr[k] + PA[i] * PB[j] + Carry;
        Carry := Product shr ConstShift;
        if Carry = 0 then
          Pr[k] := product
        else begin
          Carry := product div base;
          Pr[k] := product - carry * Base;
        end;
      end;
     Pr[k + 1] :=  Pr[k + 1] + Carry;
    end;
  end;

  procedure Karatsuba(PA, PB, PR: PInt64Array; d: integer);
  var
    AL, AR, BL, BR, ASum, BSum, X1, X2, X3: PInt64Array;
    i, halfd: integer;
  begin
    if d <= 100 then
    begin
      GradeSchoolMult(Pa, Pb, Pr, d);
      exit;
    end;
    halfd := d shr 1;
    AR := @PA[0];
    AL := @PA[halfd];
    BR := @PB[0];
    BL := @PB[halfd];
    ASum := @PR[d * 5];
    BSum := @PR[d * 5 + halfd];
    X1 := @PR[0];
    X2 := @PR[d];
    X3 := @PR[d * 2];
    for i := 0 to halfd - 1 do
    begin
      asum[i] := al[i] + ar[i];
      bsum[i] := bl[i] + br[i];
    end;
    Karatsuba(AR, Br, X1, halfd);
    Karatsuba(AL, BL, X2, halfd);
    Karatsuba(ASum, BSum, X3, halfd);
    for i := 0 to d - 1 do
      X3[i] := X3[i] - X1[i] - X2[i];
    for i := 0 to d - 1 do
      PR[i + halfd] := PR[i + halfd] + X3[i];
  end;

begin
  Alen := length(self.fDigits);
  BLen := Length(I2.fDigits);
  iMin := ALen;
  iMax := BLen;
  if ALen > BLen then
  begin
    iMin := BLen;
    iMax := ALen;
  end;
  if (iMin < 80) or (imin < Imax shr 4) then
  begin
    self.Mult(i2);
    exit;
  end;
  d := 1;
  while d < iMax do
    d := d * 2;
  setlength(self.fDigits, d);
  setlength(i2.fDigits,d);
  PA := @self.fdigits[0];
  PB := @i2.fdigits[0];
  HR := GlobalAlloc(GMEM_FIXED, d * 6 * Sizeof(int64));
  PR := GlobalLock(HR);
  FillChar(Pr[0], sizeof(int64) * 6 * d, 0);
  Karatsuba(PA, PB, PR, d);
  DoCarry(Pr, self.Base, d * 2);
  d := 2 * d;
  while (d > 0) and (pr[d - 1] = 0) do
    Dec(d);
  setlength(self.fDigits, d);
  setlength(i2.fDigits, BLen);   {Changes input variable ?????????? GDD}
  move(pr[0], self.fdigits[0], d * sizeof(int64));
  self.Trim;
  self.Sign := self.Sign * i2.Sign;
  GlobalUnlock(HR);
  GlobalFree(HR);
end;

{*************** FastSquare **************}
procedure TInteger.FastSquare;
(*
    This version is based on "fast Karatsuba multiplication", 21 Jan 1999, Carl Burch, [e-mail]
    Converted to Delphi by Charles Doumar
*)

const
  MaxArraySize = $7fffffff;
  MaxInt64ArrayElements = MaxArraySize div sizeof(int64); type
  TStaticInt64Array = array[0..MaxInt64ArrayElements - 1] of int64;
  PInt64Array = ^TStaticInt64Array;
  Int64Array = array of int64;
var
  HR: THandle;
  PA, PR: PInt64Array;
  Alen, i, d, dmemory: integer;

  procedure DoCarry(PR: PInt64Array; const base, d: integer);
  var
    i: integer;
    c: int64;
  begin
    c := 0;
    for i := 0 to d - 1 do
    begin
      Pr[i] := PR[i] + c;
      if Pr[i] < 0 then
        c := (Pr[i] + 1) div base - 1
      else
        c := PR[i] div base;
      Pr[i] := Pr[i] - c * base;
    end;
    if c <> 0 then
    begin
      pr[d - 1] := pr[d - 1] + c * base;
      ShowMessage('Error in FastSquare DoCarry');
    end;
  end;

  procedure HighSchoolSquare(PA, PR: PInt64Array; d: integer); const
  ConstShift = 48;
  var
    i, j, k,ALow: integer;
    Product, carry: int64;
  begin
    fillchar(pr[0],sizeof(int64)*2*d,0);
    Alow := d-1;
    while (Alow > 0) and (PA[ALow] = 0) do
      dec(ALow);
    // Step 1. Calculate diagonal
    for i := 0 to ALow do
    begin
      k := i * 2;
      Product := PA[i] * PA[i];
      Carry := Product shr ConstShift;
      if Carry = 0 then
        PR[k] := product
      else begin
        Carry := product div base;
        PR[k] := Product - Carry * base;
        PR[k+1] := Carry;
      end;
    end;
    // Step 2. Calculate repeating part
    k:=0;
    For i := 0 to ALow do
    begin
      Carry := 0;
      for j := i+1 to ALow do
      begin
        k := i + j;
        Product := Pr[k] + PA[i] * PA[j] * 2 + Carry;
        Carry := Product shr ConstShift;
        if Carry = 0 then
          Pr[k] := product
        else begin
          Carry := product div base;
          Pr[k] := product - carry * Base;
        end;
      end;
     Pr[k + 1] :=  Pr[k + 1] + Carry;
    end;
  end;

  procedure KaratsubaSquare(PA, PR: PInt64Array; d: integer);
  var
    AL, AR, ASum, X1, X2, X3: PInt64Array;
    i, halfd: integer;
  begin
    if d <= 100 then
    begin
      HighSchoolSquare(Pa, Pr, d);
      exit;
    end;
    halfd := d shr 1;
    AR := @PA[0];
    AL := @PA[halfd];
    ASum := @PR[d * 5];
    X1 := @PR[0];
    X2 := @PR[d];
    X3 := @PR[d * 2];
    for i := 0 to halfd - 1 do
      asum[i] := al[i] + ar[i];
    KaratsubaSquare(AR, X1, halfd);
    KaratsubaSquare(AL, X2, halfd);
    KaratsubaSquare(ASum, X3, halfd);
    for i := 0 to d - 1 do
      X3[i] := X3[i] - X1[i] - X2[i];
    for i := 0 to d - 1 do
      PR[i + halfd] := PR[i + halfd] + X3[i];
  end;

begin
  Alen := length(self.fDigits);
  i := Alen;
  d := 1;
  while d < i do
    d := d * 2;
  setlength(self.fDigits, d);
  PA := @self.fdigits[0];
  dmemory := (d shr 1) * 11 * Sizeof(int64);
  HR := GlobalAlloc(GMEM_FIXED, dmemory);
  PR := GlobalLock(HR);
  FillChar(Pr[0], dmemory, 0);
  KaratsubaSquare(PA, PR, d);
  DoCarry(Pr, self.Base, d * 2);
  d := 2 * d;
  while (d > 0) and (pr[d - 1] = 0) do
    Dec(d);
  setlength(self.fDigits, d);
  move(pr[0], self.fdigits[0], d * sizeof(int64));
  self.Trim;
  self.Sign := self.Sign * self.Sign;
  GlobalUnlock(HR);
  GlobalFree(HR);
end;

var i:integer;

initialization
  SetbaseVal(100000);
  randomize;
  Lastscratchpad:=0;
  setlength(scratchpads,20);
  for i:=0 to 19 do
  begin
    scratchpads[i]:=TInteger.create;
  end;

finalization

end.