Pseudoprimzahlen: Programme

Carmichael-Zahlen Bearbeiten

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  eine quadratfreie, aus mindestens drei Primfaktoren bestehende, natürliche Zahl der Form   ist, so dass für jedes   gilt, dass   die Zahl   teilt. Daraus folgt, dass man sich um den kleinen Fermatschen Satz   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 Bearbeiten

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 Bearbeiten

Die Carmichel-Zahlen nach Chernick sind nur ein kleiner Ausschnitt aus den Carmichael-Zahlen. Sie lassen sich durch folgende Formel generieren:  . Es gibt dabei allerdings eine Einschränkung:

  • Jeder der drei Faktoren  ,   und   muss eine Primzahl sein. Halt, hier gibt es eine Ausnahme: Wenn ein einziger der drei Faktoren selbst eine Carmichael-Zahl ist, so ist   trotzdem eine Carmichael-Zahl.

Da alle Carmichael-Zahlen quadratfrei sind, ist dann auch jeder der drei Faktoren  ,   und   quadratfrei - auch wenn einer davon eine Carmichael-Zahl ist. Da die drei Faktoren paarweise teilerfremd sind, ist auch der Ausdruck  quadratfrei.

Alle gültigen   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 Bearbeiten

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.

Modulare Arithmetik Bearbeiten

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   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   nur Zahlen bis   auftreten, ist das Ergebnis einer Multiplikation vor der Modulooperation immer kleiner als  .

Beispiel:  . Letztendlich interessiert uns ja aber nicht  , sondern vielmehr   mit sagen wir mal  . Das ist  . Nun lässt sich   auch so berechnen:

 

Dabei wird der Wert 7 niemals überschritten.

/* REXX-Programm */
say 'Only a Library!'
exit 1
/* */
/* */
m_u: procedure
  arg a,l,m
/* initialisierung */
  as = a
  ai = a
  li = (l-1)
  DO li
    a = a * ai
    a = a // m
  END
return a

Binäres Potenzieren Bearbeiten

Für Primzahltests nach dem kleinen fermatsche Satz müssen große Potenzen berechnet werden; dabei wäre es extrem zeitaufwendig, xk 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   wird das Ergebnis jeder Multiplikation mod   reduziert, wie im vorherigen Absatz beschrieben.

Insgesamt sind mit dieser Methode nur maximal 2 Multiplikationen pro Bit des Exponenten notwendig, also  .

Liste zusammengesetzter Zahlen Bearbeiten

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 Hilfmittel 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 (n² + 2k*n) von n² bis zur Obergrenze in die Menge auf
  • Sortiere die Elemente der Menge aufsteigend

Implementierung in Python Bearbeiten

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.

Fermatsche Pseudoprimzahlen Bearbeiten

Bei Primzahlen ist im Gegensatz zu fermatschen Pseudoprimzahlen und Carmichael-Zahlen die Kongruenz   für alle Basen   mit   erfüllt. Carmichael-Zahlen unterscheiden sich von anderen fermatschen Pseudoprimzahlen dadurch, dass der Anteil der Primbasen   mit  , für die   ist, größer als 97% ist. Der Einfachheit halber benutzt man als Basen   nur die Primzahlen, die kleiner als die zu testende Zahl   sind. Das Programm testet also zu jeder Primbasis  , ob eine Zahl   die Bedingung   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   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   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  , zu denen die 39 pseudoprim ist.