Algorithmen © 2000-2004 Prof. Dr. Alfred Schreiber

5.4
Das Sieb des Eratosthenes


 

  1. Eratosthenes' Sieb im Handbetrieb
  2. Umsetzung in ein Mathematica-Programm

1. Eratosthenes' Sieb im Handbetrieb

Von dem griechischen Universalgelehrten Eratosthenes von Kyrene, der im 3. Jahrhundert v. Chr. der berühmten Bibliothek von Alexandria vorstand, stammt eine interessante Methode, mit der aus einer vorgegebenen Menge ganzer Zahlen die Primzahlen herausgesiebt werden. Es ist überaus empfehlenswert, diesen Algorithmus an einem Beispiel zu beschreiben und dabei gleich einmal "von Hand" durchzurechnen.

Wir wollen uns hier mit einer kleinen Zahlenmenge, etwa den natürlichen Zahlen von 1 bis 20, begnügen.

  1. Notiere alle Zahlen von 1 bis 20:
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  2. Streiche die Zahl 1 (nicht prim):
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  3. Markiere die Zahl 2:
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  4. Streiche sämtliche echten Vielfachen von 2:
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  5. Markiere die erste Zahl, die noch nicht markiert ist (hier: 3):
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  6. Streiche sämtliche echten Vielfachen von 3:
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  7. Markiere die erste Zahl, die noch nicht markiert ist (hier: 5):
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  8. Streiche sämtliche echten Vielfachen von 5 (ohne Wirkung, da 10, 15, 20 schon gestrichen sind):
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  9. Fahre sinngemäß solange fort, bis jede der Zahlen markiert oder gestrichen ist.
        1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
  10. Die markierten Zahlen sind alle prim.

Eine umgangssprachliche Fassung könnte etwa so aussehen:

  1. Setze p = 2
  2. Streiche alle Vielfachen von p
  3. Setze p = erste nicht gestrichene Zahl > p
  4. Gehe nach 2.

Die Arbeitsweise dieses Sieb-Verfahrens lässt sich auf einer interaktiven Webseite mit eigenen Wertvortgaben anschaulich nachvollziehen.


 

2. Umsetzung in ein Mathematica-Programm

Einige Vorüberlegungen sind von Nutzen:

  1. Das Verfahren benötigt eine Tabelle (Liste) {1, ..., n}, in der alle zusammengesetzten Zahlen gestrichen werden. "Streichen" soll im folgenden bedeuten: "durch 0 ersetzen". Als Ergebnis wird die Liste der Primzahlen aus {1, ..., n} mit einer entsprechenden Anzahl von Nullen zurückgegeben.
  2. Der erste Siebdurchgang erfolgt mit der Primzahl 2; die erste zu streichende Zahl ist demnach 4. Als nächstes wird mit 3 gesiebt, und ist die nächste zu streichende Zahl. Wurde mit allen Primzahlen kleiner als p schon gesiebt, so sind alle (mit k < p) bereits gestrichen (k enthält einen Primfaktor < p); also ist die nächste zu streichende Zahl.
  3. In der Tabelle {1, ..., n} braucht nur mit den Primzahlen <= gesiebt zu werden.

Der eigentliche Algorithmus lässt sich übersichtlicher aufschreiben, wenn man das Sieben mit einer festen Primzahl p als eigenen Modul ausgliedert. Das Argument ztab steht hier für eine Tabelle aus Zahlen 1, ..., n:

5.4.1. Mathematica

sieb[p_,ztab_]:=
   Module[{n=Length[ztab],zt=ztab,i=p*p},
      While[i<=n,zt[[i]]=0;If[p>2,i=i+p+p,i=i+p]];
      Return[zt]
   ]

Um Mehrfachstreichungen zu vermeiden, sollte für p > 2 der Index i jedesmal gleich um 2p erhöht werden, denn die Zahlen sind sämtlich gerade (und wurden daher beim Sieben mit p = 2 schon vorher gestrichen!).

Damit können wir nun schon sieben:

ztab=Table[k,{k,1,30}];
ztab=sieb[2,ztab]
{1,2,3,0,5,0,7,0,9,0,11,0,13,0,15,0,17,0,19,0,21,0,23,0,25,0,27,0,29,0}

Die erste nicht gestrichene Zahl > 2 ist 3:

ztab=sieb[3,ztab]
{1,2,3,0,5,0,7,0,0,0,11,0,13,0,0,0,17,0,19,0,0,0,23,0,25,0,0,0,29,0}

... und so weiter.

Das Sieb des Eratosthenes wiederholt diesen Vorgang für alle Primzahlen p <= und lässt sich daher durch die nachstehende Funktion definieren:

5.4.2. Sieb des Eratosthenes (1. Fassung)

eratosthenes[n_]:=
  Module[{p=2,i=2,
     ztab=Table[k,{k,1,n}],m},
     ztab[[1]]=0;
     m=introot[n];
     While[p<=m,
        ztab=sieb[p,ztab];i++;
        While[ztab[[i]]==0,i++];
        p=ztab[[i]]];
     Return[ztab]
  ]
eratosthenes[30]
{0,2,3,0,5,0,7,0,0,0,11,0,13,0,0,0,17,0,19,0,0,0,23,0,0,0,0,0,29,0}

Zunächst wird die Tabelle initialisiert, die Zahl 1 gestrichen und die Schranke für die Siebzahlen mittels introot berechnet. In der While-Schleife werden dann folgende Schritte ausgeführt (solange die Siebzahl p die Schranke m nicht überschritten hat): 1. siebe mit p und rücke weiter; 2. überspringe die gestrichenen Zahlen; 3. hole die nächste noch nicht gestrichene Zahl.

Mit dem Select-Befehl können nun noch alle Elemente, die ungleich 0 sind, herausgefischt werden:

Select[%,# !=0 &]
{2,3,5,7,11,13,17,19,23,29}

Erläuterung: # steht für das erste Argument einer reinen, d.h. unbenannten Funktion; das Zeichen &  schließt den Rumpf dieser Funktion ab.

Man kann das Sieb-Verfahren noch etwas verbessern, indem von vornherein nur ungerade Zahlen betrachtet werden. In diesem Fall benötigen wir ein etwas abgeändertes sieb-Modul. Der Index (die Platznummer in der Liste) i von ist dabei nicht mehr , sondern .

5.4.3. Mathematica

sieb2[p_,ztab_]:=
    Module[{n=Length[ztab],zt=ztab,i},
       i=Quotient[p*p+1,2];
       While[i<=n,zt[[i]]=0;i=i+p];
       Return[zt]
    ]

Die Funktion sieb2 übernimmt das Sieben in folgendem neuen Hauptmodul (worin das Aussondern der Nullen aus der Ergebnisliste gleich mit aufgenommen ist):

5.4.4. Sieb des Eratosthenes (Endfassung)

eratsieb[n_]:=
    Module[{p=3,i=2, m=introot[n],ztab},
       ztab=Table[2*k-1,{k,1,Quotient[n+1,2]}];
       ztab[[1]]=2;
       While[p<=m,
          ztab=sieb2[p,ztab];i++;
          While[ztab[[i]]==0,i++];
          p=ztab[[i]]];
       Return[Select[ztab,# !=0 &]]
    ]

Hier ein Probelauf bis zur Zahl 1000 mit Angabe der Rechenzeit (CPU Pentium II / 200 MHz):

eratsieb[1000] // Timing
{0.05 Second,
{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,
 61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,
 149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,
 229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,
 313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,
 409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,
 499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,
 601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,
 691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,
 809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,
 907,911,919, 929,937,941,947,953,967,971,977,983,991,997}}

 

Stand: 28.01.2004