Miller-Rabin primzahltest

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