Modulares Potenzieren

Von liebauProgrammierung

\(
\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:

    1. [1]Anzahl der Operationen $*$ und $\%$.
  1. [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.