\(
\newcommand{\Mod}[1]{\ \mathrm{mod}\ #1}
\)
Für einen Primzahltest, basierend auf dem kleinen Fermatschen 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. Diese Operation nennt man diskrete Exponentialfunktion oder modulare Exponentiation. Die naive Berechnung, bei der zunächst $a^{p-1}$ vollständig bestimmt und anschließend die Modulo-Operation ausgeführt wird, ist jedoch aussichtslos: sind $a$ und $p$ 10-stellig, dann besitzt $a^p$ schon 100 Milliarden Stellen. Allein die Speicherung einer solchen Zahl würde je nach Darstellung zwischen 40 und 100 GB Speicherplatz benötigen.
Rechenregeln
Die Rechenregeln für Kongruenzen erweisen sich als nützlich beim Modularen Potenzieren. Durch ihre Anwendung wird die effiziente Berechnung, selbst bei großen Zahlen a und p, ermöglicht. 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}. $$
Iteratives Verfahren
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.
Ein $\mathcal{O}(\log_2 n)$-Verfahren
Nutzen wir den Zusammenhang \[ a^n = a^{\frac{n}{2}} \cdot a^{\frac{n}{2}} \quad \text{für gerades } n, \] dann ergibt sich direkt:
Divide & Conquer:
\[
\begin{aligned} a^n \mod N &= \left(a^{\frac{n}{2}} \cdot a^{\frac{n}{2}} \right) \mod N \\
&= \left( a^{\frac{n}{2}} \mod N \right) \cdot \left( a^{\frac{n}{2}} \mod N \right) \mod N \\
&= \left( \left(a^{\frac{n}{2}} \mod N \right)^2 \right) \mod N
\end{aligned} \]
Die Aufgabe wurde in zwei identische Teilprobleme halber Länge zerlegt – typisches Prinzip von Divide & Conquer.
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;
}
In der hier gezeigten Implementierung wird der Exponent bei geradem Wert halbiert. Bei ungeradem Exponenten wird er zunächst um 1 reduziert und im folgenden rekursiven Schritt halbiert. Pro Rekursionsschritt fallen nur konstant viele Operationen an (Quadrieren und ggf. eine zusätzliche Multiplikation). Formal lässt sich der Aufwand z.B. durch
\[ T(n)= \begin{cases} T(n/2)+O(1), & n \text{ gerade},\\ T(n-1)+O(1), & n \text{ ungerade} \end{cases} \]
Die Fallunterscheidung macht explizit, wie der rekursive Schritt bei geraden bzw. ungeraden Exponenten ausgeführt wird. Nach spätestens zwei rekursiven Schritten ist der Exponent mindestens halbiert, damit liegt der Gesamtaufwand in der Größenordnung $\mathcal{O}(\log_2 n)$ und das ist gegenüber der ersten Version eine wesentliche Verbesserung.
Vergleiche
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
| \(\boldsymbol{a}\) | \(\boldsymbol{N}\) | \(\boldsymbol{a^{N-1} \bmod N}\) | \(\boldsymbol{\frac{N}{\log_2(N)}}\) | \(\boldsymbol{\log_2(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.
| \(\boldsymbol{a}\) | \(\boldsymbol{N}\) | \(\boldsymbol{a^{N-1} \bmod N}\) | \(\boldsymbol{\frac{N}{\log_2(N)}}\) | \(\boldsymbol{\log_2(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.
| \(\boldsymbol{a}\) | \(\boldsymbol{N}\) | \(\boldsymbol{a^{N-1} \bmod N}\) | \(\boldsymbol{\frac{N}{\log_2(N)}}\) | \(\boldsymbol{\log_2(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).
Optimierte rekursive Variante
Die oben gezeigte rekursive Version ist bereits sehr schnell. Im ungeraden Fall kann sie jedoch noch etwas verbessert werden: Statt den Exponenten zunächst um 1 zu reduzieren, wird auch hier direkt mit der halbierten Exponentenlänge gearbeitet. Danach wird das Zwischenergebnis quadriert und nur bei ungeradem Exponenten einmal zusätzlich mit der Basis multipliziert.
Dadurch wird pro Rekursionsebene nur ein rekursiver Aufruf benötigt. Die Aufwandsabschätzung passt dann direkt zur Rekursion \[ T(n) = T\big(\lfloor n/2 \rfloor\big) + O(1) \] und der Aufwand bleibt von der Größenordnung \(\mathcal{O}(\log_2 n)\).
LongInt modular_expo_faster2(const LongInt &base,
const LongInt &exponent,
const LongInt &modulus)
{
if (modulus == 1) return 0;
if (exponent == 0) return 1;
if (exponent == 1) return base % modulus;
LongInt base_mod = base % modulus;
LongInt half = modular_expo_faster2(base_mod, exponent / 2, modulus);
LongInt result = (half * half) % modulus;
if (exponent % 2 != 0) {
result = (result * base_mod) % modulus;
}
return result;
}
Umkehrproblem: diskreter Logarithmus
Bemerkenswert ist, dass die Umkehrung der modularen Exponentiation deutlich schwieriger ist. Während sich \(a^x \bmod N\) effizient berechnen lässt, führt das Umkehrproblem \[ a^x \equiv b \pmod N \] auf den diskreten Logarithmus. Für geeignet gewählte Parameter ist dafür im Allgemeinen kein entsprechend effizientes Verfahren bekannt. Diese Schwierigkeit ist eine Grundlage vieler kryptographischer Verfahren.
SageMath-Script
Im eingebetteten SageMath-Beispiel werden kleinere Zahlenbereiche verwendet, damit die Laufzeit in SageCell praktikabel bleibt.
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.