Algorithmus zur Bestimmung der Primalität einer Ganzahl

[AlphaBug]

Hallo Leute ,

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

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]

[atreju2oo0]

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ß.

[negaH]

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

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;

[negaH]

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

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

[negaH]

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

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;