Miller-Rabin primzahltest

[qwertz543221]

Delphi-Quellcode:

function witness(a,n:ansistring):boolean;
var f,x,i,t,u:ansistring;
mathe:tmathe;
begin
result:=false;
mathe:=tmathe.Create;
t:='0';
u:=mathe.Differenz(n,'1');

while (mathe.Modulo(u,'2')='0') do
begin
u:=mathe.Quotient(u,'2');
t:=mathe.Summe(t,'1');
end;

x:=mathe.PotenzModulo(a,u,n);

i:='1';
while mathe.Vergleich(i,t,vkleiner)do
begin
f:=x;
x:=mathe.PotenzModulo(f,'2',n);

if (x='1')and(f<>'1')and(f<>mathe.Differenz(n,'1'))
  then
  begin
  result:=true;
  exit;
  end;
i:=mathe.Summe(i,'1');
end;

if not result
  then result:=x<>'1';
end;



function millerrabin(n,s:ansistring):boolean;
var i,a:ansistring;
mathe:tmathe;
begin
result:=true;
i:='1';
mathe:=tmathe.Create;
if mathe.istGerade(n)
then result:=true
else
if (mathe.istungerade(n))and(mathe.Vergleich(n,'3')>0)
  then
  begin
while mathe.Vergleich(i,s,vkleiner) do
begin
a:=mathe.Zufall(mathe,'0',mathe.Differenz(n,'1'));

if witness(a,n)=false
  then begin
  result:=false;
  exit;
  end
  else begin
  result:=true;
  exit;
  end;
end;
end
else
begin result:=false;
exit;
end;

end;

mithilfe der kleinen mathe lib(http://www.delphipraxis.net/internal...t.php?t=159592) habe ich den code von miller und rabin zur bestätigung der nicht-primalität implementiert. Das ganze liefert jedoch falsche testergebnisse, und ich weiß nicht wo sich - wiedermal- ein denkfehler versteckt.

beispielezahl/durchläufe/ergebnis)

2/50/NP
3/50/P
3/500/P
17/50/P
17/395/P
17/30/NP
991/300/NP
991/3000/P
991/20/P
991/50/NP

Oder liegt das Ganze an der im miller-rabin-algorithmus inbegriffenen fehlerwahrscheinlichkeit??

[gammatester]

Nun, es fehlen mal wieder jegliche Hinweise und Kommentare, was was ist. Liefert witness true wenn a ein Zeuge für die Nichtprimalität von n ist (das wäre die Standarddefinition)? Gleiche Frage für millerrabin.

Auf jeden Fall ist die Schleife

Delphi-Quellcode:

for i:=1 to s-1 do begin
  a:=mathe.Zufall(mathe,'0',mathe.Differenz(n,'1')); 
  if witness(a,n)=false then begin
    result:=false;
    exit;
  end
  else begin
    result:=true;
    exit;
  end;
end;

Blödsinn da sie nach dem ersten witness immer verlassen wird. (Bem: Deine boolean-Konstruktionen und Einrückungen sind auch nicht besonders übersichtlich)

Delphi-Quellcode:

result := witness(a,n);
exit;

tut genau das gleiche.


Soviel zur falschen Logik. Ein weiterer Hinweis: Normalerweise würde ich am Ende von witness erwarten, daß x mit n-1 vergleichen wird und nicht x mit 1. Also etwa:

if f<>mathe.Differenz(n,'1')) then result := (nicht prim).

Wobei für (nicht prim) steht für Deine Wahl was witness eigentlich zurückliefert.


Gammatester

PS: Quadrieren sollte man auch statt mit
x :=mathe.PotenzModulo(f,'2',n);
besser und schneller mit
x := ProduktModulo(f, f, n);

[gammatester]

Zitat von qwertz543221:

mithilfe der kleinen mathe lib(http://www.delphipraxis.net/internal...t.php?t=159592) habe ich den code von miller und rabin zur bestätigung der nicht-primalität implementiert. Das ganze liefert jedoch falsche testergebnisse, und ich weiß nicht wo sich - wiedermal- ein denkfehler versteckt.

beispiele:(zahl/durchläufe/ergebnis)

2/50/NP
3/50/P
3/500/P
17/50/P
17/395/P
17/30/NP
991/3000/P

Oder liegt das Ganze an der im miller-rabin-algorithmus inbegriffenen fehlerwahrscheinlichkeit??

Diese Zahlen zeigen, daß Du noch nicht ganz verstanden hast, was der Test eigentlich macht: In witness(a, n) ist a eine Zufallszahl mit 2 <= a <= n-2. Nun überleg mal wieviel (Zufalls-)Zahlen es zwischen 2 und 0 gibt, zwischen 2 und 1, und zwischen 2 und 17 etc.

Wenn Du Probleme hast Code von C etc nach Pascal umzusetzen, hilft Dir vielleicht meine Miller-Rabin-Implementation aus MPArith weiter

Delphi-Quellcode:

{---------------------------------------------------------------------------}
procedure mp_miller_rabin(const n: mp_int; t: word; var prime: boolean);
  {-Miller-Rabin test of n, security parameter t, from HAC p. 139 Alg.4.24}
  { if t=0, calculate t from bit_count with prob of failure < 2^-96}
label
  leave;
var
  n1, y, r: mp_int;
  s,j: longint;
const
  t0: array[0..14] of byte = (50,28,16,10,7,6,5,4,4,3,3,3,3,3,2);
begin
  {default}
  prime := false;

  if mp_error<>MP_OKAY then exit;
  {$ifdef MPC_ArgCheck}
    if mp_not_init(n) then begin
      {$ifdef MPC_HaltOnArgCheck}
        {$ifdef MPC_UseExceptions}
          raise MPXNotInit.Create('mp_miller_rabin');
        {$else}
          RunError(MP_RTE_NOTINIT);
        {$endif}
      {$else}
        set_mp_error(MP_NOTINIT);
        exit;
      {$endif}
    end;
  {$endif}

  {easy case n < 2^31 or even}
  if mp_is_longint(n,s) then begin
    if s>1 then prime := IsPrime32(s);
    exit;
  end;
  if mp_iseven(n) or (n.sign=MP_NEG) then exit;

  mp_init3(n1,r,y);
  if mp_error<>MP_OKAY then exit;

  {get n1 = n - 1}
  mp_sub_d(n,1,n1);

  {calculate r,s with n-1=2^s*r}
  mp_makeodd(n1,r,s);

  if t=0 then begin
    {calculate t from bit_count with prob of failure lower than 2^-96}
    j := mp_bitsize(n) div 128;
    if j>14 then t:=14 else t:= word(j);
    t := t0[t];
  end;

  while t>0 do begin
    {generate a in the range 2<=a<n-1, calculate y = a^r mod n}
    repeat
      mp_rand(y, n.used);
      mp_mod(y,n1,y);
      if mp_error<>MP_OKAY then goto leave;
    until mp_cmp_d(y,1)=MP_GT;
    mp_exptmod(y, r, n, y);
    if mp_error<>MP_OKAY then goto leave;

    {if y<>1 and y<>n-1 do}
    if (not mp_is1(y)) and mp_is_ne(y, n1) then begin
      j := 1;
      while (j <= s-1) and mp_is_ne(y, n1) do begin
        mp_sqrmod(y, n, y);
        {if y=1 then composite}
        if mp_is1(y) or (mp_error<>MP_OKAY) then goto leave;
        inc(j);
      end;
      {if y<>n1 then composite}
      if (mp_is_ne(y, n1)) or (mp_error<>MP_OKAY) then goto leave;
    end;
    dec(t);
  end;
  {probably prime now}
  prime := true;

leave:
  mp_clear3(n1,r,y);
end;

[himitsu]

also irgendwie komm ich mit den ganzen Pseudocodes auch nicht klar

z.B.: http://en.wikipedia.org/wiki/Miller-...primality_test

nja, so weit bin ich grad mal "schnell" gekommen ...

Delphi-Quellcode:

Type TMillerRabinResult = (mrZusammengesetzt, mrWahrscheinlichPrimzahl);

Function MillerSelfridgeRabinTest(n, k: String): TMillerRabinResult;
  Var n1, n2, s, d, a, x, r: String;

  Begin
    n := Mathe.Normalisieren(n);
    k := Mathe.Normalisieren(k);
    Result := mrZusammengesetzt;
    If Mathe.istGerade(n) or (Mathe.Vergleich(n, '2') <= 0) Then Exit;
    n1 := n;
    Mathe.Minus1(n1);
    n2 := n1;
    Mathe.Minus1(n2);
    //write n − 1 as 2^s*d with d odd by factoring powers of 2 from n − 1
    While Mathe.Vergleich(k, '0') > 0 do Begin
      Mathe.Minus1(k);
      a := Mathe.Zufallszahl('2', n2);
      x := Mathe.PotenzModulo(a, d, n);
      If (x <> '1') and (x <> n1) Then Begin
        r := '1';
        If Mathe.Vergleich(r, s) < 0 Then
          While True do Begin
            Mathe.Plus1(r);
            x := Mathe.PotenzModulo(x, '2', n);
            If (x = '1') or (Mathe.Vergleich(r, s) >= 0) Then Exit;
            If x = n1 Then Break;
          End;
      End;
    End;
    Result := mrWahrscheinlichPrimzahl;
  End;

und beim Anderem
http://www.jjam.de/Java/Applets/Prim...ler_Rabin.html

hatte ich da erstmal aufgegeben

Delphi-Quellcode:

Function MillerSelfridgeRabinTest(n, s: String): TMillerRabinResult;
  Function Witness(a, n: String): Boolean;
    Begin
      //If Mathe.Vergleich(Mathe.Differenz(n, '1'), ) = 0 Then
      let n-1 = 2^t*u, where t>=1 and u is odd
      x[0] := modular_exponentiation(a,u,n);
      for i := 1 to t do Begin
          x[i] := x[i-1]^2 mod n
          if x[i] = 1 and x[i-1] <> 1 and x[i-1] <> n-1 then return true;
      End;
      Result := x[t] <> 1;
    End;

  Var n2: String;

  Begin
    Result := mrZusammengesetzt;
    If Mathe.istGerade(n) or Mathe.istNegativ(n) Then Exit;
    Result := mrWahrscheinlichPrimzahl;
    n2 := n;
    n2 := Mathe.Minus1(n2);
    While Mathe.Vergleich(s, '0') > 0 do Begin
      Mathe.Minus1(s);
      If Witness(Mathe.Zufallszahl('1', n2), n) Then Exit;
    End;
    Result := mrZusammengesetzt;
  End;

ich kappier bei einigen "Befehlen" einfach nicht, was die da wollen

[gammatester]

Also die schnelle direkte Übertragung meiner Routine nach StringMatheLib scheint doch zu funktionieren:

Delphi-Quellcode:

{---------------------------------------------------------------------------}
procedure miller_rabin(n: string; t: longint; var prime: boolean);
  {-Miller-Rabin test of n, security parameter t, from HAC p. 139 Alg.4.24}
var
  n1,n2, y, r, x: string;
  s,j: longint;
Begin

  n := Mathe.Normalisieren(n);
  if (n='2') or (n='3')  or (n='5') then begin
    prime := true;
    exit;
  end;

  prime := false;
  if Mathe.istGerade(n) or (Mathe.Vergleich(n, '6') <= 0) then exit;

  {get n1 = n - 1}
  {get n2 = n - 2}
  n1 := n;
  Mathe.Minus1(n1);
  n2 := n1;
  Mathe.Minus1(n2);

  {calculate r,s with n-1=2^s*r}
  r := n1;
  s := 0;
  while Mathe.istGerade(r) do begin
    r := Mathe.Quotient(r,'2');
    inc(s);
  end;

  while t>0 do begin
    {generate a in the range 2<=a<n-1, calculate y = a^r mod n}
    y := Mathe.Zufallszahl('2', n2);
    y := Mathe.PotenzModulo(y,r,n);
    {if y<>1 and y<>n-1 do}
    if (y<>'1') and (y<>n1) then begin
      j := 1;
      while (j <= s-1) and (y<>n1) do begin
        y := Mathe.ProduktModulo(y,y,n);
        {if y=1 then composite}
        if y='1' then exit;
        inc(j);
      end;
      {if y<>n1 then composite}
      if y<>n1 then exit;
    end;
    dec(t);
  end;
  {probably prime now}
  prime := true;
end;

{---------------------------------------------------------------------------}
procedure TForm1.Button1Click(Sender: TObject);
  {-Testcode für miller_rabin}
var
  prime: boolean;
  n: string;
begin
  n := Mathe.Normalisieren(Edit1.Text);
  Edit1.Text := n;
  miller_rabin(n,3,prime);
  if prime then button1.Caption := 'Prime' else button1.Caption := 'NOT prime'
end;

Ein paar kleine "Einschränkungen", die aber bei StringMatheLib garantiert nicht zum Tragen kommen: der Sicherheitsparameter t ist im Longint-Bereich, außerdem muss s von n-1 = 2^s*r im Longint-Bereich sein.

In der Praxis benutzt man allerdings Miller-Rabin nicht für kleine Zahlen, ohne vorher Probedivisionen mit einigen kleine Primzahlen zu machen, mindestens bis 100 besser zB alle 16Bit-Primzahlen.

Gammatester

[qwertz543221]

danke an euch , das hat super funktioniert