ubigintsv2.pas

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

  gcdarr[a].Free;
  gcdarr[b].Free;
end;


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

{*********** IsOdd *********}
function Tinteger.IsOdd: boolean;
  {Return true if seld is an odd integer}
begin
  //  result:=(fdigits[0] mod 2)=1
  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
  {i2,i3,i4,}n_1, one: 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 := tinteger.Create;
    h  := 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;
  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;
  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;
      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;
begin
  r := tinteger.Create;

  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;

  q  := tinteger.Create;
  fn := tinteger.Create;
  fv := tinteger.Create;
  h  := tinteger.Create;
  n  := tinteger.Create;
  m  := tinteger.Create;

  n.Assign(self);
  m.Assign(I2);
  r.AssignOne;
  fv.AssignZero;
  fn.AssignOne;
  repeat
    q.Assign(m);
    q.DivideRem(n, r);
    if r.CompareZero <> 0 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.CompareZero = 0;
  h.Add(i2);
  h.Modulo(i2);
  Assign(h);
  r.Free;
  q.Free;
  fn.Free;
  fv.Free;
  h.Free;
  n.Free;
  m.Free;
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 + 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.Pow(root);
  if Compare(parta) < 0 then
  begin
    AssignOne;
    parta.Free;
    exit;
  end;
  PartB := tinteger.Create;
  PartC := tinteger.Create;
  Guess := 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;
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;

procedure TInteger.divide2;
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;

procedure TInteger.divmodsmall(d: int64; var rem: int64);
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.Assign(d);
    remh.Assign(rem);
    DivideRem(dh, remh);
    remh.ConvertToInt64(rem);
    dh.Free;
    remh.Free;
  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;
    self.Sign := self.Sign * dsign;
    if self.Sign < 0 then
      rem := -rem;
    Trim;
  end;
end;

procedure TInteger.ShiftLeftBase10(num: integer);
var
  b, i, len, len1, bb, Top, digs, add1, Mul: integer;
  Last, Carry, Sum, Mul1: int64;
begin
  bb  := self.Base;
  Top := high(fDigits);
  if Top < 0 then
  begin
    AssignZero;
    exit;
  end;
  last := fDigits[Top];
  if (num > 0) and (Last > 0) then
  begin
    b    := GetBasePower;
    len  := (num - 1) div b;
    len1 := (num) div b;
    shiftleftNum(len1);
    if len = len1 then
    begin
      Digs := 1 + system.trunc(Math.log10(Last));
      Mul  := num - b * len;  //1,2 for b=3
      case Mul of
        0: Mul1 := 1;
        1: Mul1 := 10;
        2: Mul1 := 100;
        3: Mul1 := 1000;
        4: Mul1 := 10000;
        5: Mul1 := 100000;
        6: Mul1 := 1000000;
        7: Mul1 := 10000000;
        8: Mul1 := 100000000;
        else
          Mul1 := system.trunc(Math.intpower(10, Mul));
      end;
      Add1 := (digs + mul - 1) div b;
      if add1 > 0 then
        SetLength(fDigits, Top + Len + 2);
      Carry := 0;
      for i := len1 to high(fDigits) do
      begin
        Sum   := fDigits[i] * Mul1 + Carry;
        Carry := Sum div bb;
        fDigits[i] := sum - bb * Carry;
      end;
    end;
  end;
end;

procedure TInteger.ShiftRightBase10(num: integer);
var
  b, len, len1, Top, digs, Mul: integer;
  Last: int64;
begin
  b   := GetBasePower;
  Top := High(fDigits);
  if Num > ((Top + 1) * b) then
    self.AssignZero
  else
  begin
    Last := fDigits[Top];
    if (num > 0) and (last > 0) then
    begin
      len  := (num - 1) div b;
      len1 := (num) div b;
      if len <> len1 then
        self.ShiftRightNum(len1)
      else
      begin //len1=len2
        digs := num div b;
        mul  := num - digs * b;
        shiftleftbase10(b - mul);
        shiftrightNum(len1 + 1);
      end;
    end;
  end;
end;

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

initialization
  SetbaseVal(100000);
  randomize;

finalization

  worki2.Free;
  imult1.Free;
  imult2.Free;
  isub3.Free;
  iadd3.Free;
  idiv2.Free;
  idiv3.Free;
  idivd3.Free;
  idiv4.Free;
  icomp3.Free;
  imod3.Free;
  d.Free;
  dq.Free;
  i3.Free;
  i4.Free;
  // end
end.