Miller-Rabin - Eigene Impl. findet nicht jede Primzahl...

[St.Pauli]

Problem bei meiner Implentierung des Miller-Rabin Tests

Wie schon im Titel erwähnt, findet meine Implentierung von dem Miller-Rabin Primzahltest leider nicht jede Primzahl.
Das liegt leider an meiner Implentierung des Algos, ich kann aber absolut nicht den Fehler finden.

Delphi-Quellcode:

function Bin(Int: Integer): String;
VAR i : Integer;
begin
  Result := '';
  for i := 7 downto 0 do
    Result := Result + IntToStr((Int shr i) and 1);
end;

function modular_exponentiation(a,b,n : integer) : integer;
VAR c, d, i : integer;
    b_bin : string;
begin
    c := 0;
    d := 1;
    b_bin := Bin(b);
    for i := Length(b_bin) downto 1 do
      begin
        c := 2*c;
        d := (d*d) mod n;
        IF (b_bin[i] = '1') THEN
          begin
            d := (d*a) mod n;
            Inc(c);
          end;
      end;
    Result := d;
end;

function witness(a, n : integer) : boolean;
VAR f, x, i, t, u : integer;
begin
  Result := False;

  // n-1 = 2^t * u

  t := 0;
  u := n-1;

  while ((u mod 2) = 0) do
    begin
      u := Round(u/2);
      Inc(t);
    end;

  x := modular_exponentiation(a, u, n);

  for i := t downto 1 do
    begin
      f := x;
      x := (f*f) mod n;

      IF (x = 1) AND (f <> 1) AND (f <> n-1) THEN
        begin
          Result := True;
        end;

    end;

  IF NOT Result THEN
    Result := (x <> 1);
end;

function Miller_Rabin(n, s : integer) : boolean;
VAR i : integer;
begin
  Result := True;

  Randomize;

  for i := 1 to s do
    begin
      IF (witness((Random(n-1) + 1), n)) THEN
        Result := False;
    end;

end;

Bei der Implentierung habe ich mich an diesen Pseudocode gehalten.
Ich freue mich jetzt schon über jede Antwort... THX. :cheers:

[Der_Unwissende]

Hi,
ich hab mir deinen Code + Pseudocode noch nicht angeschaut, aber mal als Mögliche Antwort vorweg (hoffe ich verwechsel hier nicht Miller-Rabin mit einem ganz anderen), es handelt sich um einen Pseudoprimzahltest. Du versuchst doch (und jetzt weißt du gleich ob ich mich irre) Zeugen für die Zerlegbarkeit (in Primfaktoren) deiner zu testenden Zahl zu finden. Hm, Basis wird zufällig gewählt und dann, weiß ich gar nicht mehr :wink:

Jedenfalls gibt es eine Fehlerrate, die von dieser Zufälligkeit (der Basis) abhängt oder anders gesagt, die Ermittlung deiner Zeugen (wie auch immer das da war), ist hier der einzigste Kritische Punkt. Der Algorithmus ist so gut, weil er recht schnell approximiert (aber eben nicht mehr). Du wirst also Fehler haben, aber die Wahrscheinlichkeit ist sehr sehr gering.

Da bei dir aber trotzdem welche auftreten (genug um bemerkt zu werden), schau ich gleich mal.
Falls ich von einem falschen Algo rede, sorry, dann vergiß alles was ich gesagt habe.

Gruß Der Unwissende

[St.Pauli]

OK, schon mal danke. Ja, du redest von dem richtigen Algorithmus!!!

Zitat:

Zitat von Der_Unwissende

Du wirst also Fehler haben, aber die Wahrscheinlichkeit ist sehr sehr gering.

Ja - Das mit der Fehlerwahrscheinlichkeit stimmt, sie ist >= 1/2 ^ s (glaub ich). Bei mir in der Praxis beträgt dies bei 20 Wiederholungen 0.00000095367431640625, dazu noch die Tatsache, dass immer die gleichen Zahlen als nicht-prim angezeigt werden und es ist fast unmöglich, dass es sich um eine Kettung von Zufällen handeln kann.

[Der_Unwissende]

Ok, ich muss sagen ich seh auch keinen Fehler.
Bei welcher Zahl klappt es denn immer nicht?

Kann mir nur vorstellen, das die Rundung beim Finden des Zeugen irgendeinen Nebeneffekt hat (geliebte Floating Points).
Aber du kannst den Algorithmus etwas umschreiben :

Delphi-Quellcode:

  t := 0; 
  u := n-1; 

  while ((u mod 2) = 0) do 
    begin 
      u := u shr 1; 
      Inc(t); 
    end;

Sollte zwar nicht wirklich etwas verändern, aber gut. Hatte erst überlegt, dass du so natürlich das Aufrunden verlieren würdest, aber ich bin mir halt nicht sicher, ob da schon ein Fehler drin steckt.

Die Idee ist es doch, dass n-1 eine gerade Zahl ist. t >= 1, damit 2^{t} >= 2 und damit wiederum bin(n-1) = bin(u) shl t, mit bin (x) = Binärdarstellung von x. Also musst du nur die Binärdarstellung von (n-1) nehmen und die 0en bis zur ersten eins zählen. Und eigentlich würde ich sagen machst du genau das (und dir war sicherlich auch die Idee klar, wollte es nur sagen um meinen Gedankengang nachvollziehbarer zu machen).
Wenn ich immer um eine Stelle nach Rechts verschiebe, hm, runde ich doch immer ab? Bei Round hingegen würde ich teilweise auch abrunden. Allerdings dürfte es (da nur solange geshiftet wird, wie die letzte Stelle 0 ist) gar keine Rundung geben. Also ist es nur aus Performancegründen sinnvoll zu shiften (anders gesagt es ist völlig für den A...).

Hm, ne, dann seh ich den Fehler erst recht nicht. Keine Ahnung warum es nicht macht was es soll. Vielleicht doch jmd. anderes?

Gruß Der Unwissende

[Amateurprofi]

Hallo St.Pauli,

im Pseudo Code steht

Delphi-Quellcode:

a:=Random(1, n-1)
if Witness(a,n) then

Die dort benutze Random Version, so denke ich, gibt eine Zufallszahl im Bereich 1 .. n-1 aus.
Der Parameter a ist also minimal 1, maximal n-1

In Deiner Version wolltest du das a:= .... sparen und ersetzt den Parameter a durch

Delphi-Quellcode:

Random(n-2) + 1

Random(n-2) ergibt eine Zufallszahl im Bereich 0 .. n-3.
Hierzu addierst du 1, erhältst also einen Wert im Bereich 1..n-2
Dein Parameter a ist minimal 1, maximal n-2

Vielleicht liegts daran?!

[St.Pauli]

Zitat:

Zitat von Der_Unwissende

Ok, ich muss sagen ich seh auch keinen Fehler.
Bei welcher Zahl klappt es denn immer nicht?

23, 47, 53, 59, 71, 79, 83, 89 zwischen 1 und 100. Erkannt hat er 5, 7, 11, 13, 17, 19, 29, 31, 37, ... Also schon mal mehr erkannt als nicht :lol: Und jeh mehr ich damit rumspiele und teste, wie die Zwischenwerte für erkannte und nicht erkannte Primzahlen aussehen, desto mehr denke ich, den Fehler nie zu finden. :|

Delphi-Quellcode:

  t := 0;
  u := n-1;

  while ((u mod 2) = 0) do
    begin
      u := u shr 1;
      Inc(t);
    end;

Bringt leider auch, wie du ja schon gesagt hast, keine Veränderung.

Edit: Warum wurde mir kein roter Kasten angezeigt, für einen neuen Post??? :gruebel:

Zitat:

Zitat von Amateurprofi

Random(n-2) ergibt eine Zufallszahl im Bereich 0 .. n-3.
Hierzu addierst du 1, erhältst also einen Wert im Bereich 1..n-2
Dein Parameter a ist minimal 1, maximal n-2

Vielleicht liegts daran?!

Danke für den Hinweis! Das ist in der Tat ein Denkfehler von mir gewesen :wall: Leider funktioniert es jetzt immer noch nicht, da ja nur durch meinen Fehler die Wahrscheinlichkeit eine zusammengesetzte Zahl als eine Primzahl auszugeben erhöt wurde!!!

[Amateurprofi]

Naja St.Pauli, dann hab ich mal deine Hausaufgaben gemacht.

In der Funktion Witness heißt es

Zitat:

for i <-- 1 to t
do blablabla
if blablabal
then return true

Ich kenne mich mit Pseudocodes nicht aus, interpretiere das aber so, daß wenn "if blablabla" True ergibt, True zurückgegeben werden soll. D.h. es müßte result=true gesetzt und die Funktion verlassen werden. Du setzt zwar result=true verläßt aber nicht die Funktion. Ich hab das mal probiert, brachte aber auch nichts. Herr MillerRabin sagte immer noch 23 sei composite.

Also weitergesucht.

In der Funktion ModularExponentiation heit es im Pseudo Code

Zitat:

Let (bk, bk-1,...b0) be the binary representation of b
for i <-- k downto 0

Bei Dir wird das so umgesetzt daß, b als Binärstring in b_bin gestellt wird.
b_bin[Length(b_bin)] ist dabei Bit0.
Auf den ersten Blick sieht das so aus, wie im Pseudo Code, aber beim zweiten Blick eben nicht.
Die Schleife for i:=k downto 0 faßt zuerst hast höchstwertige Bit an und zuletzt Bit0.
Bei Deiner Schleife for i:=Length(b_bin) downto 1 ist es genau umgekehrt.
Ich hab das korrigiert und dann für alle (ungeraden) Zahlen von 3 bis 65535 probiert. Bis 46337 bringt das Programm korrekte Ergebnisse. Für höhere Werte sagt es nur noch composite.... Hab nicht weiter geprüft, warum.

Und weil Du Dich für Primzahlen zu interessieren scheinst ist hier noch eine kleine Routine,
die Extendend-Werte daraufhin prüft, ob sie prim sind. Die zu prüfende Zahl darf natürlich nicht größer als 2^64-1 sein, also max 64 Bit Werte.

Delphi-Quellcode:

{------------------------------------------------------------------------------}
{ FIsPrimeNPA                                                                  }
{ Returns true if p is prime. NO ABORT-FUNCTIONALITY, NO PROGRESSBAR           }
{ Argument : p = the number to test.                                           }
{ Result   : 0 = p is not prime  -  then the smallest divisor is returned in r }
{            1 = p is prime                                                    }
{ p can be any extended number.                                                }
{ the routine tests if it is a positive integer and returns false if not.      }
{------------------------------------------------------------------------------}
FUNCTION FIsPrimeNPA(p:extended; var r:extended):integer;
asm
         mov      edx,eax            // save pointer to r
         sub      esp,12             // reserve space on stack
         fnstcw   [ESP].Word         // save control word
         fnstcw   [ESP+2].Word       // and once more
         fwait
         or       [ESP+2].Word,$0F00 // trunc decimals, full precision
         fldcw    [ESP+2].Word      // load new control word
         fldz                       // ST0=0
         fstp     tbyte [edx]       // as default store 0 as largest divisor
         fldz                       // ST0=0
         fld      p                 // ST0=p  ST1=0
         fld1                       // ST0=1  ST1=p  ST2=0
         fadd     st,st             // ST0=2  ST1=p  ST2=0
         fcom     st(1)             // compare 2 with p
         fstsw    ax                // status word into AX
         sahf                       // copy ah into flag register
         je       @True             // p=2, is prime
         jnc      @False2           // p<2, is not prime
         // we now know that p>2 i.e. positive
         // let's check if it is integer and odd
         fld      st(1)             // ST0=p  ST1=2  ST2=p  ST3=0
         frndint                    // truncate decimals
         fcom     st(2)             // compare with p
         fstsw    ax                // status word into AX
         sahf                       // copy ah into flag register
         jne      @False2           // p is no integer
         fdiv     st,st(1)          // divide p by 2
         frndint                    // truncate decimals
         fmul     st,st(1)          // multiply by 2
         fcomp    st(2)             // compare with p and pop from stack
         fstsw    ax                // status word into AX
         sahf                       // copy ah into flag register
         je       @False1           // p is even
         // we know now that p is an odd integer > 2
         // now let's find the maximum divisior for prime test
         fld      st(1)             // ST0=p  ST1=2  ST2=p
         fsqrt                      // ST0=Sqrt(p)  ST1=2  ST2=p
         frndint                    // truncate decimals
         fistp    qword ptr [ESP+4] // save int(sqrt(p)) and remove from stack
         fwait
         mov      ecx,[ESP+4]       // ECX=int(sqrt(p))
         sub      ecx,1             // subtract (smallest divisor - 2)
         shr      ecx,1             // ECX=number of divisors
         je       @True             // p is prime (3, 5 or 7)
         // now we can start the prime test
         fld1                       // ST0=1(divisor)  ST1=2  ST2=p  ST3=0
         // --------------------------------------------------------------
         // INCREASE DIVISOR BY 2
         // --------------------------------------------------------------
@Loop:   fadd     st,st(1)          // ST0=ST0+2
         // --------------------------------------------------------------
         // TEST IF P IS DIVISIBLE BY DIVISOR
         // --------------------------------------------------------------
         fld      st(2)       // ST0=p  ST1=divisor  ST2=2  ST3=p  ST4=0
         fprem                // ST0=p mod d  ST1=divisor  ST2=2  ST3=p  ST4=0
         fcomp    st(4)       // ST0=0 ? - ST0=divisor  ST1=2  ST2=p  ST3=0
         fstsw    ax          // control word into AX
         sahf                 // copy ah into flag register
         je       @False1     // p is divisible by d, i.e. p is not prime
         sub      ecx,1       // counter - 1
         jne      @Loop       // loop until counter = 0
@True:   mov      eax,1       // result=true
         jmp      @End
@False1: fstp     tbyte [edx]
@False2: xor      eax,eax     // result=false
         //---------------------------------------------------------------
         // finally clear fpu-stacks, reset control word and ESP
         // --------------------------------------------------------------
@End:    ffree    st(3)
         ffree    st(2)
         ffree    st(1)
         ffree    st(0)
         fldcw    [ESP].Word
         add      esp,12
end;

Gruß, Klaus

[negaH]

Hi

1.) du bewegst dich im Gebiet der Zahlentheorie er Ganzzahlen. Hier hat die Fließkommaarithmetik nichts zu suchen wenn du nicht sicherstellen kannst das sie auch 100% exakte Resultate liefert:

Zitat:

.....

Delphi-Quellcode:

 while ((u mod 2) = 0) do 
    begin 
      u := Round(u/2); 
      Inc(t); 
    end;

.....

Delphi-Quellcode:

  u := u div 2;
// oder
  U := u shr 1;

wären besser.

oder noch besser gleich optimiert

Delphi-Quellcode:

while odd(u) do 
begin
  u := u shr 1;
  inc(t);
end;

2.) der RM Test muß immer eine Primzahl als Primzahl identifzieren ! Tut er das nicht, dann ist irgendwas falsch an deinem Code. Die Aussage in denen der RM Test sich irrt bezieht sich nur darauf das er eine Zusammengesetzte Zahl als Primzahl ausgeben kann. Das sind dann meistens Carmichael Zahlen.

3.) hier im Forum solltest du ein besseren Ersatz für obige FIsPrimeNPA() finden, der schneller ist und zudem nicht mit der schlechst möglichen Komplexität arbeitet.

4.) du musst im RM Test nicht mit einem zufalls A rechnen. Es reicht zb. mit den kleinen Primzahlen beginnend mit 2 bis X als A zu arbeiten.

5.) ab diesem Moment kannst du dein Verfahren so abändern das es sequientiell die Zahlen durchgeht

[edit]

Alles schon mal da gewesen, suche !

http://www.delphipraxis.net/internal...ation&start=10
http://www.delphipraxis.net/internal...=exponentation

[/edit]

Gruß Hagen

[Amateurprofi]

Zitat:

Zitat von negaH

Hi

3.) hier im Forum solltest du ein besseren Ersatz für obige FIsPrimeNPA() finden, der schneller ist und zudem nicht mit der schlechst möglichen Komplexität arbeitet.

Hagen,

hab ich versucht - ohne Erfolg.
Wo genau ist denn der von Dir erwähnte bessere Ersatz zu finden ?

Gruß, Klaus

[St.Pauli]

@Amateurprofi, negaH und Der_Unwissende: Danke für eure Hilfe!!! :thumb:

Zitat:

Zitat von Amateurprofi

Naja St.Pauli, dann hab ich mal deine Hausaufgaben gemacht.

:wink: THX!

Zitat:

Zitat von Amateurprofi

In der Funktion ModularExponentiation heit es im Pseudo Code

Zitat:

Let (bk, bk-1,...b0) be the binary representation of b
for i &lt;-- k downto 0

Bei Dir wird das so umgesetzt daß, b als Binärstring in b_bin gestellt wird.
b_bin[Length(b_bin)] ist dabei Bit0.
Auf den ersten Blick sieht das so aus, wie im Pseudo Code, aber beim zweiten Blick eben nicht.
Die Schleife for i:=k downto 0 faßt zuerst hast höchstwertige Bit an und zuletzt Bit0.
Bei Deiner Schleife for i:=Length(b_bin) downto 1 ist es genau umgekehrt.

Ja stimmt. Danke - ich hatte mich so auf die Variabeln u und t konzentriert, um den Fehler dort zu finden, dass ich den ganzen Rest vollkommen vergessen hatte...

Zitat:

Zitat von Amateurprofi

Ich hab das korrigiert und dann für alle (ungeraden) Zahlen von 3 bis 65535 probiert. Bis 46337 bringt das Programm korrekte Ergebnisse. Für höhere Werte sagt es nur noch composite.... Hab nicht weiter geprüft, warum.

Ich vermute, dieser Fehler hat wieder mit der Binär-Darstellung von b zu tun. Werde mich gleich nochmal in Ruhe damit auseinander setzten :coder:

@ negaH: Danke für deine Optimierung. Zu zweitens: Davon haben wir ja gesprochen. Es muss sich um einen Fehler in der Implentierung handeln! Und zu den Suchvorschlägen: Da hatte ich wohl die falschen Suchbegriffe (Miller & Rabin) - und aus den Ergebnissen wurde mir nicht geholfen... :oops: