Die Computerseite von Eckart Winkler
Primzahltests und Primfaktorzerlegung - Teil 2

 

Dieser Artikel ist im Heise-Verlag in der Zeitschrift "c't", Ausgabe 05/1989 erschienen. Er ist hier aus Gründen der besseren Übersicht dreigeteilt. Grundlage ist der Artikel "Lange Integer-Zahlen in C". Links:


Monte-Carlo-Tests

Nachdem unsere Ergebnisse, die für Micro-Rechner doch beachtlich sind, (hoffentlich) weiteres Interesse geweckt haben, wollen wir uns nun mit der Frage beschäftigen, wie wir für beliebige Zahlen die Eigenschaft "prim" oder "nicht prim" entscheiden können.

Hier muß beim momentanen Stand der Dinge jedoch ein Dämpfer erfolgen. Zwar sind die Bemühungen weiterhin im Gange, konkret verwertbare Ergebnisse sind jedoch selten. So hat J.Pollard einen Test konstruiert, der sich auf eine gewisse Zahl c stützt, deren Wert leider niemand kennt. Ein anderer Test von Gary Miller kann nur als korrekt angesehen werden, wenn die sog. "Riemannsche Vermutung" bewiesen wird, was aber bisher nicht geschehen ist. Beide sind daher in der Praxis nicht anwendbar.

So kam man auf die Idee der "probabilistischen Primzahltests", auch "Monte-Carlo-Tests" genannt. Was haben wir darunter zu verstehen? Nun, ein Primzahltest, wie wir ihn vor Augen haben, liefert immer eine Aussage der Form: "n ist prim" oder "n ist nicht prim". Ein probabilistischer Test hingegen sagt aus: "n ist mit einer Wahrscheinlichkeit w eine Primzahl". D.h. wir können nicht ganz sicher davon ausgehen, daß n prim ist. Die Wahrscheinlichkeit w wird aber sehr klein sein, so daß auch das Risiko, daß wir eine falsche Antwort erhalten, sehr gering ist.

Eine Anmerkung noch zur Namensgebung: "Probabilistisch" kommt von dem lateinischen Wort "probabilis", was "wahrscheinlich" heißt. Die Bezeichnung "Monte-Carlo-Test" spielt auf die praktische Anwendung derartiger Tests an. Hierfür wird nämlich die Bildung einer Zufallsfolge erforderlich. Dabei könnte neben der Definition von Pseudo-Zufallszahlen auch ein Roulette behilflich sein. Und solche findet man unter anderem in der Spielbank von Monte Carlo.

Der Test, den wir nun besprechen, basiert auf dem schon erwähnten Test von Gary Miller. Die Gültigkeit von diesem ist ja noch nicht erwiesen, und so führte M.O.Rabin einige Änderungen durch. Der im folgenden angegebene Satz ist daher als Satz von Miller und Rabin bekannt: 

Satz von Miller und Rabin:



Es gelte n-1 = 2^k*u, u ungerade. Wenn für m zufällig gewählte 
Werte a<n die Bedingungen
a^u mod n = 1  oder  a^(2^j*u) mod n = n-1  für ein j mit 0<=j<k 
erfüllt sind, kann man mit einer Fehlerwahrscheinlichkeit von bis 
zu 4^(-m) die Zahl n als Primzahl ansehen.

Was haben wir also bei der Durchführung des Miller-Rabin-Tests zu tun? Zunächst spalten wir von n-1 alle Zweier ab. Der ungerade Anteil von n-1 heißt u. Die Anzahl der in n-1 enthaltenen Zweier wird in k festgehalten. Nun wählen wir uns eine gewisse Anzahl von Zufallszahlen a. Für jede führen wir den angegebenen Test aus.

Zwei Bedingungen sind gegeben, nur eine muß erfüllt sein. Gilt z.B. a^u mod n = 1, hat das momentan aktuelle a den Test erfüllt, und wir können zum nächsten Zufallswert übergehen. Ist die Bedingung nicht erfüllt, ist noch nichts verloren. Dann prüfen wir die Gültigkeit der zweiten Gleichung. Hierfür haben wir sogar noch mehrere Möglichkeiten, da das hier vorkommende j mehrere Werte annehmen kann. Dies hängt von der Größe von k ab.

Dann würden wir also prüfen, ob a^u mod n = n-1 ist. Ist dies nicht der Fall, würde uns auch die Gültigkeit von a^(2*u) mod n = n-1 genügen. Oder a^(2^2*u) mod n = n-1, usw. Wenn auch nur eine einzige dieser Gleichungen richtig ist, ist dieser Teil des Tests mit positivem Ergebnis abgeschlossen. Sind jedoch alle falsch, wissen wir, daß n nicht prim sein kann.

Der Test beruht auf der Tatsache, daß mindestens eine der Bedingungen für jede Primzahl gilt, für eine zusammengesetzte Zahl jedoch nur mit der Wahrscheinlichkeit 1/4. Führen wir den Test also einmal mit Erfolg durch, wissen wir: Diese Zahl ist mit Wahrscheinlichkeit 3/4 eine Primzahl. Nach dem nächsten erfolgreichen Durchlauf beträgt die Fehlerwahrscheinlichkeit nur noch 1/16, nach dem dritten 1/64 usw. Eine vernünftige Anzahl von Durchläufen ist 30. Hat eine Zahl alle 30 Tests bestanden, ist sie mit einer Fehlerwahrscheinlichkeit von 4^(-30)=8.67E-19 eine Primzahl. Rabin selbst ist der Meinung, dies sei klein im Vergleich zur Häufigkeit der praktisch auftretenden Computerfehler.

Wie schon beim Brillhart-Selfridge-Test können wir wieder einiges an Zeit einsparen, indem wir alte Ergebnisse ausnutzen und nicht alles von Grund auf neu berechnen. So werden wir zunächst der Variablen w den Ausdruck a^u mod n zuweisen. Nachdem getestet wurde, ob dies gleich 1 ist, wird der Wert in der zweiten Bedingung schon wieder benötigt, nämlich für den Fall j=0. Es wird also sofort gefragt, ob w gleich n-1 ist.

Erst für die nächste Abfrage brauchen wir einen neuen Wert. Der Unterschied zu unserer bisherigen Zahl besteht aber nur in einer 2 im Exponenten. Wir brauchen also lediglich w^2 mod n zu berechnen. Dies wird in einer Schleife k-mal durchgeführt. Tritt irgendwann der Wert n-1 auf, ist der Test erfüllt. Geht die Schleife zu Ende, ohne daß dies passiert wäre, ist die Zahl unwiderruflich nicht prim.

Dies ist die Realisierung des Miller-Rabin-Tests: 

#define MaxNr 50
#include "longcalc.c"

void main()
{
  longint z;
  int i;
  for (i=1;i<=10;i++) {
    scanl(z); mill_rab(z,30); }
}

/* User-Version von mr. */
void mill_rab(longint n, int p)
{
  if (mr(n,p)) printf(" - prim\n");
    else printf(" - nicht prim\n");
}

/* Unterzieht die Zahl n dem Primzahltest von Miller-Rabin. Der Test
   wird p-mal durchgefuehrt. Ausgabe: 1, wenn n prim; 0, wenn nicht. */
void mr(longint n, int p)
{
  int i,k;
  longint u,a;
  transl(n,u); decl(u); k=0;
  while (!oddl(u)) { div2l(u); k++; }
  longl(2L,a);
  for (i=1;i<=p;i++) {
    if (!mr_check(k,u,n,a)) return 0;
    abmodn(a,a,n,a); decl(a); }
  return 1;
}

/* Es sei n=2^k*u+1, u ungerade. Die Funktion prueft, ob die Zahl a
   das Miller-Rabin-Kriterium erfuellt. Ausgabe: 1: ja; 0: nein. */
void mr_check(int k, longint u, longint n, longint a)
{
  int i;
  longint w,n1;
  ahsmodn(a,u,n,w);
  if (eql(w,einsl)) return 1;
  transl(n,n1); decl(n1);
  for (i=0;i<k;i++) {
    if (eql(w,n1)) return 1;
    abmodn(w,w,n,w); }
  return 0;
}

Was wir uns eben überlegt haben, wird in der Routine mr_check durchgeführt. Die Variablen in der Argumentliste tragen dieselben Bezeichnungen wie im Satz. n1 wird gleich n-1 gesetzt. Dieser Wert soll ja erreicht werden, wenn n prim ist. Ist eine der Bedingungen erfüllt, wird 1 zurückgegeben, sonst 0.

Die Prozedur mr führt den Test p-mal aus. Zu Beginn wird von der Zahl n der ungerade Anteil abgespalten und die Zahl k ermittelt. Dann geht es schon in die Schleife, die p-mal ausgeführt wird. Die Zufallsfolge ist hier sehr einfach gewählt. Wir setzen als Anfangswert a=2, in jedem Schleifendurchlauf ändern wir dies in a^2-1 mod n. Dies ist nur als Vorschlag zu verstehen, hier kann durchaus eine andere Bildungsvorschrift eingebaut werden, die Folge muß auch gar nicht rekursiv definiert sein. Liefert irgendwann mr_check den Wert 0, wird auch dieses Unterprogramm mit derselben Ausgabe abgebrochen. Erst wenn alle p Tests mit Erfolg beendet sind, kann eine 1 zurückgegeben werden.

mill_rab ist zur bequemen Anwendung des Tests gedacht, da die Resultate nicht über Codezahlen zurückgegeben werden, sondern im Klartext auf dem Bildschirm erscheinen. Das Hauptprogramm fragt in einer Schleife jeweils eine longint-Zahl ab. Diese wird dem Miller-Rabin-Test unterworfen, das Ergebnis wird durch Benutzung von mill_rab direkt auf den Bildschirm geschrieben.

Mit diesem Monte-Carlo-Test sind wir also in der Lage, beliebige Zahlen als Primzahlen zu erkennen, wenn auch mit einem gewissen Restrisiko. Beispielsweise können wir so feststellen, daß die am Anfang des Artikels genannte Zahl prim ist. Bei den durch den Brillhart-Selfridge-Test überprüften Zahlen kommen wir hier zum selben Ergebnis. Allerdings ist es nicht zu empfehlen, den Miller- Rabin-Test zur Berechnung von Rekordprimzahlen zu verwenden, das würde viel zu lange dauern.

Erwähnt werden soll auch die Tatsache, daß besonders Zahlentheoretiker derartigen Verfahren mit Skepsis begegnen. Ihnen nutzt eine Zahl, die mit einer Fehlerwahrscheinlichkeit von 8.67E- 19 prim ist, wenig. Denn dann müßten Sätze, die sich aus der Prim- Eigenschaft einer Zahl ergeben, mit derselben Fehlerwahrscheinlichkeit versehen werden, und das wäre mehr als ungewöhnlich. Für praktische Belange, wie sie vom RSA-Verfahren vorgegeben werden, mag dies indes genügen.


Primfaktorzerlegung

Wir wissen ja bereits, daß die Sicherheit dieses Verschlüsselungsverfahrens darauf beruht, daß eine Zahl n nicht in ihre Primfaktoren zerlegt werden kann. Trotzdem gibt es hierfür Verfahren, und eines davon wollen wir uns nun ansehen. Es wurde von J.Pollard entdeckt und basiert ebenfalls auf einer Zufallsfolge. Wenn auch der Zeitpunkt, zu dem sich ein Resultat einstellt, rein zufällig ist, ist dieses Verfahren doch von einer anderen Art als der Miller-Rabin-Test. Denn die Korrektheit des Ergebnisses kann sofort leicht überprüft werden. Das Verfahren lautet: 

Verfahren von Pollard:



Definiere eine rekursive Zufallsfolge, z.B. durch x(0)=2, 
x(i+1)=x(i)^2-1 mod n. Berechne im i-ten Schritt das Produkt aller 
Faktoren (x(2*j)-x(j)) von j=1 bis i. Dann berechne den ggT dieses 
Produkts und n. Ist er ungleich 1, hat er den Wert eines Faktors 
von n.

Was ist hier zu tun? Wir werden eine Variable q benutzen, in der das vorkommende Produkt gespeichert wird. In jedem Schritt wird q auf den neuesten Stand gebracht. Irgendwann steckt ein Teiler von n auch in q. Um ihn zu ermitteln, müssen wir nur den größten gemeinsamen Teiler von n und q berechnen.

Gedanken müssen wir uns über die Berechnung von q machen. Im i-ten Schritt ist q=(x(2)-x(1))*(x(4)-x(2))*...*(x(2*i)-x(i)) mod n. Das unterscheidet sich von dem Produkt im Schritt vorher nur durch den letzten Faktor (x(2*i)-x(i)). Dies muß natürlich mod n gerechnet werden. Wir müssen also nur den bisherigen Wert von q nehmen und erhalten durch q=q*(x(2*i)-x(i)) mod n den neuen Wert. Wir initialisieren q selbstverständlich mit 1.

Ein Problem taucht auf. Wir müssen in jedem Schritt zwei Werte der Zufallsfolge in die Rechnung einbeziehen, das ist der zuletzt berechnete, also x(2*i), und einer, der schon lange vorher ermittelt wurde, x(i). Wie bekommen wir dies in den Griff? Nun, wir könnten alle Werte der Folge in einem Vektor abspeichern. Dann könnten wir sie sogar mit den im Verfahren angegebenen Indizes ansprechen. Allerdings würden wir unseren Algorithmus damit auf eine Höchstzahl an durchzuführenden Schritten beschränken. Denn longint-Zahlen verbrauchen viel Speicherplatz, oft ist jedoch die Anwendung mehrerer 1000 Schritte notwendig. Beides läßt sich nicht vereinbaren.

Wir führen daher lieber zwei Variablen x und y mit. In beiden wird jeweils ein Element der Zufallsfolge gespeichert, in x das mit Index 2*i, in y das mit Index i. x muß in jedem Schritt also zweimal der Berechnungsvorschrift der Folge unterworfen werden, y einmal. Auf diese Weise werden zwar die Werte von y alle doppelt berechnet, aber der Speicher kann nie knapp werden, wir können beliebig viele Schritte durchführen.

Hier ist die Realisierung: 

#define MaxNr 50
#include "longcalc.c"

void main()
{
  longint z;
  int i;
  for (i=1;i<=10;i++) {
    scanl(z); pollard(z,10000,25); }
}

/* User-Version von poll. */
void pollard(longint n, int st, int ch)
{
  longint a,b;
  int f;
  poll(n,st,ch,a,&f);
  if (f) { 
    divl(n,a,b); printf("Faktoren: "); printl(a);
    printf("  "); printl(b); }
  else printf("Vermutlich prim");
  printf("\n");
}

/* Berechnet einen Faktor a von n nach der Methode von Pollard. In f
   wird 1 zurueckgegeben, falls die Suche erfolgreich war, sonst 0.
   s ist die Maximalzahl der Schritte, alle c Schritte wird der ggT
   von q und n berechnet. */
void poll(longint n, int st, int ch, longint a, int *f)
{
  longint x,y,q,r;
  int i;
  longl(2L,x); longl(2L,y); longl(1L,q);
  for (i=1;i<=st;i++) {
    abmodn(x,x,n,x); decl(x); abmodn(x,x,n,x); decl(x);
    abmodn(y,y,n,y); decl(y);
    if (gtl(x,y)) subl(x,y,r);
      else { subl(n,y,r); addl(r,x,r); }
    abmodn(q,r,n,r);
    if (eq0l(r)) {
      ggtl(n,q,a); if (nel(a,einsl)) *f=1; else *f=0; return; }
    else transl(r,q);
    if (i%ch==0) {
      printf("Schritt %d\n",i); ggtl(n,q,a);
      if (nel(a,einsl)) { *f=1; return; } } }
  *f=0; return;
}

Das wesentliche Unterprogramm ist poll. In der Parameterliste stehen fünf Variablen. n ist die zu zerlegende Zahl. In a wird ein Faktor zurückgeliefert, falls einer existiert. In f finden wir eine Codezahl, die über Erfolg oder Mißerfolg des Verfahrens Aufschluß gibt. Bei Erfolg wird sie auf 1 gesetzt, ansonsten auf 0.

Die Variablen st und ch beeinflussen die Laufzeit des Verfahrens. In st können wir die Zahl der durchzuführenden Schritte begrenzen, in ch wird angegeben, wie oft der ggT-Test gemacht werden soll. Denn wir müssen in jedem Schritt zwar x,y und q aktualisieren, eine Prüfung, ob schon ein Faktor gefunden ist, ist aber nicht in jedem Schritt zu empfehlen, denn auch sie kostet Zeit. Und ist der gesuchte Faktor einmal in q enthalten, geht er aufgrund der Definition des Verfahrens auch nicht mehr verloren. Schreiben wir also als drittes Argument eine 25, wird diese Überprüfung nur alle 25 Schritte durchgeführt.

Der Ablauf ist nun klar. x und y werden mit 2 initialisiert, q mit 1. Dann erfolgt die Schleife, die st-mal ausgeführt wird. Zu Beginn werden x und y auf den neuesten Stand gebracht. Nun müssen wir als neuen Faktor für q den Wert (x-y) nehmen. Sollte y größer als x sein, ziehen wir zunächst y von n ab und addieren hierzu x. Da wir mod n rechnen, ist dies korrekt. Unsere Prozedur subl funktioniert jedoch nur, wenn y nicht größer als x ist, deshalb diese Fallunterscheidung.

Der Faktor wird also in r gespeichert. Durch abmodn(q,r,n,r) aktualisieren wir q. Das Ergebnis wird jedoch zunächst in r abgelegt, da ein weiterer Sonderfall eintreten kann. Denn da unsere Zufallsfolge rekursiv definiert ist, müssen irgendwann x und y gleich sein. Dann ist der Faktor (x-y) gleich 0, und alle weiteren Werte von q wären 0. In diesem Fall muß das Verfahren abgebrochen werden.

Es könnte aber sein, daß mittlerweile der gesuchte Faktor von n gefunden ist. Deshalb führen wir den ggT-Test ein letztes Mal durch. Der vorherige Wert von q ist ja noch vorhanden. Das Resultat wird ausgegeben. Bei kleinen Zahlen n kann es hier leicht passieren, daß Faktoren nicht gefunden werden. Deshalb ist eine Anwendung nur für große Zahlen zu empfehlen.

Weiter im else-Teil der Abfrage: Ist also r ungleich 0, können wir wie gewohnt weitermachen. Das heißt in diesem Fall, wir übertragen den Wert von r nach q. Nun kann eine reguläre Überprüfung des ggT erfolgen, dies aber nur alle ch Schritte. Als Bedingung, wann es wieder Zeit hierfür ist, dient der Ausdruck i%ch==0, % berechnet den Rest der Division i/ch. Ist der ggT gleich 1, geht es weiter im Verfahren, ansonsten ist ein Faktor gefunden. Als Service für den Benutzer wird an dieser Stelle ausgegeben, wieviele Schritte bereits durchgeführt wurden.

Zum bequemen Gebrauch ist die Prozedur pollard angegeben, die die Resultate direkt auf den Bildschirm schreibt. Wurde ein Faktor gefunden, kann durch Division auch ein zweiter Faktor ausgerechnet werden. Das Hauptprogramm ist fast identisch zu dem in Listing 2. In einer Schleife werden durch scanl longint-Werte abgefragt. Von diesen werden jeweils zwei Faktoren ermittelt und auf den Bildschirm geschrieben.

Bei der praktischen Erprobung zeigte es sich, daß das Verfahren umso länger läuft, je größer die Primfaktoren von n sind. Die eigentliche Größe von n hat hiermit nur wenig zu tun. So wurde bei der Mersenneschen Zahl 2^87-1 schon bei der ersten Überprüfung nach 25 Schritten der Faktor 1798993 gefunden. Später stellte sich heraus, daß hierin sogar die 7 enthalten ist, dies war der Grund für die Schnelligkeit des Verfahrens.

Wesentlich länger dauerte die Faktorisierung der bereits erwähnten Zahl 2^67-1. Sie zerfällt in die Primfaktoren 193707721 und 761838257287. Obwohl die Zahl selbst nur über 21 Stellen verfügt, waren 18875 Schritte des Verfahrens notwendig, um diese zu entdecken. Die Faktorisierung einer 100-stelligen Zahl mit zwei etwa gleichgroßen Primfaktoren ist daher auch mit dieser Methode undenkbar.

 

Übersicht Programmiersprache C Index