\(
\newcommand{\Mod}[1]{\ \mathrm{mod}\ #1}
\)
Für einen Primzahltest (kleiner Fermatscher Satz) wird die Auswertung von
$$
a^{p-1} \Mod{p} , \qquad 0< a < p , \quad a,p \in \mathbb{N}
$$ für große Zahlen $a$ und $p$ benötigt.
Die naive Berechnung, bestimme erst $a^{p-1}$ und dann die Restklasse, ist aussichtslos: sind $a$ und $p$ 10-stellig, dann besitzt $a^p$ schon 100 Milliarden Stellen. Zum Vergleich: die Anzahl der Atome im Universum wird auf $10^{85}$ geschätzt (nach diversen Quellen im Netz).
Die Rechenregeln für Kongruenzen helfen hier weiter. Für die Addition modulo $N$ gilt
\begin{equation}
(a + b) \Mod{N} = ( (a \Mod{N}) + (b \Mod{N}) ) \Mod{N}\end{equation} und genauso hat man für die Multiplikation
\begin{equation} \label{modmult}
(a * b) \Mod{N} = ( (a \Mod{N}) * (b \Mod{N}) ) \Mod{N} .
\end{equation} Es gilt also $$a^{n} \Mod{N} = (a \, a^{n-1}) \Mod{N} = ( (a \Mod{N}) \, (a^{n-1} \Mod{N}) ) \Mod{N}. $$ Die wiederholte Anwendung von \eqref{modmult} zur Bestimmung von $a^n \Mod{m}$ ergibt das iterative Verfahren:
typedef long long LongInt;
// berechnet a^n mod m
LongInt modular_expo(const LongInt &base, const LongInt &exponent,const LongInt &modulus){
LongInt result = 1;
for (LongInt k=0; k < exponent; ++k)
result = ( result * base) % modulus;
return result;
}
Der Aufwand[1] ist $\mathcal{O}(n)$. Diese Komplexität wäre i. Allg. für viele (andere) Verfahren schon sehr gut. Für das vorliegende Problem geht es aber noch besser.
Aus $a^n = a^{\frac{n}{2}} \, a^{\frac{n}{2}}$ für gerades $n$ folgt sofort
\begin{eqnarray*}
a^n \Mod{N} & = & (a^{\frac{n}{2}} \, a^{\frac{n}{2}})\Mod{N} \\
& = & ( ( a^{\frac{n}{2}} \Mod{N}) * ( a^{\frac{n}{2}} \Mod{N}) ) \Mod{N} .
\end{eqnarray*}
Die Aufgabe ist jetzt durch zwei (identische) Probleme halber Länge ersetzt worden. Der Aufwand beträgt nur noch $\mathcal{O}(n/2)$:
if ( exponent % 2 ==0) {
LongInt zw = modular_expo(base, exponent/2,modulus);
result = (zw*zw) % modulus;
}
Das Vorgehen kann rekursiv fortgesetzt werden.
LongInt modular_expo_faster(const LongInt &base, const LongInt &exponent,const LongInt &modulus){
if ( exponent == 2) return (base*base) % modulus;
if ( exponent == 0) return 1;
if ( modulus == 1) return 0;
LongInt result = 1;
if ( exponent % 2 ==0) { // exponent gerade
LongInt zw = modular_expo_faster(base, exponent/2,modulus) % modulus;
result = (zw*zw) % modulus;
}
else { // exponent ungerade
result = (base % modulus);
result = (result * ( modular_expo_faster(base, exponent-1,modulus) % modulus) ) % modulus;
}
return result;
}
Der Aufwand für diesen Algorithmus ist von der Größenordnung $\mathcal{O}(\log_2 N)$. Das ist gegenüber der ersten Version eine wesentliche Verbesserung.
Die nachfolgenden Tabellen zeigen, welchen enormen Vorteil bzgl. der Rechenzeit die „divide and conquer“-Variante hat. Für die Beispiele aus den Tabellen wurden zufällig Zahlen aus dem Bereich $(2,10^6)$ für die Basis $a$ und aus unterschiedlichen Intervallen ungerade Zahlen für den Exponenten $N$ gewählt. Ein möglicher integer overflow wird in beiden Algorithmen nicht berücksichtigt.
$0< N < 10^{8}$: der zweite Algorithmus ist 100000-fach(!) schneller
a | N | a^{N-1} mod N | N/log2(N) | log2(N) | time1 [s] | time2 [s] | speedup |
---|---|---|---|---|---|---|---|
877154 | 22570019 | 1 | 923944.14 | 24.43 | 0.30 | 1.64e-06 | 1.85e+05 |
125791 | 88727699 | 18308659 | 3360530.93 | 26.40 | 1.12 | 1.70e-06 | 6.57e+05 |
467305 | 50750945 | 5658675 | 1982696.45 | 25.60 | 0.65 | 1.70e-06 | 3.80e+05 |
324550 | 54588773 | 48406499 | 2123903.14 | 25.70 | 0.69 | 1.77e-06 | 3.90e+05 |
545898 | 37471999 | 16407209 | 1489389.00 | 25.16 | 0.47 | 2.55e-06 | 1.86e+05 |
Mit 22570019 wird hier eine Primzahl[2] „gefunden“.
Sind für $N$ Zahlen bis $10^{9}$ zugelassen, dann ist das rekursive Verfahren sogar um den Faktor $10^6$ schneller.
a | N | a^{N-1} mod N | N/log2(N) | log2(N) | time1 [s] | time2 [s] | speedup |
---|---|---|---|---|---|---|---|
938109 | 422731183 | 288375025 | 14752355.42 | 28.66 | 5.42 | 2.09e-06 | 2.59e+06 |
761798 | 178398297 | 107859217 | 6508386.34 | 27.41 | 2.27 | 1.75e-06 | 1.30e+06 |
464304 | 289374179 | 10691741 | 10294950.30 | 28.11 | 3.66 | 2.11e-06 | 1.73e+06 |
578052 | 883364873 | 1 | 29724475.60 | 29.72 | 11.28 | 1.92e-06 | 5.86e+06 |
612409 | 562081813 | 137882531 | 19337984.27 | 29.07 | 7.17 | 1.76e-06 | 4.08e+06 |
Für $N<10^{10}$ benötigt eine Auswertung des iterativen Verfahrens schon mehr als eine Minute Rechenzeit.
a | N | a^{N-1} mod N | N/log2(N) | log2(N) | time1 [s] | time2 [s] | speedup |
---|---|---|---|---|---|---|---|
28985 | 4775379709 | 1864993043 | 148520648.66 | 32.15 | 60.79 | 2.17e-06 | 2.81e+07 |
17020 | 6586984007 | 5103444642 | 201949593.26 | 32.62 | 83.84 | 2.11e-06 | 3.98e+07 |
17495 | 6775105577 | 1868531966 | 207458795.14 | 32.66 | 86.25 | 2.60e-06 | 3.32e+07 |
Alle Testrechnungen wurden auf einem PC mit i7-5820K, 3.30GHz, CPU durchgeführt (g++ 8.2, C++, 1 Thread).
Links zu dem Thema:
- Modular exponentiation (wikipedia)
- Kleiner Fermatscher Satz (wikipedia)
- Binäre Exponentiation, Square-and-Multiply (wikipedia)
- Restklassenring
- integer overflow (stackoverflow)
-
- [1]Anzahl der Operationen $*$ und $\%$.
- [2]Ist das Ergebnis $1$, dann ist $N$ ein Kandidat für eine Primzahl. Aus $a^{N-1}\Mod{N} \not=1$ folgt, dass $N$ keine Primzahl ist.