Hallo Leute ,

Ich hab ein Problem mit der Umsetzung eines
Algorithmus zur Bestimmung der Primalität einer Ganzahl.
kann mir Jemand helfen ? :

markieren

Code:

4 The Algorithm

Input: integer n > 1

 1. if (n is of the form a^b, b>1) output COMPOSITE; // n ist das Ergebnis einer Exponierung
 2. r = 2;
 3. while(r<n) {
 4.  if (gcd(n,r) <>(nicht gleich) 1) output COMPOSITE; // gcd = größter gemeinsamer Teiler
 5.  if (r is prime)
 6.    let q be the largest prime factor of r
 7.    if (q >= (4 Sqrt(r) log n)) and (n((r-1)/q) <>(nicht deckungsgleich(kongruent)) 1(mod r))
 8.      break;
 9.  r = r + 1;
10. }
11. for a = 1 to (2 Sqrt(r) log n)
12. if ( ((x-a)^n) <>(nicht deckungsgleich(kongruent)) (((x^n)-a)(mod (x^r)-1, n)) ) output COMPOSITE;
13. output PRIME;

Theorem 4.1. The algorithm above returns PRIME if and only if n is prime.

Zusätzlich ist noch das dazugehörige PDF-File angehängt.

Den GCD hab ich schon umgesetzt, aber die Zeile 1 ist mein 1. Hauptproblem.
Dabei wird festgestellt, ob n das Ergebnis einer Exponierung (a^b) ist, wobei b > 1.
Nur hab ich keinen Plan, wie ich das Umsetzen kann.

Mathematik kann so schön sein ... wenn man sie versteht !

[edit=sakura] Titel geändert. Mfg, sakura[/edit]

Da a und B element der natürlichen Zahlen sein müssen wird es etwas einfacher aber afaik gibt es keinen schnellen Algo um das zu berechnen. Da bleibt wohl nur ne For-Schleife übrig die dann durchlaufen wird.
Damit wird die Zeit für große Zahlen natürlich sehr groß.

Hast du dir das Dokument ganz genau durchgelesen ?

Tja, der schöne neue Indische Primality Test, AKS. Er ist theoretisch der schnellste Prime Test Algorithmus, WENN man die modulare Arithmetik in sehr großen Polynomen effizient programmieren KÖNNTE. Dem ist aber nicht so und somit ist es noch keinem Mathematiker gelungen tatsächlich die theoretischen Behauptungen auch praktisch zu programmieren.

Dein Problem ist momentan das die nur mit Ganzzahlen in normalen modularen Ringen arbeitest. Das ist aber grundsätzlich falsch da eben der Formelteil

((x - a)^n != (x^a - n)(mod x^r -1, n))

Polynomarithmetik mit modularen Polynomen benötigt !! Du musst dir das Dokument ganz genau durchlesen.

Vieleicht kannste ja was mit meinem damaligem Testsource was anfangen. IInteger ist ein LargInteger Typ, den du auch benötigst da die modularen Ringe der Polynome 64Bit weit überschreiten. IPoly ist das modulare Polynom, sprich sowas wie (17x^5 + 1234x^4 + 0x^3 + -123x^2 + x^1 + 1) mod p.
Gerade aber in diesem Punkt versucht ja AKS seinen Trick auszuspielen. Er arbeitet mit solchen rießigen Polynomen mit weit über 1000 Koeffizienten und das alles modular. Diese Polynome stellen sozusagen ein riesiges Set aus möglichen Teilern einer zusammengesetzten Zahl dar. Während der eigentlichen Überprüfung der Zahl auf prime wird dies also mit modularen Polynomen quasi in parallel durch sehr viele kleinerer potentielle Teiler dividiert. Sollte EINER davon ein echter Teiler sein so wird das Resultat ein Polynom sein das <> 1 ist. So kann man AKS sehr sehr vereinfacht erklären.


Gruß Hagen

zusammenfalten · markieren

Delphi-Quellcode:

(*
Input: Integer n > 1

if (n has the form ab with b > 1) then output COMPOSITE

r := 2
while (r < n) {
   if (gcd(n,r) is not 1) then output COMPOSITE
   if (r is prime greater than 2) then {
      let q be the largest factor of r-1
      if (q > 4sqrt(r)log n) and (n(r-1)/q is not 1 (mod r)) then break
   }
   r := r+1
}

for a = 1 to 2sqrt(r)log n {
   if ( (x-a)p is not (xp-a) (mod xr-1,p) ) then output COMPOSITE
}

*)

function LargestFactor(N: Cardinal): Cardinal;
var
  P: Cardinal;
  I: Integer;
begin
  for I := 1 to Primes.MaxIndex do
  begin
    P := Primes[I];
    if N mod P = 0 then
    begin
      N := N div P;
      while N mod P = 0 do
        N := N div P;
      if N = 1 then
      begin
        Result := P;
        Exit;
      end else
        if IsPrime(N) then
        begin
          Result := N;
          Exit;
        end;
    end;
  end;
end;


procedure NMix(var A: IInteger; const B: IPoly; C: Integer); overload;
// compose A as linear array of all coefficients in B where each
// chunk in A is C digits long
type
  PDigit = ^TDigit;
  TDigit = type Cardinal;

  PDigits = ^TDigits;
  TDigits = array[0..MaxInt shr 5 -1] of TDigit;

  PInt = ^TInt;
  TInt = packed record
    VTable: Pointer;
    RefCount: Integer;
    Next: Pointer;

    FCount: Cardinal;
    FSize: Cardinal;
    FNum: PDigits;
    FNeg: Boolean;
// FFlags: array[0..2] of Byte;
  end;

var
  I: Integer;
  P: PDigit;
  T: IInteger;
begin
  NSet(A, 0);
  P := (A as IDigitAccess).Alloc(C * Length(B));
  for I := 0 to High(B) do
  begin
    if NSgn(B[I]) >= 0 then NCpl(T, B[I], 32)
      else NSet(T, B[I]);
  // NMod2k1(T, B[I], 32);
    with PInt(T)^ do Move(FNum^, P^, FCount * 4);
{    if B[I] <> nil then
      with PInt(B[I])^ do Move(FNum^, P^, FCount * 4);}
    Inc(P, C);
  end;
  NNorm(A);
end;

procedure NMix(var A: IPoly; const B: IInteger; C: Integer); overload;
// decompose B of C digits chunks into each coefficent of A

  procedure MSet(var A: IInteger; P: Pointer; C: Integer);
  begin
    NSet(A, P^, C * 4);
    NCut(A, 64);
    if NSize(A) <= 48 then
    begin
      NNot(A, 32, True);
      NDec(A);
    end else NCut(A, 32);
  end;

type
  PDigit = ^TDigit;
  TDigit = type Cardinal;
var
  I,J: Integer;
  P: PDigit;
begin
  if NSgn(B) = 0 then
  begin
    A := nil;
    Exit;
  end;
  with B as IDigitAccess do
  begin
    J := Count mod C;
    I := (Count + C -1) div C;
    SetLength(A, I);
    P := Digits;
  end;
  for I := 0 to High(A) -1 do
  begin
    MSet(A[I], P, C * 4);
    Inc(P, C);
    NCalcCheck;
  end;
  I := High(A);
  MSet(A[I], P, J * 4);
  NNorm(A);
end;

function NMaxSize(const A: IPoly; Piece: TPiece = piBit): Integer; overload;
var
  S,I: Integer;
begin
  Result := 0;
  for I := 0 to High(A) do
  begin
    S := NSize(A[I], Piece);
    if S > Result then Result := S;
  end;
end;

procedure NSqrX(var A: IPoly; const B: IPoly); overload;
// A = B^2
var
  X: IInteger;
  S: Integer;
begin
  if Length(B) = 0 then
  begin
    A := nil;
    Exit;
  end;
  S := NMaxSize(B);
  Inc(S, 32 - S mod 32);
  Inc(S, NLog2(Length(B)));
  S := (S + S + 31) div 32;
  NMix(X, B, S);
  NSqr(X);
  NMix(A, X, S);
end;

procedure NSqrX(var A: IPoly); overload;
// A = A^2
begin
  NSqrX(A, A);
end;

procedure NMulX(var A: IPoly; const B,C: IPoly); overload;
// A = B * C
var
  X,Y: IInteger;
  S,T: Integer;
begin
  if (Length(B) = 0) or (Length(C) = 0) then
  begin
    A := nil;
    Exit;
  end;
  S := NMaxSize(B);
  T := NMaxSize(C);
  if T > S then S := T;
  Inc(S, 32 - S mod 32);
  Inc(S, S + NLog2(Length(B)) + NLog2(Length(C)));
  S := (S + 31) div 32;
  NMix(X, B, S);
  NMix(Y, C, S);
  NMul(X, Y);
  NMix(A, X, S);
end;

procedure NMulX(var A: IPoly; const B: IPoly); overload;
// A = A * B
begin
  NMulX(A, A, B);
end;

procedure NModxk1(var A: IPoly; K: Integer; const M: IInteger); overload;
// A := (A mod (x^k -1)) mod M
var
  Deg: Integer;
begin
  NMod(A, M);
  Deg := High(A);
  while Deg >= K do
  begin
    NAddMod(A[Deg - K], A[Deg], M);
    Dec(Deg);
    while (NSgn(A[Deg]) = 0) and (Deg >= 0) do Dec(Deg);
  end;
  SetLength(A, Deg +1);
  NNorm(A);
end;

procedure NPowModxkm1(var A: IPoly; const E: IInteger; K: Integer; const M: IInteger); overload;
// A = (A^E mod (x^k -1)) mod M
var
  I: Integer;
  T: IPoly;
  Inv2k: Cardinal;
begin
  if NSgn(E) = 0 then
  begin
    SetLength(A, 1);
    NSet(A[0], 1);
    Exit;
  end;

  NModxk1(A, K, M);
  if NCmp(E, 1) > 0 then
  begin
    NSet(T, A);
    for I := NHigh(E) -1 downto 0 do
    begin
      NSqrX(A);
      NModxk1(A, K, M);
      if NBit(E, I) then
      begin
        NMulX(A, T);
        NModxk1(A, K, M);
      end;
    end;
  end;
end;

procedure AKS;
var
  N,T: IInteger;
  R,Q,S: Cardinal;
  I,J,LogN: Integer;
  P,M,K: IPoly;
  Min,C: Double;
begin

  NSet(N, 10);
  NPow(N, 9);
  NMakePrime(N, [1,2]);

  WriteLn( NStr(N) );


  LogN := NSize(N);
  for I := 2 to Primes.MaxIndex do
  begin
    R := Primes[I];
    Q := LargestFactor(R -1);
{    if Q >= Sqrt(R) * 4 * NLn(N, 2) then
    begin
      NPowMod(T, N, NInt((R -1) div Q), NInt(R));
      if (NCmp(T, 1) <> 0) and (NCmp(T, R -1) <> 0) then Break;
    end;}

    Min := 2 * Round(Sqrt(R)) * NLn(N);
    C := 0;
    S := 0;
    while (S <= Q) and (C < Min) do
    begin
      Inc(S);
      C := C + Ln(Q + S -1) - Ln(S);
    end;
    if C >= Min then
    begin
      NPowMod(T, N, NInt((R -1) div Q), NInt(R));
      if (NCmp(T, 1) <> 0) and (NCmp(T, R -1) <> 0) then Break;
    end;
  end;
  WriteLn('R: ', R:10, ' S: ', S:10, ' Q: ', Q:10);

// Bernstein
// n=1000000007, r=3623, s=1785, q=1811

  SetLength(M, R+1); // M = x^r -1
  NSet(M[R], 1);
  NDec(N);
  NSet(M[0], N);
  NInc(N);

  C := 0;
  for I := 1 to S do
  begin
    NSet(P, [1, -I]);
    StartTimer;
    NPowModxkm1(P, N, R, N);
    C := C + Ticks;
    Write( #29, C / I / 1000 * S:10:2, #20 );
  end;
end;

Wichtig ist dieser Teil, er ist die letzte Schleife aus dem Dokument, benutzt aber eine Algorithmus der durch Daniel J. Bernstein verbesserte Eingangsparameter beechnet.

zusammenfalten · markieren

Delphi-Quellcode:

// Bernstein
// n=1000000007, r=3623, s=1785, q=1811

  SetLength(M, R+1); // M = x^r -1
  NSet(M[R], 1);
  NDec(N);
  NSet(M[0], N);
  NInc(N);

  C := 0;
  for I := 1 to S do
  begin
    NSet(P, [1, -I]);
    StartTimer;
    NPowModxkm1(P, N, R, N);
    C := C + Ticks;
    Write( #29, C / I / 1000 * S:10:2, #20 );
  end;
end;

Bei der Zahl n=1000000007 muss also r=3623, s=1785 und q=1811 sein. Dies bedeutet das in der Schleife mit dem Polynom 1000000007x^3623 + 1000000008x^0 angefangen wird. Dieses Polynom hat für die klitzekleine Zahl 1000000007
also schon 3623 Koeffizienten, baut sich also im laufe der Schleife aus 3623 eigenen Zahlen zusammen.

Wollte man zb. Zahlen > 2^256 so überprüfen so würden die Polynomberechnungen über Vektoren mit mehr als 1Mio Elementen rechnen müssen.

Nach einigen eigenen Tests und einigen Diskusssionen mit Mathematikern dachte ich das AKS ein schlechter mathematischer Witz sein muß. Wie gesagt ohne diese aufwendige Polynomarithmetik gehts nicht.Nebenbei bemerkt,auch wenn obige Codeteile ziemlich verwurschtelt aussehen so sind sie weit effizienter als die Versuche die die Leute von Miracl oder GMP mit ihren Libs angestellt haben.

Gruß Hagen

Zum Schritt "1.) if (n is of the form a^b, b>1) output COMPOSITE; // n ist das Ergebnis einer Exponierung"

Diesen kannste erstmal weglassen, wenn du die nachflogenden Schritte nach D.J. Bernstein benutzt.
Im Grunde gibt es für diesen 1. Schritt ebenfalls nur probalistische Verfahren. Meistens wird dieser Schritt auf die Überprüfung ob N ein perfektes Quadarat ist reduziert. Zusätzlich baut man aber noch eine Trialdivision zu allen kleinen Primzahlen mit ein, sprich man dividiert modular testweise N zu allen Primzahlen bis 2^16.

Es gibt aber, ebenfalls von D.J.Bernstein einige gute Dokumente für effiziente Algortihmen im Schritt 1.
Suche mal nach "Detecting perfect powers in essentially linear time" Daniel J. Bernstein, file: powers.ps.
Bernstein benutzt Operationen die auf sogenannten 2-adic Zahlen arbeiten.

Im nachfolgendem Auszug meiner Sourcen ist das dann die Funktion NIsPerfectPower(). Der Algo. von Bernstein ist in diesem Zusammenhang der schnellste bisher bekannte Weg.

Gruß Hagen

zusammenfalten · markieren

Delphi-Quellcode:

function NRootMod2k_2adic(var A: IInteger; const B,E: IInteger; K: Cardinal): Boolean;
// A = B^(1/E) mod 2^k in 2-adic metric such that (a^e * b) mod 2^k == 1 mod 2^k
var
  R,S,T,U: IInteger;

  procedure BRootMod2k_2adic(K: Cardinal);
  // recursive newton method
  begin
    if K = 1 then NSet(R, 1) else
    begin
      BRootMod2k_2adic((K +1) shr 1);

      NPowMod2k(T, R, S, K); // T = R^(E +1) * B mod 2^k
      NMulMod2k(T, B, K);
      NMulMod2k(R, S, K); // R = R * (E +1) - T
      NSub(R, T);
      NMulMod2k(R, U, K); // R = R div E mod 2^k = R * C^-1 mod 2^k
    end;
  end;

resourcestring
  sNRootMod2k = 'NRootMod2k(), Parameter must be > 0 and K > 0 and E > 0 and E = 2 or odd';
var
  I: Integer;
begin
  Result := False;
  I := NSgn(E, True);
  if (K = 0) or (I <= 0) or (NSgn(B) <= 0) then NRaise(@sNRootMod2k);
  if I = 1 then Result := NInvMod2k(A, B, K) else
    if I = 2 then Result := NSqrtMod2k_2adic(A, B, K) else
      if not NOdd(E) then NRaise(@sNRootMod2k) else
        if K = 1 then
        begin
          Result := True;
          NSet(A, 1);
        end else
        begin
          Result := NOdd(B);
          if Result then
          begin
            NInvMod2k(U, E, K); // U = E^-1 mod 2^k
            NSet(S, E);
            NInc(S); // S = E +1
            BRootMod2k_2adic(K);
            NSwp(R, A);
          end else NSet(A, 0);
        end;
end;

function NRootMod2k(var A: IInteger; const E: IInteger; K: Cardinal): Boolean;
begin
  Result := NRootMod2k(A, A, E, K);
end;

function NRootMod2k(var A: IInteger; const B,E: IInteger; K: Cardinal): Boolean;
// A = B^(1/E) mod (2^k)
begin
  if NSgn(E, True) = 1 then
  begin
    NCut(A, B, K);
    if PInt(A).FNeg then NCpl(A, K, True);
    Result := PInt(A).FCount <> 0;
  end else
    if NSgn(E, True) = 2 then Result := NSqrtMod2k(A, B, K)
      else Result := NRootMod2k_2adic(A, B, E, K) and NInvMod2k(A, K);
end;

function NIsPower(const N: IInteger; B,E: Integer): Boolean; overload;
// result := N = |B^E|
begin
  if (Abs(B) = 2) and (E >= 0) then Result := NIsPowerOf2(N) = E
    else Result := NIsPower(N, NInt(B), NInt(E));
end;

function NIsPower(const N,B: IInteger; E: Integer): Boolean; overload;
begin
  if (NCmp(B, 2, True) = 0) and (E >= 0) then Result := NIsPowerOf2(N) = E
    else Result := NIsPower(N, B, NInt(E));
end;

function NIsPower(const N,B,E: IInteger): Boolean; overload;
// result := N = |B^E|
var
  K,S: Integer;
  R,T: IInteger;
begin
  Result := False;
  K := 64;
  S := NSize(N);
  if NSgn(B) < 0 then NSet(T, B, True) else T := B;
  repeat
    if K > S then K := S;
    if not NPowMod2k(R, T, E, K) or (NCmp(R, N, K, True) <> 0) then Exit;
    if K >= S then Break;
    K := K * 8;
  until False;
  Result := True;
end;

function NIsPowerOf2(const A: IInteger): Integer;
var
  H: Cardinal;
begin
  Result := NSize(A) -1;
  if Result >= 0 then
  begin
    H := PInt(A).FNum[PInt(A).FCount -1];
    if (H and (H -1) <> 0) or (NLow(A) <> Result) then
      Result := -1;
  end;
end;

function NCheckSqrTable(R: Byte): Boolean;
const
  T256: array[0..7] of Cardinal = ($02030213,$02020212,$02020213,$02020212,
                                    $02030212,$02020212,$02020212,$02020212);
begin
  Result := Odd(T256[R div 32] shr (R mod 32));
end;

function NCheckSqr(const A: IInteger): Boolean;
const
  T193: array[0.. 6] of Cardinal = ($645AAC20,$3638B3EE,$4F85F6D4,$ADBE87C8,
                                    $DF3471B0,$10D56899,$FFFFFFFE);
  T97 : array[0.. 3] of Cardinal = ($74BAE4A0,$9F9867E4,$149D74B8,$FFFFFFFE);
  T673: array[0..21] of Cardinal = ($C05A8C20,$7A08BB4E,$6545B0DE,$3DCE2A68,
                                    $6722F3FE,$16BB3891,$266C16EA,$2DF1F0D4,
                                    $DDADF754,$9E9DA604,$07484B83,$8196E5E7,
                                    $ABBED6EC,$AC3E3ED0,$5DA0D990,$247375A1,
                                    $FF3D139A,$5951CEF1,$EC368A98,$CB744179,
                                    $10C5680D,$FFFFFFFE);
  T257: array[0.. 8] of Cardinal = ($191854E8,$81E9ABE2,$6CED7CA6,$E0897EE3,
                                    $1DFA441C,$94FADCDB,$1F565E05,$5CA86261,
                                    $FFFFFFFE);
  T241: array[0.. 7] of Cardinal = ($94EA6880,$C3985CEE,$B37056F6,$D02DE3E8,
                                    $3B345F1E,$670DBDA8,$5CA5DCE8,$FFFE0459);
  T17 = $FFFE5CE8;
  T13 = $FFFFE9E4;
  T9 = $FFFFFF6C;
  T7 = $FFFFFFE8;
  T5 = $FFFFFFEC;

  P1 = $08A19417; // 193 * 97 * 17 * 13 * 7 * 5
  P2 = $165C5F19; // 9 * 673 * 257 * 241

{$IFDEF ExpensiveSqrTest}
  T41 : array[0.. 1] of Cardinal = ($7D4AF8C8,$FFFFFE4C);
  T37 : array[0.. 1] of Cardinal = ($A1DEE164,$FFFFFFE9);
  T61 : array[0.. 1] of Cardinal = ($F5A60DC4,$E8EC196B);
  T59 : array[0.. 1] of Cardinal = ($C1846D44,$FDD49DE7);
  T53 : array[0.. 1] of Cardinal = ($CCFC512C,$FFED228F);
  T47 : array[0.. 1] of Cardinal = ($E4D8AC20,$FFFFFBCA);
  T43 : array[0.. 1] of Cardinal = ($7C5C11AC,$FFFFFCA7);
  T83 : array[0.. 2] of Cardinal = ($015CE164,$57F4EC8D,$FFFD978C);
  T79 : array[0.. 2] of Cardinal = ($7902D0C8,$BF618AAE,$FFFFECF4);
  T73 : array[0.. 2] of Cardinal = ($F472ECA0,$DD38BD86,$FFFFFE14);
  T71 : array[0.. 2] of Cardinal = ($94E26880,$E9B8D68E,$FFFFFFFE);
  T67 : array[0.. 2] of Cardinal = ($D81439AC,$A63D7E45,$FFFFFFFC);
  T251: array[0.. 7] of Cardinal = ($650C4D44,$6BE4DD27,$848471C0,$C110A94F,
                                    $E0D6AF77,$9FC71DED,$91B44D82,$FDD4DCF5);
  T239: array[0.. 7] of Cardinal = ($14A86080,$8B30CEE8,$D254F662,$48ED8F82,
                                    $D5B4BE0E,$F32EB990,$EAD7E88C,$FFFFFEF9);
  T233: array[0.. 7] of Cardinal = ($09721C68,$2A61BF8C,$E5DDEA7A,$A4CC948B,
                                    $5EEE9F44,$F6195179,$E13A40C7,$FFFFFE58);
  T229: array[0.. 7] of cardinal = ($F1E425C4,$885483CD,$37D0A72E,$9DBF6E65,
                                    $3942FB29,$F04A845D,$E909E3EC,$FFFFFFE8);
  T227: array[0.. 7] of Cardinal = ($8156E164,$35DC66E9,$E949011C,$F8EC8E05,
                                    $77F6D685,$899C453C,$978957E6,$FFFFFFFD);
  T223: array[0.. 6] of Cardinal = ($0DF03C68,$2A597508,$BDB1A8C8,$2C6699CB,
                                    $ECEA7242,$EF5165AB,$E9C3F04F);
  T211: array[0.. 6] of Cardinal = ($BCC6958C,$B20507CB,$5F602D18,$9059D546,
                                    $FB2E74BF,$CC22C1F5,$FFFCE569);
  T199: array[0.. 6] of Cardinal = ($496A9848,$18C116E4,$A1BC7EB8,$81C27A2B,
                                    $977CE7E2,$E6A96DD8,$FFFFFFED);
  T197: array[0.. 6] of Cardinal = ($C836792C,$07157049,$CAD5EFBC,$7DEAD4CC,
                                    $83AA380F,$279B04E4,$FFFFFFED);
  T191: array[0.. 5] of Cardinal = ($B0684880,$E7A0966A,$EB9C16C4,$DC97C628,
                                    $A996FA18,$FEEDE9F2);
  T181: array[0.. 5] of Cardinal = ($D5EE05C4,$A66C8309,$BF7875B4,$994B6B87,
                                    $EAE4304D,$FFE8E81D);
  T179: array[0.. 5] of Cardinal = ($55A40DC4,$C4E4136F,$5C50C3A0,$8DCFA3CF,
                                    $A550937D,$FFFDC4FD);
  T173: array[0.. 5] of Cardinal = ($5C1E19AC,$EC257481,$68C59DF6,$A90DDBEE,
                                    $1E0EA04B,$FFFFED66);
  T167: array[0.. 5] of Cardinal = ($4492A420,$18B86BAC,$9C4DC6F8,$29E2E7E0,
                                    $DAB6DDCA,$FFFFFFFB);
  T163: array[0.. 5] of Cardinal = ($F89E39AC,$88153421,$5245DB5C,$BD357EEC,
                                    $A6386E07,$FFFFFFFC);
  T157: array[0.. 4] of Cardinal = ($35F481E4,$F8E42A45,$D9B9E766,$289509C7,$E9E04BEB);
  T151: array[0.. 4] of Cardinal = ($5D80F0C8,$B379420A,$328CAACE,$45AFBD61,$FFECF0FE);
  T149: array[0.. 4] of Cardinal = ($08A4FD0C,$5F9D1B45,$7E98EDC6,$4428B62E,$FFEC2FC9);
  T139: array[0.. 4] of Cardinal = ($0CEED50C,$7D250983,$F5B41F50,$488CF3E6,$FFFFFCF5);
  T137: array[0.. 4] of Cardinal = ($ADB03468,$46F9EF0A,$DE7D88CC,$B036D543,$FFFFFE58);
  T131: array[0.. 4] of Cardinal = ($E5CE4544,$034C8521,$B5ECD3FC,$D5D8C587,$FFFFFFFD);
  T127: array[0.. 3] of Cardinal = ($399054E8,$8FE96982,$BE69680E,$E8D5F663);
  T113: array[0.. 3] of Cardinal = ($29BA1468,$0CC1EDEE,$7651DEDE,$FFFE58A1);
  T109: array[0.. 3] of Cardinal = ($418E6D44,$4FFC97A3,$9C60B17A,$FFFFE8AD);
  T107: array[0.. 3] of Cardinal = ($957681E4,$9CCC6841,$E91567DE,$FFFFFD87);
  T103: array[0.. 3] of Cardinal = ($89701C68,$4269BDA8,$C7F16EEA,$FFFFFFE9);
  T101: array[0.. 3] of Cardinal = ($3C049D8C,$FAAD57CD,$6E480F2C,$FFFFFFEC);
  T89 : array[0.. 2] of Cardinal = ($FD88F0C8,$FD594A6A,$FE4C3C46);

  T31 = $EDE2B848;
  T29 = $EC2EDD0C;
  T23 = $FFFACCA0;
  T19 = $FFFCF50C;
  T11 = $FFFFFDC4;
  T3 = $FFFFFFFC;

  I3 = $AAAAAAAB; // modular inverse of 3 to 2^32, eg 3^-1 mod 2^32
  I11 = $BA2E8BA3;
  I19 = $286BCA1B;
  I23 = $E9BD37A7;
  I29 = $4F72C235;
  I31 = $BDEF7BDF;
  I37 = $914C1BAD;
  I41 = $C18F9C19;
  I43 = $2FA0BE83;
  I47 = $677D46CF;
  I53 = $8C13521D;
  I59 = $A08AD8F3;
  I61 = $C10C9715;
  I67 = $07A44C6B;
  I71 = $E327A977;
  I73 = $C7E3F1F9;
  I79 = $613716AF;
  I83 = $2B2E43DB;
  I89 = $FA3F47E9;
  I101 = $7C32B16D;
  I103 = $D3431B57;
  I107 = $8D28AC43;
  I109 = $DA6C0965;
  I113 = $0FDBC091;
  I127 = $EFDFBF7F;
  I131 = $C9484E2B;
  I137 = $077975B9;
  I139 = $70586723;
  I149 = $8CE2CABD;
  I151 = $BF937F27;
  I157 = $2C0685B5;
  I163 = $451AB30B;
  I167 = $DB35A717;
  I173 = $0D516325;
  I179 = $D962AE7B;
  I181 = $10F8ED9D;
  I191 = $EE936F3F;
  I197 = $3D137E0D;
  I199 = $EF46C0F7;
  I211 = $6E68575B;
  I223 = $DB43BB1F;
  I227 = $9BA144CB;
  I229 = $478BBCED;
  I233 = $1FDCD759;
  I239 = $437B2E0F;
  I251 = $9A020A33;

  P3 = $6A917C99; // 41 * 37 * 31 * 29 * 23 * 19 * 3
  P4 = $FCC0D7AD; // 61 * 59 * 53 * 47 * 43 * 11
  P5 = $87B81DA9; // 83 * 79 * 73 * 71 * 67
  P6 = $BEC8D631; // 251 * 239 * 233 * 229
  P7 = $7EB0F0B1; // 227 * 223 * 211 * 199
  P8 = $48A9A435; // 197 * 191 * 181 * 179
  P9 = $2C11944D; // 173 * 167 * 163 * 157
  P10 = $19899AC9; // 151 * 149 * 139 * 137
  P11 = $0C36CCA9; // 131 * 127 * 113 * 109
  P12 = $05E7A779; // 89 * 101 * 103 * 107


  IP3 = $10BB17A9; // P3^-1 mod 2^32
  IP4 = $FC8EA425; // P4^-1 mod 2^32
  IP5 = $75B3D699; // P5^-1 mod 2^32
  IP6 = $0138C2D1; // P6^-1 mod 2^32
  IP7 = $C5775851; // P7^-1 mod 2^32
  IP8 = $BD478E1D; // P8^-1 mod 2^32
  IP9 = $701FC485; // P9^-1 mod 2^32
  IP10 = $36BB9F79; // P10^-1 mod 2^32
  IP11 = $E9489799; // P11^-1 mod 2^32
  IP12 = $4D4F12C9; // P12^-1 mod 2^32
{$ENDIF}

{prime count    probability  cum. prop.
  256  44      0.17187500  0.171875
        0.17187500000000000
  193  97      0.50259067  0.08638277202072538860103626943
  97  49      0.50515464  0.043636658031088082901554404145
  17    9      0.52941176  0.023101760134105455653764096312
  13    7      0.53846154  0.012439409302979860736642205707
    7    4      0.57142857  0.007108233887417063278081260404
    5    3      0.60000000  0.004264940332450237966848756242
        0.02481419829789229
    9    4      0.44444444  0.001895529036644550207488336108
  673  337      0.50074294  0.000949172786551580118757160874
  257  129      0.50194553  0.000476433032938341771671882306
  241  121      0.50207469  0.000239204966744976574158911863
        0.05608635715837826
  41  21      0.51219512  0.00012251961711328068432529632
  37  19      0.51351351  0.00006291547905817116222109811
  31  16      0.51612903  0.000032472505320346406307663541
  29  15      0.51724138  0.00001679612344155848602120528
  23  12      0.52173913  0.000008763194839073992706715798
  19  10      0.52631579  0.00000461220781003894352985042
    3    2      0.66666667  0.00000307480520669262901990028
        0.01285426991142190
  61  31      0.50819672  0.000001562605924712647534703421
  59  30      0.50847458  0.000000794545385447108915950892
  53  27      0.50943396  0.00000040476840390701774963536
  47  24      0.51063830  0.000000206690248803583531728695
  43  22      0.51162791  0.000000105748499387879946465844
  11    6      0.54545455  0.000000057680999666116334435915
        0.01875923702111852
  83  42      0.50602410  0.000000029187975734661277666366
  79  40      0.50632911  0.000000014778721890967735527274
  73  37      0.50684932  0.000000007490585068024742664509
  71  36      0.50704225  0.000000003798043133082968111582
  67  34      0.50746269  0.000000001927365172012252474534
        0.03341421236054701
  251  126      0.50199203  0.000000000967521958858740286021
  239  120      0.50209205  0.00000000048578508394581102227
  233  117      0.50214592  0.000000000243934999234591800882
  229  115      0.50218341  0.000000000122500108785930380356
        0.06355832852268207
  227  114      0.50220264  0.000000000061519878421128032425
  223  112      0.50224215  0.000000000030897876157696590276
  211  106      0.50236967  0.000000000015522155794861794167
  199  100      0.50251256  0.000000000007800078288875273451
        0.06367405193497382
  197  99      0.50253807  0.000000000003919836297455086658
  191  96      0.50261780  0.00000000000197017950029156188
  181  91      0.50276243  0.000000000000990532234953216194
  179  90      0.50279330  0.00000000000049803296729491317
        0.06384973956033545
  173  87      0.50289017  0.000000000000250455885287037259
  167  84      0.50299401  0.000000000000125977810563539699
  163  82      0.50306748  0.000000000000063375340283498499
  157  79      0.50318471  0.000000000000031889502435645742
        0.06403090664631079
  151  76      0.50331126  0.00000000000001605034559674885
  149  75      0.50335570  0.000000000000008079033018497743
  139  70      0.50359712  0.000000000000004068577779099583
  137  69      0.50364964  0.000000000000002049137713561104
        0.06425743762218754
  131  66      0.50381679  0.000000000000001032389993091854
  127  64      0.50393701  0.000000000000000520259524077785
  113  57      0.50442478  0.000000000000000262431795331272
  109  55      0.50458716  0.00000000000000013241971324055
        0.06462216392983359
  107  54      0.50467290  0.000000000000000066828640327007
  103  52      0.50485437  0.000000000000000033738731038877
  101  51      0.50495050  0.000000000000000017036388940423
  89  45      0.50561798  0.000000000000000008613904520439
        0.06505001641855890

1036 Bytes needed for all Lookuptable
  252 Bytes needed for no expensive test

  Follow code examine in avg. ~2.5 Cycles per Digit if A is NOT a perfect Sqr.
  on 2048 Bit A ~ 3000 times faster as NSqrt() }

var
  R,T: Cardinal;
  S: Long96;
begin
  Result := (A = nil) or (PInt(A).FCount = 0);
  if not Result and not PInt(A).FNeg and NCheckSqrTable(A[0]) then
  begin
// prop: 0.171875 pass above
// ~12 Cycles
    S := NSum(A); // S = A mod 2^96-1 // A.Count x Additions + 1 DIV
    R := NMod(S, P1);
    T := R mod 193; if Odd(T193[T shr 5] shr (T and 31)) then Exit;
    T := R mod 97; if Odd(T97 [T shr 5] shr (T and 31)) then Exit;
    if Odd(T17 shr (R mod 17)) then Exit;
    if Odd(T13 shr (R mod 13)) then Exit;
    if Odd(T7 shr (R mod 7)) then Exit;
    if Odd(T5 shr (R mod 5)) then Exit;
    NCalcCheck;

// ~ 4 Cycles
    R := NMod(S, P2); // 3 x 32Bit DIV's
    if Odd(T9 shr (R mod 9)) then Exit;
    T := R mod 673; if Odd(T673[T shr 5] shr (T and 31)) then Exit;
    T := R mod 257; if Odd(T257[T shr 5] shr (T and 31)) then Exit;
    T := R mod 241; if Odd(T241[T shr 5] shr (T and 31)) then Exit;
// prop: 0.0002392 pass above, 1 of 4181 Candidates pass

{$IFDEF ExpensiveSqrTest}
    R := DModDInv2k32(PInt(A).FNum, P3, PInt(A).FCount, IP3);
    T := DInvMul2k32(R, 41, I41); if Odd(T41[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 37, I37); if Odd(T37[T shr 5] shr (T and 31)) then Exit;
    if Odd(T31 shr DInvMul2k32(R, 31, I31)) then Exit;
    if Odd(T29 shr DInvMul2k32(R, 29, I29)) then Exit;
    if Odd(T23 shr DInvMul2k32(R, 23, I23)) then Exit;
    if Odd(T19 shr DInvMul2k32(R, 19, I19)) then Exit;
    if Odd(T3 shr DInvMul2k32(R, 3, I3)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P4, PInt(A).FCount, IP4);
    T := DInvMul2k32(R, 61, I61); if Odd(T61[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 59, I59); if Odd(T59[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 53, I53); if Odd(T53[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 47, I47); if Odd(T47[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 43, I43); if Odd(T43[T shr 5] shr (T and 31)) then Exit;
    if Odd(T11 shr DInvMul2k32(R, 11, I11)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P5, PInt(A).FCount, IP5);
    T := DInvMul2k32(R, 83, I83); if Odd(T83[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 79, I79); if Odd(T79[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 73, I73); if Odd(T73[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 71, I71); if Odd(T71[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 67, I67); if Odd(T67[T shr 5] shr (T and 31)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P6, PInt(A).FCount, IP6);
    T := DInvMul2k32(R, 251, I251); if Odd(T251[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 239, I239); if Odd(T239[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 233, I233); if Odd(T233[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 229, I229); if Odd(T229[T shr 5] shr (T and 31)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P7, PInt(A).FCount, IP7);
    T := DInvMul2k32(R, 227, I227); if Odd(T227[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 223, I223); if Odd(T223[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 211, I211); if Odd(T211[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 199, I199); if Odd(T199[T shr 5] shr (T and 31)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P8, PInt(A).FCount, IP8);
    T := DInvMul2k32(R, 197, I197); if Odd(T197[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 191, I191); if Odd(T191[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 181, I181); if Odd(T181[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 179, I179); if Odd(T179[T shr 5] shr (T and 31)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P9, PInt(A).FCount, IP9);
    T := DInvMul2k32(R, 173, I173); if Odd(T173[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 167, I167); if Odd(T167[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 163, I163); if Odd(T163[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 157, I157); if Odd(T157[T shr 5] shr (T and 31)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P10, PInt(A).FCount, IP10);
    T := DInvMul2k32(R, 151, I151); if Odd(T151[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 149, I149); if Odd(T149[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 139, I139); if Odd(T139[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 137, I137); if Odd(T137[T shr 5] shr (T and 31)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P11, PInt(A).FCount, IP11);
    T := DInvMul2k32(R, 131, I131); if Odd(T131[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 127, I127); if Odd(T127[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 113, I113); if Odd(T113[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 109, I109); if Odd(T109[T shr 5] shr (T and 31)) then Exit;

    R := DModDInv2k32(PInt(A).FNum, P12, PInt(A).FCount, IP12);
    T := DInvMul2k32(R, 107, I107); if Odd(T107[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 103, I103); if Odd(T103[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 101, I101); if Odd(T101[T shr 5] shr (T and 31)) then Exit;
    T := DInvMul2k32(R, 89, I89 ); if Odd(T89 [T shr 5] shr (T and 31)) then Exit;
// 1 of 116091372690195275 candidates come here
{$ENDIF}
    Result := True;
  end;
end;

function NCheckSqrt(const A: IInteger): Boolean;
asm // check if A is a perfect Square, by doing a NSqrt(nil, A)
       MOV EDX,EAX
       XOR EAX,EAX
       JMP NSqrt // Result := NSqrt(nil, A);
end;

function NIsPerfectSqr(const A: IInteger; FullTest: Boolean): Boolean;
// test if A is a perfect square
begin
  Result := NCheckSqr(A) and (not FullTest or NCheckSqrt(A));
end;

function NIsPerfectSqr(const N: Int64): Boolean;
var
  R: Cardinal;
begin
  Result := N = 0;
  if not Result and (N > 0) and NCheckSqrTable(N) then
    Result := (DSqrt64(@R, @R, @N, 2) = 0) and (R = 0);
end;

function NIsPerfectPower(var B: IInteger; const N: IInteger; Bound: Cardinal = 0): Cardinal;
// check if N is a perfect power, if true returns Base in B and Exponent (smallprime) in Result
// implementation of
// "Detecting perfect powers in essentially linear time"
// Daniel J. Bernstein, file: powers.ps
var
  X,Y: IInteger;
  I,K,L: Integer;
  P: TSmallPrimeSieve;
label
  Error;
begin
  Result := 0;
  if Bound = 0 then Bound := TSmallPrimeSieve.MaxPrime;
  I := NSgn(N, True);
  if I < 0 then goto Error;
  if I < 3 then
  begin
    if @B <> nil then NSet(B, I);
    Result := 1;
    Exit;
  end;
  if NCheckSqr(N) then
    if @B <> nil then
    begin
      if NSqrt(X, N) then
      begin
        NSwp(B, X);
        Result := 2;
        Exit;
      end;
    end else
      if NCheckSqrt(N) then
      begin
        Result := 2;
        Exit;
      end;

  if not NOdd(N) then
  begin
    // to slow,
{    I := 2;
    L := Primes[I];
    while L < K do
    begin
      if NRoot(X, N, L) then
      begin
        if @B <> nil then NSwp(B, X);
        Result := L;
        Exit;
      end;
      if L >= Bound then goto Error;
      Inc(I);
      L := Primes[I];
    end; }
    goto Error;
  end else
    if Bound > 2 then
    begin
      K := NSize(N);
      L := (K +1) shr 1;
      if NInvMod2k(Y, N, L +1) then // Y = N^-1 mod 2^(Round(K/2) +1)
      begin
      // first check to exponent 2,
      // suppressed: NIsPerfectSqr() is always faster
   {    if NSqrtMod2k_2adic(X, Y, L) then  // X = Y^1/2 mod 2^(Round(K/2)
        begin
          if NIsPower(N, X, 2) then
          begin
            if @B <> nil then NSwp(B, X);
            Result := 2;
            Exit;
          end;
          NCpl(X, L);                      // X = 2^Round(K/2)  - X
          if NIsPower(N, X, 2) then
          begin
            if @B <> nil then NSwp(B, X);
            Result := 2;
            Exit;
          end;
        end;}

     // now check to smallprime powers
        P := Primes;
        NAutoRelease(P);
        I := 2;
        L := P[I];
        while L < K do
        begin
          if NRootMod2k_2adic(X, Y, NInt(L), (K + L -1) div L) and NIsPower(N, X, L) then
          begin
            if @B <> nil then NSwp(B, X);
            Result := L;
            Exit;
          end;
          if Cardinal(L) >= Bound then goto Error;
          Inc(I);
          L := P[I];
        end;
        if @B <> nil then NSet(B, N);
        Result := 1;
      end else
      begin
     Error:
        if @B <> nil then NSet(B, 0);
        Result := 0;
      end;
    end;
end;