4. Solovay-Strassen-Algorithmus

4.1 Allgemeine Beschreibung

Ein weiteres Beispiel für einen probabilistischen Primzahltest ist der Algorithmus von Robert Solovay und Volker Strassen. In diesem Algorithmus wird das Jacobi-Symbol J_(a,p) benutzt, um zu testen, ob p prim ist. Das Jacobi-Symbol J_(a,n) ist ein mathematisches Symbol, das für alle ganzen Zahlen a und alle ungeraden ganzen Zahlen n definiert ist und beim Test bezüglich der Primzahleigenschaften von Zahlen eine Rolle spielt. Es kann über mehrere Formeln berechnet werden.

Algorithmus zur rekursiven Berechnung des Jacobi-Symbols

1. J_(1,n) = 1
2. J_{(a *b,n)} = J_(a,n) *J_(b,n)
3. J_(2,n) = 1, falls (n²-1)/8 gerade ist, sonst -1
4. J_(a,n) = J_{((a mod n),n)}
5. J_(a,b1*b2) = J_(a,b1) *J_(a,b2)
6. Ist ggT_(a,b) = 1 und a und b sind ungerade, so gilt:
   • J_(a,b) = J_(b,a), falls (a-1)(b-1)/4 gerade ist
   • J_(a,b) = -J_(b,a), falls (a-1)(b-1)/4 ungerade ist

Dieser Algorithmus ist in der C-Funktion long jacobi(long a,long b) implementiert.

Der Solovay-Strassen Primzahltest

Der Algorithmus testet folgendermaßen, ob eine Zahl p prim ist:

1. Wähle eine ungerade Zufallszahl als Kandidat für die Primzahl.
2. Wähle eine Zufallszahl a kleiner als p.
3. Falls ggT_(a,p) ≠ 1, ist p keine Primzahl.
4. Berechne j=a(p-1)/2 mod p
5. Berechne das Jacobi-Symbol J_(a,p).
6. Falls j ≠ J_(a,p), so ist p definitiv nicht prim, sonst liegt die Wahrscheinlichkeit, daß p nicht prim ist, bei höchstens 50 Prozent.

Struktogramm Solovay-Strassen-Primzahltest

4.2 Sourcecode

#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <time.h>


#define DURCHLAEUFE 10


/****************************************************************/
/* long random32(long p)                                        */
/* übergebener Parameter: long p                                */
/* Rückgabewert: long                                           */
/* Die Funktion random() kann nur mit 16-Bit-Integerzahlen      */
/* arbeiten.                                                    */
/* Die Funktion long random32(long) vergrößert diesen           */
/* Zahlenbereich auf 32 Bit und benutzt dafür die               */ 
/* random-Funktion.                                             */
/****************************************************************/
long random32(long p)
{
  int int1,int2;
  long a;

  int1=p>>16;
  int2=random(int1);
  a=int2;
  a<<=16;
  int1=p & 0x0000FFFF;
  int2=random(int1);
  a|=int2;
  return (a);
}

/****************************************************************/
/* long ggt(long x,long y)                                      */
/* übergebene Parameter: long x, long y                         */
/* Rückgabewert: long                                           */
/* Die Funktion ggt berechnet den größten gemeinsamen Teiler    */
/* zweier eingegebener 32-Bit-Integerzahlen mit dem Euklidschen */
/* Algorithmus und gibt ihn zurück.                             */
/****************************************************************/
long ggt(long x,long y)
{
  long g;
  if (x<0)
    x=-x;
  if (y<0)
    y=-y;
  if (x+y==0)
  {
    printf("\nFehler!\n");
    exit(1);
  }
  g=y;
  while (x>0)
  {
    g=x;
    x=y%x;
    y=g;
  }
  return (g);
}

/****************************************************************/
/* long exponent(long x, long y, long n)                        */
/* übergebene Parameter: long x, long y, long n                 */
/* Rückgabewert: long                                           */
/* Die Funktion berechnet x hoch y mod n und arbeitet mit dem   */
/* Prinzip der Verkettung oder binärem Quadrieren und           */
/* Multiplizieren.                                              */
/* Das Ergebnis wird als Rückgabewert zurückgegeben.            */
/****************************************************************/
long exponent(long x,long y, long n)
{
  long s,t,u;
  int i;

  s=1;
  t=x;
  u=y;

  /* Schleife ausführen, solange u ungleich 0 ist */
  while (u)
  {
    /* Wenn das unterste Bit der Zahl 1 ist */
    if (u&1)
      s=(s*t)%n;

    /* Die Bits von u um 1 nach rechts schieben */
    u>>=1;
    t=(t*t)%n;
  }
  return (s);
}

/****************************************************************/
/* long jacobi(long a, long b)                                  */
/* übergebene Parameter: long a, long b                         */
/* Rückgabewert: long                                           */
/* Die Funktion jacobi berechnet das Jacobi-Symbol für zwei     */
/* 32-Bit-Integerzahlen und gibt es als Rückgabewert zurück.    */
/****************************************************************/
long jacobi(long a,long b)
{
  long g;

  if ((b%2)!=1)
    printf("\nFehler! b gerade!\n");

  if (a>=b) a%=b;/* Regel 4 */
  if (a==0) return (0);        /* Definition 2 */
  if (a==1) return (1);  /* Regel 1 */

  if (a<0)
    if (((b-1)/2)%2 == 0)
      return (jacobi(-a,b));
  else
    return (-jacobi(-a,b));

  if (a%2==0)                    /* a ist gerade */
    if (((b*b-1)/8)%2==0)
      return (jacobi(a/2,b));
  else
    return (-jacobi(a/2,b));/* Nach den Regeln 3 und 2 */

  g=ggt(a,b);

  if (g==a)                    /* a geht in b auf */
    return (0);                    /* nach Regel 5 */
  else
    if (g!=1)
      return (jacobi(g,b)*jacobi(a/g,b));     /* nach Regel 2 */
    else
      if ((((a-1)*b-1)/4)%2==0)
        return (jacobi(b,a));                   /* nach Regel 6a */
    else
  return (-jacobi(b,a));                  /* nach Regel 6b */

}

/****************************************************************/
/* int solovay_strassen(long p,int t)                           */
/* übergebene Parameter: long p, int t                          */
/* Rückgabewert: int                                            */
/* Die Prozedur testet die Long-Int-Zahl p mit dem              */
/* probabilistischen Primzahltest von Robert Solovay und Volker */
/* Strassen auf Primzahleigenschaften. Der Algorithmus wird     */
/* iterativ mehrmals durchlaufen.                               */
/* Die übergebene Variable t gibt die Zahl der Durchläufe an.   */
/* Ist die Zahl keine Primzahl, wird 0 zurückgegeben, wenn die  */
/* Zahl die  Tests bestanden hat, also wahrscheinlich prim ist, */
/* wird 1 zurückgegeben.                                        */
/****************************************************************/
int solovay_strassen(long p,int t)
{
  int i,int1,int2;
  long a,j,jac;

/* Die Schleife wird t=i - mal wiederholt */
  for (i=t; i>=0; i--)
  {
    /* a ist eine Zufallszahl kleiner p */
    a=random32(p);

    if (ggt(a,p)!=1)
      return (0);
    j=exponent(a,(p-1)/2,p);

    /* Das Jacobi-Symbol von a und p wird berechnet */
    jac=jacobi(a,p);

    /* Wenn j ungleich jac ist p definitiv nicht prim */
    if (j!=jac)
      return (0);
  }
  /* p hat alle Testdurchläuufe bestanden 
     und ist wahrscheinlich prim */
  return (1);
}

/****************************************************************/
/* Hauptprozedur                                                */
/****************************************************************/
void main(void)
{
long zahl;
randomize();

/* Eine zufällige Zahl kleiner 2147483647 wird erzeugt. Weil
random() einen int-Wert zurückgibt, wird die 32-Bit-Zahl zahl
mit der Funktion random32() aus dem Bitmuster zweier int-Zahlen
erzeugt */

zahl=random32(0x7FFFFFFF);

/* Das höchstwertigste Bit der Zahl wird auf Null gesetzt,
damit die signed long - Zahl nicht negativ ist.
Dann wird das niederwertigste und das zweit-höchstwertige Bit auf 1
gesetzt, damit sichergestellt ist, daß die Zahl groß und 
ungerade ist */

zahl &= 0x7FFFFFFF;
zahl |= 0x40000001;

/* Der Solovay-Strassen-Test wird auf die Zahl angewendet.
Der Test wird DURCHLAEUFE-mal durchlaufen. */

if (solovay_strassen(zahl,DURCHLAEUFE))
  printf("\n%ld ist vermutlich eine Primzahl. \n",zahl);
else
  printf("\n%ld ist keine Primzahl. \n",zahl);

getch();
}