Pseudoprimzahlen: Programme

Modulare Arithmetik

Da man es bei der Berechnung von fermatschen Pseudoprimzahlen häufig mit großen Potenzen großer Zahlen zu tun bekommt, hat man ein kleines Problem. Zum Beispiel für den Nachweis, dass 341 eine fermatsche Pseudoprimzahl zur Basis 2 ist, muss man die Zahl $2^{340}$ berechnen. Diese Zahl heißt ausgeschrieben:

2239744742177804210557442280568444278121645497234649534899989100963791871180160945380877493271607115776

Das ist für ein Programm wie MuPad wahrscheinlich kein Problem, für die meisten herkömmlichen Programmiersprachen aber sehr wohl. Durch eine Kombination von Multiplikation und Modulo kann man aber auch wesentlich größere Potenzen kleinhalten.
Da bei Rechenoperationen Modulo $n$ nur Zahlen bis $n-1$ auftreten, ist das Ergebnis einer Multiplikation vor der Modulooperation immer kleiner als $n^2$.

Beispiel: $5^6=15625$. Letztendlich interessiert uns ja aber nicht $5^6=15625$, sondern vielmehr $5^6\mathrm{mod}n$ mit sagen wir mal $n=7$. Das ist $5^6\mathrm{mod}7=1$. Nun lässt sich $5^6\mathrm{mod}7$ auch so berechnen:

$(((((((((5\ast 5)\mathrm{mod}7)\ast 5)\mathrm{mod}7)\ast 5)\mathrm{mod}7)\ast 5)\mathrm{mod}7)\ast 5)\mathrm{mod}7$

Nach der Modulo-Operation ist der Wert immer kleiner als 7.

// Programm in Pseudocode zur Veranschaulichung: 
funktion powermod(Basis, Exponent, N)
    Produkt = 1
    k = 0
    while k < Exponent 
        Produkt = (Produkt * Basis) mod N
        k = k + 1
    return Produkt

Binäres Potenzieren

Für Primzahltests nach dem kleinen fermatschen Satz müssen große Potenzen berechnet werden; dabei wäre es extrem zeitaufwendig, x^k durch k - 1 Multiplikationen zu berechnen. Ein effiziente Berechnung ist durch binäre Exponentiation möglich:

Gegeben seien die Basis x und der Exponent k.

Berechnungsschritte:

setze P = 1 (Produkt) und b = x
wiederhole folgende Schritte solange k > 0 ist:
	wenn k ungerade ist, multipliziere P mit b → P und vermindere k um 1 → k	P = P * b k = k - 1
	quadriere b → b	b = b * b
	halbiere k → k	k = k / 2
Ergebnis: P = xk, wenn k = 0 ist

Bei Berechnung Modulo $n$ wird das Ergebnis jeder Multiplikation mod $n$ reduziert, wie im vorherigen Absatz beschrieben.

In Pseudocode:

ohne Modulo	Modulo N
function Potenz(x, k) { Produkt = 1 Faktor = x Exponent = k while Exponent > 0 { if Exponent mod 2 == 1 // wenn der Exponent ungerade ist: { Produkt = Produkt * Faktor Exponent = Exponent - 1 } Exponent = Exponent / 2 Faktor = Faktor * Faktor } return Produkt }	function Potenz(x, k, N) { Produkt = 1 Faktor = x Exponent = k while Exponent > 0 { if Exponent mod 2 == 1 // wenn der Exponent ungerade ist: { Produkt = (Produkt * Faktor) mod N Exponent = Exponent - 1 } Exponent = Exponent / 2 Faktor = (Faktor * Faktor) mod N } return Produkt }

Statt mod 2 kann in Zeile 9 auch das (schnellere) bitweise Und mit 1 verwendet werden, also Exponent = Exponent & 1. Die Division durch 2 in Zeile 15 kann durch Schieben um ein Bit nach rechts ersetzt werden: Exponent = Exponent >> 1.

Insgesamt sind mit dieser Methode nur maximal 2 Multiplikationen pro Bit des Exponenten notwendig, also $2\cdot \lceil lo{g}_{2}k\rceil$.

Liste zusammengesetzter Zahlen

Um Primzahltests zu prüfen bzw. Pseudoprimzahlen zu finden, wäre eine Funktion, die eine Liste (ungerader) zusammengesetzter Zahlen erstellt, nützlich. Sie spart zum einen die Überprüfung aller Zahlen, die den Test bestehen, mit einem im untersuchten Zahlenbereich sicheren Primzahltest, um Pseudoprimzahlen von Primzahlen zu unterscheiden, und zum anderen werden Primzahlen von vornherein nicht getestet.

Mit einem Algorithmus, der wie das Sieb des Eratosthenes funktioniert, lässt sich ein solches Hilfsmittel realisieren:

• Erstelle eine leere Menge für das Ergebnis
• Für alle ungeraden n von 3 bis Obergrenze:
  • wenn n nicht in der Menge enthalten ist:
    • nimm alle ungeraden Vielfachen von $n$ ab $n^2$ bis zur Obergrenze, also $(n^2+2k\cdot n)$, in die Menge auf
• Sortiere die Elemente der Menge aufsteigend

Implementierung in Python

def odd_composites(lim):
    result = set()
    for n in range(3, lim, 2):
        if n not in result:
            for m in range(n*n, lim, 2*n):
                result.add(m)
    return sorted(result)

Die Funktion lässt sich leicht so erweitern, dass auch gerade Zahlen zurückgegeben werden; diese sind aber für Primzahltests weniger interessant.

Carmichael-Zahlen

Ein Programm, das eine Liste von Carmichael-Zahlen ausgibt, ist geradezu trivial. Dass dem so ist, verdanken wir dem Kriterium von Korselt. Dieses Theorem besagt, dass eine Carmichel-Zahl $c$eine quadratfreie, aus mindestens drei Primfaktoren bestehende, natürliche Zahl der Form $c=p_1\cdot p_2\cdot \mathrm{...}\cdot p_n$ ist, so dass für jedes $p_i$ gilt, dass $p_i-1$ die Zahl $c-1$ teilt. Daraus folgt, dass man sich um den kleinen Fermatschen Satz $a^{n-1}\equiv 1(\mathrm{mod}n)$ nicht im geringsten zu kümmern braucht. Neben Korselts Kriterium gibt es auch noch Formeln zum Erzeugen spezieller Carmichael-Zahlen. Auch bei diesen Beispielen braucht man sich um den kleinen Fermatschen Satz nicht zu kümmern.

zweiter Ansatz

Das folgende Programm benutzt drei ineinander verschachtelte Schleifen, durch die sichergestellt wird, dass drei zueinander teilerfremde Primzahlen ausgewählt werden. Aus diesen drei Primzahlen wird ein Produkt gebildet, das die zu testende Zahl darstellt. An dieser Zahl wird nun Korselts Kriterium getestet, was einfach ist, da dem Programm die Primteiler der zu testenden Zahl bekannt sind. Ist Korselt Kriterium erfüllt, dann handelt es sich bei der Zahl um eine Carmichael-Zahl, und die Zahl wird, samt Primzahlfaktoren als Carmichael-Zahl ausgegeben.

Details: Die Variable p.0 enthält die Anzahl der unterschiedlichen Primzahlen. In p.1 bis p.(p.0) sind die einzelnen Primzahlen abgelegt.

p.1  = 3
p.2  = 5
p.3  = 7
p.4  = 11
p.4  = 13
p.5  = 17
p.6  = 19
p.7  = 23
p.8  = 29
p.9  = 31
p.10 = 37
.
.
.
p.549 = 4001
p.550 = 4003
p.551 = 4007
p.552 = 4013
p.553 = 4019
p.554 = 4021
p.555 = 4027
p.556 = 4049
p.557 = 4051
p.558 = 4057
i=558
do a=1 to (i-2)
  do b=a+1 to (i-1)
    do c=b+1 to i
      t = 0
      z=p.a*p.b*p.c
      ax=p.a-1
      bx=p.b-1
      cx=p.c-1
      zax=(z/p.a)-1
      zbx=(z/p.b)-1
      zcx=(z/p.c)-1
      pz=(zax // ax) + (zbx // bx) + (zcx // cx)
      if (pz = 0) then t = 1
      if (t = 1) then do
        say z||"="||p.a||"*"||p.b||"*"||p.c
        lineout("rcrmn___.txt",z||" = "||p.a||"*"||p.b||"*"||p.c)
      end
      t=0
    end
  end
end

Das Programm kann man modifizieren, indem man die Anzahl der Primzahlen erhöht. Für Carmichael-Zahlen mit mehr als drei Primfaktoren muss man die Anzahl der Schleifen erhöhen, und die Prüfroutinen erweitern.

Carmichael-Zahlen nach Chernick

Die Carmichel-Zahlen nach Chernick sind nur ein kleiner Ausschnitt aus den Carmichael-Zahlen. Sie lassen sich durch folgende Formel generieren: $M_3(n)=(6n+1)\cdot (12n+1)\cdot (18n+1)$. Es gibt dabei allerdings eine Einschränkung:

• Jeder der drei Faktoren $(6n+1)$, $(12n+1)$ und $(18n+1)$ muss eine Primzahl sein. Halt, hier gibt es eine Ausnahme: Wenn ein einziger der drei Faktoren selbst eine Carmichael-Zahl ist, so ist $M_3(n)$ trotzdem eine Carmichael-Zahl.

Da alle Carmichael-Zahlen quadratfrei sind, ist dann auch jeder der drei Faktoren $(6n+1)$, $(12n+1)$ und $(18n+1)$ quadratfrei - auch wenn einer davon eine Carmichael-Zahl ist. Da die drei Faktoren paarweise teilerfremd sind, ist auch der Ausdruck $M_3(n)=(6n+1)\cdot (12n+1)\cdot (18n+1)$quadratfrei.

Alle gültigen $n$ für Carmichael-Zahlen nach Chernick enden im Dezimalsystem mit einer 0, einer 1 einer 5 oder einer 6.

l.1  = 0
l.2  = 1
l.3  = 5
l.4  = 6
do i = 0 to 1000
  do j = 1 to 4
    n=i+l.j
    m=(6*n+1)*(12*n+1)*(18*n+1)
    if((primzahl(6*n+1)=true) && (primzahl(12*n+1)=true) && (primzahl(18*n+1)=true)) then do
      say m||"="||(6*n+1)||"*"||(12*n+1)||"*"||(18*n+1)
      lineout("rcrmn___.txt",m||" = "||(6*n+1)||"*"||(12*n+1)||"*"||(18*n+1))
    end
  end
end

Das Programm lässt sich entsprechend der Formeln für Carmichael-Zahlen nach Chernick, mit mehr als drei Primzahlfaktoren, nach der entsprechenden Formel erweitern.

Absolute Eulersche Pseudoprimzahlen

Das Programm, mit dem oben Carmichael-Zahlen erzeugt werden, kann man mit einer kleinen Modifikation dazu bringen, absolute eulersche Pseudoprimzahlen zu erzeugen:

      zax=(z-1)/2
      pz=(zax // ax) + (zax // bx) + (zax // cx)

statt

      zax=(z/p.a)-1
      zbx=(z/p.b)-1
      zcx=(z/p.c)-1
      pz=(zax // ax) + (zbx // bx) + (zcx // cx)

Das hierbei alle absoluten eulerschen Pseudoprimzahlen ausgegeben werden, ist eine andere Frage. Aber die Zahlen, die ausgegeben werden, sind absolute Eulersche Pseudoprimzahlen.

Fermatsche Pseudoprimzahlen

Bei Primzahlen ist im Gegensatz zu fermatschen Pseudoprimzahlen und Carmichael-Zahlen die Kongruenz $a^{n-1}\equiv 1(\mathrm{mod}n)$ für alle Basen $a$ mit $2\le a<n$ erfüllt. Carmichael-Zahlen unterscheiden sich von anderen fermatschen Pseudoprimzahlen dadurch, dass der Anteil der Primbasen $a$ mit $2\le a<n$, für die $a^{n-1}\equiv 1(\mathrm{mod}n)$ ist, größer als 97% ist. Der Einfachheit halber benutzt man als Basen $a$ nur die Primzahlen, die kleiner als die zu testende Zahl $n$ sind. Das Programm testet also zu jeder Primbasis $a$, ob eine Zahl $n$ die Bedingung $a^{n-1}\equiv 1(\mathrm{mod}n)$ erfüllt. Erfüllt sie sie, wird die Variable tm1 = 1 gesetzt und gleichzeitig die Zählvariable tm3 um eins erhöht. Erfüllt die Zahl $n$ die Bedingung nicht, wird die Variable tm2 = 2 gesetzt. Nach Ablauf der Tests wird die Variable gtm = tm1 + tm2 gesetzt. Aus gtm und dtm lässt sich ablesen, ob es sich bei der Zahl $n$ um eine Primzahl, eine fermatsche Pseudoprimzahl oder eine Carmichael-Zahl handelt.

gtm	dtm	Art der Zahl
1	-	Primzahl
3	tm3 < dtm	fermatsche Pseudoprimzahl
3	tm3 >= dtm	Carmichael-Zahl

/* Ein Programm zur Ermittlung von Primzahlen, Pseudoprimzahlen */
/* und Carmichaelzahlen   */
/* Ein langsames Programm */

#include <stdio.h>

int primfeld[400000];
int tst[400000];

unsigned long modup(unsigned long base, unsigned long y)
{
  unsigned long exponent, product = 1;
  exponent = y-1;
  while (exponent)
  {
    /* wenn Exponent gerade ist, Basis quadrieren und Exponent halbieren bis der Exponent ungerade ist: */
    while (exponent & 1 == 0)
    {
      exponent >>= 1;
      base = (base * base) % y;
    }
    /* Exponent ungerade (immer nach obiger Schleife): Produkt *= Basis ... %= y, Exponent dekrementieren: */
    product = (product * base) % y;
    exponent--;
  }
  return(product);
} 

int main()
{
  unsigned long index, index2, anzp=1, m, dtm;
  int tm1, tm2, tm3, gtm;
  FILE *prim;
  FILE *pspr;
  FILE *cnbr;
  prim = fopen("prim.dat","w");
  pspr = fopen("pspr.dat","w");
  cnbr = fopen("cnbr.dat","w");
  primfeld[1] = 2;
  for(index=3;index<=4000000;index++)
  {
    tm1 = 0;
    tm2 = 0;
    tm3 = 0;
    /* faktor$ = "" */
    for(index2=1;index2<=anzp;index2++)
    {
      m = modup(primfeld[index2], index);
      tst[index2] = m;
      if (m == 1)
      {
        tm1 = 1;
        tm3++;
      }
      if ( m != 1) { tm2 = 2; }
    }
    gtm = tm1 + tm2;
    if (gtm == 1)
    {
      anzp++;
      primfeld[anzp] = index;
      fprintf(prim,"%d \n",index);
    }
    if (gtm == 3)
    {
      dtm=anzp/2;
      if (tm3 < dtm)
      {
        fprintf(pspr,"%d: ",index);
        for(index2=1;index2<=anzp;index2++)
        {
          m = modup(primfeld[index2], index);
          if (m == 1)
          {
            fprintf(pspr,"%d, ",primfeld[index2]);
          }
        }
        fprintf(pspr,"\n",0);
      }
      if (tm3 >= dtm)
      {
        fprintf(pspr,"%d: ",index);
        for(index2=1;index2<=anzp;index2++)
        {
          m = modup(primfeld[index2], index);
          if (m != 1)
          {
            fprintf(cnbr,"N%d, ",primfeld[index2]);
          }
        }
        fprintf(cnbr,"\n",0);
      }
    }
  }
  fclose(prim);
  fclose(pspr);
  fclose(cnbr);
  return 0;
}

Nachteil des Programms: Es wird schnell sehr langsam. Außerdem ist es für einige Pseudoprimzahlen blind. So gibt es keine Primzahlen $p<39$, zu denen die 39 pseudoprim ist.