5. Arithmetische Algorithmen

5.3 Multiplikation

Schriftliche Multiplikation - die Schulmethode

Wie multipliziert man zwei natürliche Zahlen? Die "schriftliche Multiplikation" haben wir alle in der Grundschule gelernt. Daher ist Ihnen vielleicht gar nicht klar: die schriftliche Multiplikation ist ein Algorithmus! Und wahrscheinlich der erste "richtige", der Ihnen im Leben überhaupt begegnet ist. Man kann sogar noch weiter gehen und unsere Dezimaldarstellung als Datenstruktur zur Darstellung von natürlichen Zahlen bezeichnen, die sich gegen die römischen Zahlen) durchgesetzt hat, unter Anderem weil sie effiziente und vergleichsweise einfache Multiplikationsalgorithmen unterstützt. Rufen wir uns diesen "Uralgorithmus" ins Gedächtnis.

Der Schulalgorithmus zur schriftlichen Multiplikation benötigt eine Subroutine (bzw. eine Helfer-Funktion), die eine beliebig lange Zahl mit einer einstelligen Zahl multipliziert. Nennen wir diese Funktion integerTimesDigit

Sobald wir dies verstanden haben, können wir beliebig lange Zahlen multiplizieren:

Strenggenommen müssen wir noch zwei Subroutinen definieren: eine zum Addieren beliebig langer Zahlen (aber das überlasse ich Ihnen), und eine zur Multiplikation zweier einstelliger Zahlen. Für letzteres gibt es keinen Trick. Das müssen Sie auswendig lernen (unter dem Namen kleines Einmaleins). Machen Sie sich klar, dass der gerade beschriebene Algorithmus unabhängig von der Basis unseres Zahlensystems ist - nur das kleine Einmaleins ist für jede Basis anders. Für Basis 2 ist der obige Algorithmus besonders einfach.

Laufzeitanalyse

Wie sieht es jetzt nun mit der Laufzeit der Schulmethode aus? Um eine \(m\)-stellige Zahl mit einer \(n\)-stelligen zu multiplizieren, müssen wir die Subroutine integerTimesDigit insgesamt \(n\)-mal aufrufen. Dies hat eine Laufzeit von \(O(mn)\). Das Ergebnis ist eine Liste von \(m\) Zahlen

\begin{align*} a_0, a_1, a_2, \dots, a_{m-1}\ , \end{align*}

die jeweils wieviele Stellen haben? Jede höchstens \(m+1\); aber sie sind ja "verschoben"; wir haben oben einfach die Nullen am Ende weggelassen. Schreiben wir sie nun der Korrektheit halber hin:

Die längste Zahl, \(a_{m-1}\), kann also bis zu \(n+1 + m-1 = n + m\) Stellen haben. Wir müssen nun noch \(n+1\) Additionen ausführen. Die Zahlen bleiben dabei durchweg höchstens \(n+m\)-stellig (überprüfen Sie es: \(999\times 99\) hat 5 Stellen), was insgesamt \(O(n(n+m))\) Schritt benötigt.

Die Asymmetrie zwischen \(n\) und \(m\) im Ausdruck \(O(n(n+m))\) sollte uns stutzig machen. Wenn wir \(n\) sehr viel größer ist als \(m\), sagen wir \(n = m^2\), dann ist dies \(O(m^2 (m^2 + m)) = O(m^2)\); würden wir die Zahlen vertauschen, dann hätten wir eine Laufzeit von \(O(m(n+m)) = O(m (m^2 + m) = O(m^3))\), also deutlich schneller. Wir können natürlich eine "Wrapper"-Funktion schreiben, die sicherstellt, dass die links stehende Zahl die längere ist (so wie im Beispiel \(24719 \times 8907\). Es geht aber auch anders:

Divide-and-Conquer

Erinnern Sie sich, was Mergesort so einfach und erfolgreich gemacht hat: Wir haben das unsortierte Array in zwei (fast) gleichgroße Teilarrays zerlegt; jedes dieser zwei rekursiv sortiert; dann die beiden sortieren Teilarrays mittels der recht einfachen Funktion merge zu einem ganzen sortieren Array zusammengefügt. Dieses Prinzip Zerlegen - rekursiv lösen - zusammenfügen nennt man in der Algorithmik Divide and Conquer. Können wir für Multiplikation einen Divide-and-Conquer-Algorithmus entwickeln? Wie zerlegt man "eine natürliche Zahl in zwei gleichgroße Hälften"? Damit meine ich nicht \(8000 = 4000 + 4000\), denn da haben die "Hälften" ja gleich viele Stellen wie die ursprüngliche Zahl; nein, ich will eine \(n\)-stellige Zahl irgendwie sinnvoll in zwei \(n/2\)-stellige zerlegen. In unserer Dezimaldarstellung ist eine natürliche Zahl ja einfach ein Array aus Ziffern, 73036417 ist also irgendwie das Array \([7, 3, 0, 3, 6, 4, 1, 7]\). Das zerlegen wir jetzt in \([7,3,0,3]\) und \([6,4,1,7]\). Arithmetisch geht folgendes vor sich:

\begin{align*} 73036417 = 7303 \cdot 10^4 + 6417 \end{align*}

Wie verwenden wir dies nun, um zwei natürliche Zahlen zu multiplizieren? Um unser Gehirn zu schonen und damit wir uns auf die wirklich wichtigen Dinge konzentrieren können, gehe ich im folgenden davon aus, dass beide Zahlen \(n\)-stellig sind und \(n\) auch noch gerade ist. Um nun also zwei \(n\)-stellige Zahlen \(X, Y\) zu multiplizieren, schreiben wir

\begin{align*} X & = A \cdot 10^{n/2} + B \\ Y & = C \cdot 10^{n/2} + D \end{align*}

und schauen dann, was bei der Multiplikation geschieht:

\begin{align*} X \cdot Y & = (A \cdot 10^{n/2} + B) + (C \cdot 10^{n/2} + D) \\ & = AC \cdot 10^n + (AD + BC) \cdot 10^{n/2} + CD \end{align*}

An einem konkreten Beispiel:

\begin{align*} 73,036,417 \cdot 43,826,247 & = 7303 \cdot 4382 \cdot 10^8 + (7303\cdot 6247 + 6417 \cdot 4382) \cdot 10^4 + 6417\cdot 6247 \\ & = 46,863,351 \cdot 10^8 + 73,741,135 \cdot 10^4 + 40,086,999 \\ & = 3,200,912,051,436,999 \end{align*}

Hier ist nun unser Pseudo-Code für den naiven Divide-and-Conquer-Algorithmus:

def naiveMult(x, y):
	# x, y are two n-bit integers, and n is even 	
	a = first(x)
	b = secondHalf(x)
	c = firstHalf(y)
	d = secondHalf(y)

	ac = naiveMult (a, c)
	ad = naiveMult (a, d)
	bc = naiveMult (b, c)
	bd = naiveMult (b, d)

	result = ac * 10**n + (ad + bc) * 10**(n/2) + cd                         
  return result

Da fehlt natürlich noch die Rekursionsbasis, beispielsweise wenn \(n=1\) ist; da können wir nichts mehr zerlegen sondern wenden einfach das kleine Einmaleins an. Um die Laufzeit zu analysieren, beachten Sie, dass die Zeieln 3,4,5,6 und 13 jeweils \(O(n)\) Schritte benötigen. Lassen Sie sich nicht von den * und 10**n in Zeile 13 einschüchtern: hier wird nicht wirklich etwas multipliziert; mit \(10^n\) multiplizieren heißt ja einfach, hinten \(n\) Nullen anzuhängen, und das geht in \(O(n)\) Zeit.

Die Zahl \(k\) hängt von der genauen Implementierung ab und davon, was wir als "Schritt" verstehen. Ein Takt auf dem Computer beispielsweise. Aber genau das versuchen wir ja mit der \(O\)-Notation unter den Teppich zu kehren. Um nun aber die Laufzeit von naiveMult zu analysieren, müssen wir dieses \(k\) wieder unter dem Teppich hervorholen.

Sei \(T(n)\) also die Anzahl der Schritte, die naiveMult im Worst-Case macht, wenn es zwei \(n\)-stellige Zahlen als Eingabe bekommt. Wir gehen von nun an auch davon aus, dass \(n\) eine Zweierpotenz ist, dass also alle Zahlenarrays sich in zwei gleich große Hälften unterteilen lassen, bis sie irgendwann einstellig sind und der Basisfall (das kleine Einmaleins) erreicht ist. Also: \(n = 2^d\). Was ist \(T(n)\)? Es gilt

\begin{align} T(n) & \leq kn + 4 T(n/2) \label{T-rec} \end{align}

Warum? Der Term \(kn\) ist die Anzahl der Schritte, die in den Zeilen 3, 4, 5, 6, 13 getan wird. Der Term \(4 T(n/2)\) steht für die Anzahl der Schritte, die durch die rekursiven Aufrufe in Zeilen 8, 9, 10, 11 verursacht werden. Was ist \(T(1)\)? Im Basisfall wenden wir das kleine Einmaleins an. Wenn wir beispielsweise als Basis 256 wählen, dann können wir in dem Fall einfach zwei Bytes multiplizieren. Die Logik dafür ist auf dem Prozessorchip fest implementiert; die Multiplikation braucht also nur einen Schritt. Da aber das ganze drum herum (der Test \(n=1\), der Funktionsaufruf etc.) auch etwas braucht, sagen wir einfach: \(T(1) \leq k\). Das gleiche \(k\) wie oben. Wenn wir \(k\) richtig wählen, dann sind wir auf der sicheren Seite. Ersetzen wir nun das \(\leq\) in (\ref{T-rec}) durch ein \(=\), definieren uns also eine Funktion \(S(n)\) per

\begin{align*} S(n) := \begin{cases} k & \textnormal{ if $n=1$,} \\ kn + 4 S(n/2) & \textnormal{ else.} \end{cases} \end{align*}

Es sollte klar sein, dass \(T(n) \leq S(n)\) gilt. Um für \(S(n)\) eine geschlossene Form zu berechnen (und somit die Laufzeit von naiveMult zu verstehen), müssen wir den Rekursionsbaum zeichnen. Ich nehme, wie oben bereits erwähnt, an, dass \(n\) eine Zweierpotenz ist.

Ein Knoten auf Ebene \(i\) repräsentiert also einen rekursiven Aufruf von naiveMult mit \(n/2^i\)-stelligen Zahlen. Er verursacht also \(kn/2^i\) Schritte. Ebene \(i\) trägt also insgesamt

\begin{align*} 4^i kn / 2^i = kn2^i \end{align*}

zum Gesamtwert \(S(n)\) bei. Wir sehen: die Ebenen werden immer "teurer", im Gegensatz zu Mergesort, wo jede Ebene ungefähr gleichviel Arbeit verursacht hat. Was ist die letzte Ebene? Erinnern wir uns: wir nehmen an, dass \(n = 2^d\) ist. Dann haben wir Ebenen \(0, 1, \dots, d\), und Ebene \(d\) trägt \(kn 2^d\) bei. Daher: für \(n = 2^d\) gilt

\begin{align*} S(n) & = kn (1 + 2 + 4 + 8 + \cdots + 2^d) \\ & = kn (2^{d+1} - 1) \\ & = kn (2n - 1) \ . \end{align*}

Wir haben also einen ganz gehörigen Aufwand betrieben und am Ende doch wieder nur einen \(O(n^2)\)-Algorithmus bekommen, auch nicht besser als die Schulmethode.

Die Karatsuba-Multiplikation

Umsonst war der Aufwand doch nicht, weil er die Bühne bereitet für einen tatsächlich viel schnelleren Algorithmus: die Karatsuba-Multiplikation. Der fast schon magische Trick ist, die vier Multiplikationen in Zeilen 8, 9, 10, 11 durch drei Multiplikationen (und mehrere Additionen) zu ersetzen. Die zentrale Beobachtung ist

\begin{align*} (A+B) (C+D) = AC + AD + BC + BD \end{align*} und somit \begin{align} AD + BC = (A+B) (C+D) - AC - BD \label{karatsuba-trick} \end{align}

Wir berechnen nun also rekursiv die Produkte \(AC\), \(BD\) und \((A+B)(C+D)\). Mittels (\ref{karatsuba-trick}) erhalten wir nun \(AD + BC\) durch eine einfache (und billige) Subtraktion. Der Pseudo-Code ist hier:

def karatsuba (x, y):
	# x, y are two n-bit integers, and n is even 	
	a = first(x)
	b = secondHalf(x)
	c = firstHalf(y)
	d = secondHalf(y)

	ac = karatsuba (a, c)
	bd = karatsuba (b, d)
	a_plus_b_times_c_plus_d = karatsuba (a+b, c+d)

	ad_plus_bc = a_plus_b_times_c_plus_d - ac - bd 

	result = ac * 10**n + ad_plus_bc * 10**(n/2) + cd  ac * 10**n + (ad + bc) * 10**(n/2) + cd
	return result

Die Implementierung ist nun etwas schwieriger, weil wir zusätzlich Subtraktion als Subroutine benötigen. Aber auch dies geht in \(O(n)\). Und klar: wenn wir eine Bibliothek für Big-Integer-Arithmetik schreiben, dann wird die sowieso eine Funktion für Subtraktion enthalten. Als Rekursionsbasis legen wir nun den Fall fest, dass die Zahlen \(3\) oder weniger Stellen haben (und da wenden wir dann z.B. naiveMult an); der Basisfall verursacht höchstens \(2k\) Schritte, für passend gewähltes \(k\) (warum Faktor 2, wird weiter unten klar). Für größere \(n\) haben wir die "Netto-Arbeit" aus den Zeilen 3, 4, 5, 6, 12, 14. Da wird nur zerlegt, addiert und subtrahiert, und daher benötigen wir auch hier maximal \(kn\) Schritte, für passend gewähltes \(k\). Unangenehm ist jetzt jedoch, dass der rekursive Aufruf in Zeile 10 es nicht notwenigerweise mit \(n/2\)-stelligen Zahlen zu tun hat; durch die Addition kann nämlich eine Stelle hinzukommen; es sind also maximal \(\ceil{n/2}+1\)-stellige Zahlen. Die rekursiven Aufruf ein Zeile 8 und 9 bekommen \(ceil{n/2}\)-stellige Zahlen. Um uns die Laufzeitabschätzung zu erleichert, sind wir konservativ und sagen: alle Zahlen in den rekursiven Aufrufen sind maximal \(\ceil{n/2}+1\)-stellig. Das stimmt dann auf jeden Fall. Wir können die Laufzeit abschätzen durch

\begin{align*} R(n) := \begin{cases} 3k & \textnormal{ if $n \leq 3$} \\ kn + 3 S(\ceil{n/2}+1) & \textnormal{ else.} \end{cases} \end{align*}

Das geht einfach zu analysieren, wenn \(n\) die spezielle Form \(n = 2^{d}+2\) hat. Dann gilt nämlich \(\ceil{n/2}+1 = 2^{d-1}+2\), und die Anzahl der Stellen in den rekursive Aufrufen hat wiederum die gleiche spezielle Form. Wir zeichnen den Rekursionsbaum:

Auf Ebene \(i\) haben wir also \(3^i\) Knoten, von denen jeder einen Aufruf mit \(2^{d-i} + 2\)-stelligen Zahlen darstellt. Für Ebenen \(i = 0, 1, \dots, d-1\) ist die Rekursionsbasis noch nicht erreicht, und der Beitrag von Ebene \(i\) zu \(R(n)\) ist

\begin{align} 3^i k (2^{d-i} + 2) = k2^d \pfrac{3}{2}^i + 3^i k \label{contrib-layer-i} \end{align}

Für \(i=d\) wird dieser Ausdruck zu \(2k3^d\). Da haben wir aber die Rekursionsbasis bereits erreicht, müssen also anders rechnen! Auf Ebene \(d\) haben wir \(3^d\) Knoten. Jeder davon trägt \(2k\) zu \(S(n)\) bei. Also stimmt der Ausdruck (\ref{contrib-layer-i}) auch für \(i=d\). Wir erhalten \(S(n)\), indem wir über alle Ebenen summieren:

\begin{align} S(2^d + 2) & = \sum_{i=0}^d \left( k2^d \pfrac{3}{2}^i + 3^i k \right) \nonumber \\ & = k 2^d \sum_{i=0}^d \pfrac{3}{2}^i + k \sum_{i=0}^d 3^i \label{S-geom-sum} \end{align}

Wir brauchen nun die Formel für die geometrische Summe:

\begin{align*} 1 + x + x^2 + x^3 + \cdots + x^d = \frac{x^{d+1} - 1}{x-1} \ . \end{align*}

Somit gilt

\begin{align*} (\ref{S-geom-sum}) & = k2^d \frac{\pfrac{3}{2}^{d+1}-1}{\frac{3}{2}-1} + k \frac{3^{d+1} - 1}{2} \\ & \leq 2 k 2^{d} \pfrac{3}{2}^{d+1} + \frac{k}{2} 3^{d+1} \\ & = k 3^{d+1} + \frac{k}{2} 3^{d+1} = \frac{9k}{2} 3^{d} \ . \end{align*}

Bevor wir nun eine schöne Form für \(S(n)\) angeben, können wir aus dem obigen Ausdruck schon etwas hilfreiches ablesen: wenn wir \(n\) ungefähr verdoppeln, also von \(2^d + 2\) auf \(2^{d+1} + 2\), dann verdreifacht sich die Laufzeit. Das deckt sich wunderbar mit unseren empirischen Daten über die Laufzeit der Multiplikation in Python.