Algorithmus zur Bestimmung der Primalität einer Ganzahl

[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;