ubigintsv2.pas

Delphi的大数运算演示 pudn上大多是VC的 所以传个Delphi的

    SetLength(fDigits, 20);
    n := system.abs(i2);
    i := 0;
    repeat
      // start add by Charles Doumar
      nn := n div Base;
      fDigits[i] := n - nn * Base;
      n  := nn;
      // end
      //      fdigits[i]:=n mod base;
      //      n:=n div base;
      Inc(i);
    until n = 0;
    if i2 < 0 then
      Sign := -1
    else if i2 = 0 then
      Sign := 0
    else if i2 > 0 then
      Sign := +1;
    SetLength(fDigits, i);
    Trim;
  end;
end;

{************* Assign   (String type *********}
procedure TInteger.Assign(const i2: string);
{Convert a  string number}
var
  i, j:    integer;
  zeroval: boolean;
  n, nn:   int64;
  pos:     integer;
begin
  n := length(I2) div GetBasePower + 1;
  SetLength(fDigits, n);
  for i := 0 to n - 1 do
    fDigits[i] := 0;
  Sign := +1;
  j   := 0;
  zeroval := True;
  n   := 0;
  pos := 1;
  for i := length(i2) downto 1 do
  begin
    if i2[i] in ['0'..'9'] then
    begin
      n   := n + pos * (Ord(i2[i]) - Ord('0'));
      pos := pos * 10;
      if pos > Base then
      begin
        // start add by Charles Doumar
        nn  := n div Base;
        fDigits[j] := n - nn * Base;
        n   := nn;
        // end
        //          fdigits[j]:= n mod base;
        //          n:= n div base;
        pos := 10;
        Inc(j);
        zeroval := False;
      end
      else;
    end
    else if i2[i] = '-' then
      Sign := -1;
  end;
  fDigits[j] := n;  {final piece of the number}
  if zeroval and (n = 0) then
    Sign := 0;
  Trim;
end;

{************ Add *********}
procedure TInteger.Add(const I2: TInteger);
{add - TInteger}
var
  ii: TInteger;
begin
  if i2 = self then
  begin
    ii := TInteger.Create;
    ii.Assign(i2);
    if Sign <> ii.Sign then
      AbsSubtract(ii)
    else
      AbsAdd(ii);
    ii.Free;
  end
  else
  begin
    if Sign <> i2.Sign then
      AbsSubtract(i2)
    else
      AbsAdd(i2);
  end;
end;


{**************** AbsAdd ***************}
procedure tinteger.AbsAdd(const i2: tinteger);
{add values ignoring signs}
var
  i: integer;
  n, Carry: int64;
begin
  I3.Assign(self);
  SetLength(fDigits, max(length(fDigits), length(i2.fDigits)) + 1);
  {"add" could grow result by two digit}
  i     := 0;
  Carry := 0;
  while i < min(length(i2.fDigits), length(i3.fDigits)) do
  begin
    n     := i2.fDigits[i] + i3.fDigits[i] + Carry;
    // Start add by Charles Doumar
    Carry := n div Base;
    fDigits[i] := n - Carry * Base;
    // end
    //      fdigits[i]:= n mod base;
    //      if n >= base then carry:=1 else carry:=0;
    Inc(i);
  end;
  if length(i2.fDigits) > length(i3.fDigits) then
    while i <{=}length(i2.fDigits) do
    begin
      n     := i2.fDigits[i] + Carry;
      // Start add by Charles Doumar
      Carry := n div Base;
      fDigits[i] := n - Carry * Base;
      // end
      //      fdigits[i]:= n mod base;
      //      if n >= base then carry:=1 else carry:=0;
      Inc(i);
    end
  else if length(i3.fDigits) > length(i2.fDigits) then
  begin
    while i <{=}length(i3.fDigits) do
    begin
      n     := i3.fDigits[i] + Carry;
      // Start add by Charles Doumar
      Carry := n div Base;
      fDigits[i] := n - Carry * Base;
      // end
      //      fdigits[i]:= n mod base;
      //      if n >= base then carry:=1 else carry:=0;
      Inc(i);
    end;
  end;
  fDigits[i] := Carry;
  Trim;
end;


{************* Add (int64) ********}
procedure TInteger.Add(const I2: int64);
{Add - Int64}
begin
  IAdd3.Assign(I2);
  Add(IAdd3);
end;

{*************** AbsSubtract *************}
procedure TInteger.AbsSubtract(const i2: Tinteger);
{Subtract values ignoring signs}
var
  c:  integer;
  i3: TInteger;
  i, j, k: integer;
begin {request was subtract and signs are same,
         or request was add and signs are different}
  c  := AbsCompare(i2);
  i3 := TInteger.Create;
  if c < 0 then {abs(i2) larger, swap and subtract}
  begin
    i3.Assign(self);
    Assign(i2);
  end
  else if c >= 0 then
    {self is bigger} i3.Assign(i2);
  for i := 0 to high(i3.fDigits) do
  begin
    if fDigits[i] >= i3.fDigits[i] then
      fDigits[i] := fDigits[i] - i3.fDigits[i]
    else
    begin  {have to "borrow"}
      j := i + 1;
      while (j <= high(fDigits)) and (fDigits[j] = 0) do
        Inc(j);
      if j <= high(fDigits) then
      begin
        for k := j downto i + 1 do
        begin
          Dec(fDigits[k]);
          fDigits[k - 1] := fDigits[k - 1] + Base;
        end;
        fDigits[i] := fDigits[i] - i3.fDigits[i];
      end
      else
        ShowMessage('Subtract error');
    end;
  end;
  i3.Free;
  Trim;
end;


{*************** Mult  (Tinteger type) *********}
procedure TInteger.Mult(const I2: TInteger);
{Multiply - by Tinteger}
const
  ConstShift = 48;
var
  Carry, n, product: int64;
  xstart, ystart, i, j, k: integer;
begin
  xstart := high(self.fDigits);
  ystart := high(i2.fDigits);
  imult1.AssignZero;
  imult1.Sign := i2.Sign * Sign;
  SetLength(imult1.fDigits, xstart + ystart + 3);
  // long multiply ignoring base
  for i := 0 to xstart do
  begin
    Carry := 0;
    for j := 0 to ystart do
    begin
      k     := i + j;
      product := i2.fDigits[j] * self.fDigits[i] + imult1.fDigits[k] + Carry;
      Carry := product shr ConstShift;
      if Carry = 0 then
        imult1.fDigits[k] := product
      else
      begin
        Carry := product div Base;
        imult1.fDigits[k] := product - Carry * Base;
      end;
    end;
    imult1.fDigits[ystart + i + 1] := Carry;
  end;
  // place in proper base
  xstart := length(imult1.fDigits) - 1;
  Carry  := 0;
  for i := 0 to xstart - 1 do
  begin
    n     := imult1.fDigits[i] + Carry;
    Carry := n div Base;
    imult1.fDigits[i] := n - Carry * Base;
  end;
  imult1.fDigits[xstart] := Carry;
  Assign(imult1);
  //  Trim; trim in assign
end;

{*************** Mult  (Int64 type) *********}
procedure TInteger.Mult(const I2: int64);
{Multiply - by int64}
var
  Carry, n, d: int64;
  i:     integer;
  ITemp: TInteger;
begin
  d := system.abs(i2);
  if d > $ffffffff then
  begin
    itemp := TInteger.Create;
    itemp.Assign(i2);
    self.Mult(itemp);
    itemp.Free;
    exit;
  end;
  Carry := 0;
  for i := 0 to high(fDigits) do
  begin
    n     := fDigits[i] * d + Carry;
    Carry := n div Base;
    fDigits[i] := n - Carry * Base;
  end;
  if Carry <> 0 then
  begin
    i := high(fDigits) + 1;
    SetLength(fDigits, i + 11 div GetBasePower + 1);
    while Carry > 0 do
    begin
      n     := Carry;
      Carry := n div Base;
      fDigits[i] := n - Carry * Base;
      Inc(i);
    end;
  end;
  Trim;
end;


{************ Divide *************}
procedure TInteger.Divide(const I2: TInteger);
{Divide - by TInteger}
var
  dummy: int64;
begin
  //if i2.CompareZero = 0 then
  //  exit;
  //add by hk
  if high(i2.fDigits) = 0 then
    divmodsmall(i2.Sign * i2.fDigits[0], dummy)
  else
    {IDiv3 holds the remainder (which we don't need)}
    DivideRem(I2, idiv3);
end;

{************* Divide (Int64) **********}
procedure TInteger.Divide(const I2: int64);
{Divide - by Int64}
var
  dummy: int64;
begin
  //add by hk
  if i2 = 0 then
    exit;
  if system.abs(i2) < Base then
    divmodsmall(i2, dummy)
  else
  begin
    idiv2.Assign(i2);
    Divide(idiv2);
  end;
end;

{***************** Modulo *************}
procedure Tinteger.Modulo(const i2: TInteger);
{Modulo (remainder after division) - by Tinteger}
var
  k: int64;
begin
  //add by hk
  if high(i2.fDigits) = 0 then
  begin
    divmodsmall(i2.Sign * i2.fDigits[0], k);
    assignsmall(k);
  end
  else
  begin
    DivideRem(i2, imod3);
    Assign(imod3);
  end;
end;

{************ Modulo (Int64) ************}
procedure TInteger.Modulo(const I2: int64);
{Modulo - by Int64}
var
  i3: Tinteger;
  k:  int64;
begin
  //add by hk
  if i2 = 0 then
    AssignZero
  else if system.abs(i2) < Base then
  begin
    divmodsmall(i2, k);
    assignsmall(k);
  end
  else
  begin
    i3 := TInteger.Create;
    i3.Assign(i2);
    Modulo(i3);
    i3.Free;
  end;
end;


{**************** Dividerem ***************}
procedure TInteger.DivideRem(const I2: TInteger; var remain: TInteger);
    (*
    This version is based on a paper "Multiple-length Division Revisited: a Tour
    of the Minefield", by Per Brinch Hansen, Software - Practice and Experience,
    Vol 24(6), June 1994.

    Efficient implementation of long division
    *)
{Product}
  procedure product(var x: TInteger; y: TInteger; k: integer);
  var
    {carry,} i:     integer;
    Carry, m, temp: int64;
  begin
    // multiple-length division revisited
    m := y.GetLength;
    x.AssignZero;
    Carry := 0;
    if length(x.fDigits) <= m then
      SetLength(x.fDigits, m + 1);
    for i := 0 to m - 1 do
    begin
      temp  := y.fDigits[i] * k + Carry;
      // start add by Charles Doumar
      Carry := temp div BaseVal;
      x.fDigits[i] := temp - Carry * BaseVal;
      // end
      //      x.fdigits[i] := temp mod baseval;
      //      carry := temp div baseval
    end;
    (*if m <= x.getlength{w} then*) x.fDigits[m] := Carry;
    //else if carry <> 0 then showmessage('Product Overflow');
  end;

  procedure Quotient(var x: TInteger; y: TInteger; k: integer);
  var
    i, m: integer;
    temp, Carry: int64;
  begin
    m := y.GetLength;
    x.AssignZero;
    Carry := 0;
    SetLength(x.fDigits, m);
    for i := m - 1 downto 0 do
    begin
      temp  := Carry * BaseVal + y.fDigits[i];
      x.fDigits[i] := temp div k;
      // start add by Charles Doumar
      Carry := temp - x.Digits[i] * k;
      // end
      //carry := temp mod k
    end;
  end;

  procedure Remainder(var x: TInteger; y: TInteger; k: integer);
  var
    Carry, n, temp: int64;
    i, m: integer;
  begin
    m := y.GetLength;
    x.AssignZero;
    Carry := 0;
    SetLength(x.fDigits, M);
    for i := m - 1 downto 0 do
      //      carry := (carry * baseval + y.fdigits[i]) mod k;
      //      Start update by Charles Doumar
    begin
      n     := (Carry * BaseVal + y.fDigits[i]);
      temp  := n div k;
      Carry := n - temp * k;
    end;
    //       end
    x.fDigits[0] := Carry;
  end;

  function Trial(r, d: TInteger; k, m: integer): int64;
  var
    d2, r3: int64;
    km:     integer;
  begin
    {2 <= m <= k+m <= w}
    km := k + m;
    if length(r.fDigits) < km + 1 then
      SetLength(r.fDigits, km + 1);
    r3     := (r.fDigits[km] * BaseVal + r.fDigits[km - 1]) * BaseVal +
      r.fDigits[km - 2];
    d2     := d.fDigits[m - 1] * BaseVal + d.fDigits[m - 2];
    Result := min(r3 div d2, BaseVal - 1);
  end;

  function Smaller(r, dq: TInteger; k, m: integer): boolean;
  var
    i, j: integer;
  begin
    {0 <= k <= k+m <= w}
    i := m;
    j := 0;
    while i <> j do
      if r.fDigits[i + k] <> dq.fDigits[i] then
        j := i
      else
        Dec(i);
    Result := r.fDigits[i + k] < dq.fDigits[i];
  end;

  procedure Difference(var r: TInteger; dq: TInteger; k, m: integer);
  var
    borrow, diff, i: integer;
    acarry: int64;
  begin
    {0 <= k <= k+m <= w}
    if length(r.fDigits) < m + k + 1 then
      SetLength(r.fDigits, m + k + 1);
    if length(dq.fDigits) < m + 1 then
      SetLength(dq.fDigits, m + 1);
    borrow := 0;
    for i := 0 to m do
    begin
      diff   := r.fDigits[i + k] - dq.fDigits[i] - borrow + BaseVal;
      // start add by Charles Doumar
      acarry := diff div BaseVal;
      r.fDigits[i + k] := diff - acarry * BaseVal;
      borrow := 1 - acarry;
      // end
      //      r.fdigits[i + k] := diff mod baseval;
      //      borrow := 1 - diff div baseval
    end;
    if borrow <> 0 then
      ShowMessage('Difference Overflow');
  end;

  procedure LongDivide(x, y: TInteger; var q, r: TInteger; const n, m: integer);
  var
    f, k: integer;
    qt:   int64;
  begin
    {2 <= m <= n <= w}
    f := BaseVal div (y.fDigits[m - 1] + 1);

    product(r, x, f);
    //r.assign(x); r.mult(f);
    product(d, y, f);
    //d.assign(y); d.mult(f);
    q.AssignZero;
    SetLength(q.fDigits, n - m + 1);
    for k := n - m downto 0 do
    begin
      {2 <= m <= k+m <=n <= w}
      qt := trial(r, d, k, m);
      product(dq, d, qt);
      //dq.assign(d); dq.mult(qt);
      if length(dq.fDigits) < M + 1 then
        SetLength(dq.fDigits, M + 1);
      if smaller(r, dq, k, m) then
      begin
        qt := qt - 1;
        product(dq, d, qt);
        //dq.assign(d); dq.mult(qt);
      end;
      if k > high(q.fDigits) then
        SetLength(q.fDigits, k + 1);
      q.fDigits[k] := qt;
      difference(r, dq, k, m);
    end;
    idiv4.Assign(r);
    Quotient(r, idiv4, f);
    r.Trim;
  end;

  procedure Division(x, y: TInteger; var q, r: TInteger);
  var
    m, n, y1: integer;
  begin
    m := y.GetLength;
    if m = 1 then
    begin
      y1 := y.fDigits[m - 1];
      if y1 > 0 then
      begin
        Quotient(q, x, y1);
        Remainder(r, x, y1);
      end
      else
        ShowMessage('Division Overflow');
    end
    else
    begin
      n := x.GetLength;
      if m > n then
      begin
        q.AssignZero;
        r := x;
      end
      else {2 <= m <= n <= w}
        longdivide(x, y, q, r, n, m);
    end;
  end;

var
  signout: integer;

begin
  if Sign <> i2.Sign then
    signout := -1
  else
    signout := +1;
  if not self.IsZero then
    Sign := +1;
  if not i2.IsZero then
    i2.Sign := +1;
  //  if i2.compare(0) = 0 then
  if i2.CompareZero = 0 then
    // start add by Charles Doumar
  begin
    remain.AssignZero;
    exit;
  end;

  if AbsCompare(i2) >= 0 then  {dividend>=divisor}
  begin
    idivd3.Assign(self);
    division(idivd3, i2, self, remain);