Forum: Knobelecke

Hier also mal die Lösung:

#include <stdio.h>

static const long long m = 10000000000LL;

int main(void)
{
        int i;
        double res = 1;

        for (i = 0; i < 7830457; i++) {
                res *= 2;
                res = (long long)res % m;

        }

        res *= 28433.;
        res = (long long)res % m;
        res += 1.;

        printf("%lli\n", (long long)res);

        return 0;
}

res enthält das (Zwischen-)Ergebnis der Rechnung; zunächst wird (2^7830457)%10000000000 berechnet. Dabei nutzen wir aus, dass (a*b)%m das Gleiche ist wie ((a%m)*(b%m))%m. Somit bleibt das Zwischenergebnis immer in einem Bereich, der von einem double noch dargestellt werden kann.

Nach der Schleife kommt dann der triviale Rest: Multiplikation mit 28433 und die abschliessende Addition mit 1.

Vor der Ausgabe wird das Ergebnis nochmal in einen integer-Typ gecastet, damit da keine Nachkommastellen ausgegeben werden.


Der "rätselhafte" Teil bei diesem Schritt bestand darin, dass (a*b)%m das Gleiche ist wie ((a%m)*(b%m))%m. Das ist immer dann nützlich, wenn man Berechnungen durchführt, die den Wertebereich des zugrundeliegenden Datentyps überschreiten, man aber sowieso nur an einem Teil (und insbesondere dem "hinteren" Teil) des Ergebnisses interessiert ist.


Auf meinem Rechner (schäbige 1.4GHz) dauert diese Berechnung eine knappe Sekunde. Daher meine ursprüngliche Frage, wie sich die Berechnung optimieren lässt.

Dazu erstmal ein leicht nachvollziehbares Beispiel: angenommen, wir sollen 3^128 ausrechnen. Das können wir entweder (wie oben) stur in einer Schleife machen und das Zwischenergebnis immer mit 3 multiplizieren, oder wir stellen folgende Überlegung an:

3^128 = 3^64 * 3^64
3^64 = 3^32 * 3^32
...
3^2 = 3^1 * 3^1

Das lässt sich ebenfalls als Schleife darstellen, jedoch deutlich effizienter, da sich bei jedem Schleifendurchlauf die Anzahl der benötigten Schritte halbiert statt nur um 1 verringert. Die Laufzeit ist von O(n) auf O(log n) gefallen. Wenn der Exponent eine Zweierpotenz ist, geht das auch alles schön auf. Für alle anderen Fälle müssen wir unsere Formulierung etwas verfeinern:

x^n = x^(n/2) * x^(n/2), falls n gerade
x^n = x^(n/2) * x^(n/2) * x, falls n ungerade

Ich habe 3 generös durch x ersetzt um zu zeigen, dass die Wahl der Basis dabei keine Rolle spielt. Die Division (n/2) soll dabei natürlich ganzzahlig sein mit der gewohnten C-Semantik.

In Pseudocode sieht das dann ungefähr so aus:

// computes a^b
pow(a, b)
{
        if (b == 0) return 1;

        if (b%2 == 0) return mpow(a*a, b/2)
        if (b%2 == 1) return a * mpow(a*a, b/2)
}

Angewendet auf unser ursprüngliches Beispiel von der Berechnung der letzten 10 Stellen der Primzahl könnte eine Lösung also so aussehen:

#include <stdio.h>
#include <math.h>

static const long long m = 10000000000LL;

// computes (a^b) % m
double mpow(double a, int b, double m)
{
	if (b < 0)  return 0;
	if (b == 0) return 1;
	if (a == 0) return 0;

	// aa = (a^2) % m
	double aa = fmod((a*a), m);

	if (b%2 == 0) return fmod(    mpow(aa, b/2, m), m);
	if (b%2 == 1) return fmod(a * mpow(aa, b/2, m), m);
}

int main(void)
{
	double res = mpow(2, 7830457, m);
	res = fmod(res*28433, m);
	res += 1.;

	printf("%lli\n", (long long)res);

	return 0;
}

mpow() tut den Trick: es berechnet (a^b)%m mit der angegebenen Optimierung. Und wie vorhergesagt braucht mein treuer Rechner nun nur noch knapp 0.001 Sekunden.

Statt der vorher benötigten 7830457 Iterationen brauchen wir nunmehr nur noch maximal 23 rekursive Aufrufe. In diesem Fall ergibt das eine Beschleunigung um den Faktor ~340454. Bam!


Was es hier mitzunehmen gab war der Modulo-Trick, den ich vorschnellerweise als bekannt vorausgesetzt hatte, und den Potenzen-Trick, der sich aus den guten alten Potenzgesetzen aus der 8. Klasse ableitet:

Potenzen werden multipliziert, indem man die Basis beibehält und die Exponenten addiert.
Potenzen werden dividiert, indem man die Basis beibehält und die Exponenten subtrahiert.
Potenzen werden potenziert, indem man die Basis beibehält und die Exponenten multipliziert.

Story am Rande: der Schreiber dieser Zeilen hatte in der Schule versehentlich mal im Rahmen einer inoffiziellen Schwammschlacht während der kleinen Pause seiner unvermutet in den Klassenraum eintretenden Mathelehrerin einen nassen Schwamm an den Kopf geworfen und wurde damit belohnt, obige Potenzgesetze 10 Mal abschreiben zu dürfen. Ich weiß nicht mehr viel Stoff aus meiner Schulzeit, aber die Potenzgesetze sind unauslöschlich in mein Gedächtnis eingebrannt.


War das Rätsel zu schwer? War euch der hier vorgestellte Mathekram völlig unbekannt? Jetzt wo ihr die Lösung kennt: hätte man darauf kommen können, auch ohne die Tricks vorher zu kennen?

Gruss,
Philip

--
The C Programming Quiz - bitte Fragen einreichen :)