Schnelle Rechenoperationen mit großen zahlen

volkard

SpR schrieb:

Bin mittlerweile dabei E = M * c² auszurechnen mit meinem großen Datentyp dabei ist mir aufgefallen ,dass ich noch gar nicht weiß wie man mit großen Zahlen multipliziert. Wäre sehr freundlich, wenn mir jemand die Theorie erklärt oder ein paar Seiten dazu verlinkt.

Lahm: Siehe Grundschule.

Ich würde aber gleich weitergehen und Karatsuba machen. Ist nichtmal schwieriger, aber schon ganz ordentlich schnell. Die noch schnelleren werden ab hier schlagartig sehr schwierig.

rkhb

SpR schrieb:

Bin mittlerweile dabei E = M * c² auszurechnen mit meinem großen Datentyp dabei ist mir aufgefallen ,dass ich noch gar nicht weiß wie man mit großen Zahlen multipliziert. Wäre sehr freundlich, wenn mir jemand die Theorie erklärt oder ein paar Seiten dazu verlinkt.

Das hast Du schon in der Grundschule gelernt, nur mit Dezimalzahlen. Rechne mal "auf dem Papier" aus:

1234 * 5678
-----------
       9872
      8638
     7404
    6170
===========
    7006652

Eigentlich hast Du die Zahlen auseinandergerissen und mit einzelnen Ziffern gerechnet. Die Zahl 5678 ist also eigentlich nur eine zusammengeklebte Folge der Dezimalzahlen '5', '6', '7', '8', '9'. Nun stelle Dir vor, da wären keine Dezimalzahlen (Basis 10), sondern Dwords (Basis 2^32). Dann ständen da oben nicht zwei Dezimalzahlen mit 4 Dezimalziffern, sondern zwei BigInts mit 4 Dwords. Nennen wir die erste Dezimalzahl x und die zweite Dezimalzahl y und nummerieren jede Ziffer von rechts ab 0. Abstrakt würde dann aus '1234 * 5678' ein 'x3x2x1x0 * y3y2y1y0'. Die vier Zeilen unter dem Strich hast Du wie folgt ermittelt:

1. Zeile: x * y0
2. Zeile: x * y1 und um eins nach links verschoben
3. Zeile: x * y3 und um zwei nach links verschoben
4. Zeile: x * y4 und um drei nach links verschoben

und alles addiert.

Die einzelnen Multiplikationen kannst Du auf dieselbe Weise aufdröseln. Wenn man sich die Ziffern als Elemente in einem Array vorstellt (z.B. y0 als Element 0, y1 als Element 1, y3 als Element 3 usw.), dann wird die Verschieberei einfach. Die jeweiligen Ergebnisse landen im Ziel am entsprechenden Index. In Zeile 4 z. B. landet die Null im Index 4. Es ist kein Zufall, dass die Operation mit y4 im Index 4 landet. Du verstehst das besser, wenn Du die Rechnung tatsächlich auf einem Stück Karopapier machst, und Dir bei jedem kleinsten Schritt klarmachst, was Du hier eigentlich machst.

Ein Problem ist, dass das Ergebnis sehr viel größer sein kann, als die einzelnen Zahlen. Die Multiplikation zweier 4-stelliger Zahlen kann 8-stellig sein (Addition der Stellen). Die Multiplikation zweier BigInts (mit fester Stellenanzaahl) kann die Grenzen eines einzelnen BigInts sprengen. Diese "Sprengung" kannst Du bei nachfolgendem Programm gut erkennen:

#include <stdio.h>

typedef union
{
	struct
	{
		unsigned int dwords[4];
	} bigint;
	struct
	{
		unsigned long long ll;
		unsigned long long big;
	} ll;
} big_union;

void bigint_mul (void* big1, void* big2, void* erg)
{
	_asm
	{
		; erg auf Null setzen
		mov edi, erg
		xor eax, eax
		mov ecx, 4
		rep stosd

		mov esi, big1
		mov ebx, big2
		mov edi, erg

		; x0 * y0
		mov eax, [esi+0]
		mul dword ptr [ebx+0]
		mov [edi+0], eax
		mov [edi+4], edx

		; x1 * y0
		mov eax, [esi+0]
		mul dword ptr [ebx+4]
		add [edi+4], eax
		adc [edi+8], edx				// Übertrag
		adc [edi+12], 0					// Weiterer Übertrag

		; x2 * y0
		mov eax, [esi+0]
		mul dword ptr [ebx+8]
		add [edi+8], eax
		adc [edi+12], edx

		; x3 * y0
		mov eax, [esi+0]
		mul dword ptr [ebx+12]
		add [edi+12], eax
		; Übertrag wäre außerhalb des BigInt-Wertebereichs

		; x0 * y1
		mov eax, [esi+4]
		mul dword ptr [ebx+0]
		add [edi+4], eax
		adc [edi+8], edx
		adc [edi+12], 0

		; x1 * y1
		mov eax, [esi+4]
		mul dword ptr [ebx+4]
		add [edi+8], eax
		adc [edi+12], edx
		; Weiterer Übertrag wäre außerhalb des BigInt-Wertebereichs

		; x2 * y1
		mov eax, [esi+4]
		mul dword ptr [ebx+8]
		add [edi+12], eax            // Edit: Bug beseitigt
		; Übertrag wäre außerhalb des BigInt-Wertebereichs

		; x3 * y1
		; mov eax, [esi+4]
		; mul dword ptr [ebx+12]
		; Ergebnis ist außerhalb des BigInt-Wertebereichs

		; x0 * y2
		mov eax, [esi+8]
		mul dword ptr [ebx+0]
		add [edi+8], eax
		adc [edi+12], edx
		; Weiterer Übertrag wäre außerhalb des BigInt-Wertebereichs

		; x1 * y2
		mov eax, [esi+8]
		mul dword ptr [ebx+4]
		add [edi+12], eax
		; Übertrag wäre außerhalb des BigInt-Wertebereichs

		; x2 * y2
		mov eax, [esi+8]
		mul dword ptr [ebx+8]
		; Ergebnis ist außerhalb des BigInt-Wertebereichs

		; x3 * y2
		; mov eax, [esi+8]
		; mul dword ptr [ebx+12]
		; Ergebnis ist außerhalb des BigInt-Wertebereichs

		; x0 * y3
		mov eax, [esi+12]
		mul dword ptr [ebx+0]
		add [edi+12], eax
		; Übertrag wäre außerhalb des BigInt-Wertebereichs

		; x1 * y3
		; mov eax, [esi+12]
		; mul dword ptr [ebx+4]
		; Ergebnis ist außerhalb des BigInt-Wertebereichs

		; x2 * y3
		; mov eax, [esi+12]
		; mul dword ptr [ebx+8]
		; Ergebnis ist außerhalb des BigInt-Wertebereichs

		; x3 * y3
		; mov eax, [esi+12]
		; mul dword ptr [ebx+12]
		; Ergebnis ist außerhalb des BigInt-Wertebereichs
	}
}

unsigned int bigint_div10 (void* big)
{
	_asm
	{
		mov esi, big
		mov ebx, 10

		mov eax, [esi+12]
		xor edx, edx
		div ebx
		mov [esi+12], eax

		mov eax, [esi+8]
		div ebx
		mov [esi+8], eax

		mov eax, [esi+4]
		div ebx
		mov [esi+4], eax

		mov eax, [esi+0]
		div ebx
		mov [esi+0], eax

		mov eax, edx					// Rückgabe: Rest
	}
}

void bigint2dez (void* big, char* buf)
{
	big_union big1;

	_asm								// http://dcla.rkhb.de/umwandlung/int2dez.html
	{
		mov esi, big					// big1 = big
		lea edi, big1
		mov ecx, 4
		rep movsd

		xor cl, cl
		lea esi, big1

		Schleife_1:
		push esi
		call bigint_div10
		add esp, 4
		push eax
		add cl, 1
		mov eax, [esi]
		or eax, [esi+4]
		or eax, [esi+8]
		or eax, [esi+12]
		jnz Schleife_1

		mov edi, buf
		Schleife_2:
		pop eax
		or al, 00110000b
		stosb
		loop Schleife_2
		mov byte ptr [edi], 0
	}
}

int main(void)
{
	big_union big1, big2, big3;
	char dez[40] = "buf";
	int i;

	big1.ll.ll = 1000; big1.ll.big=0;
	big2.ll.ll = 1000; big2.ll.big=0;
	bigint2dez (&big1, dez);
	printf ("big1\t0x%016I64X%016I64X\t%s\n",big1.ll.big,big1.ll.ll,dez);

	for (i=0; i<15; ++i)
	{
		bigint_mul (&big1, &big2, &big3);
		bigint2dez (&big3, dez);
		printf ("*1000\t0x%016I64X%016I64X\t%s\n",big3.ll.big,big3.ll.ll,dez);
		big1 = big3;
	}

	return 0;
}

viele grüße
ralph

volkard schrieb:

SpR schrieb:

Bin mittlerweile dabei E = M * c² auszurechnen mit meinem großen Datentyp dabei ist mir aufgefallen ,dass ich noch gar nicht weiß wie man mit großen Zahlen multipliziert. Wäre sehr freundlich, wenn mir jemand die Theorie erklärt oder ein paar Seiten dazu verlinkt.

Lahm: Siehe Grundschule.

Ich würde aber gleich weitergehen und Karatsuba machen. Ist nichtmal schwieriger, aber schon ganz ordentlich schnell. Die noch schnelleren werden ab hier schlagartig sehr schwierig.

Ich danke dir, du hast recht das eigentlich Malrechnen lahm ist, allerings frage ich lieber nach ob e´s einen guten Algorithmus dafür gibt, anstatt irgend nen mist zu proggen^^. Vielen Dank Ralph dein Post hat mir auch sehr geholfen, könnte jemand mal den Namen oder nen Link zu den Algorithmen geben, welcher schneller als Karatsuba sind, heißt nicht das ich sie verstehe würde mir aber gerne einmal das Prinzip ansehen.

volkard

SpR schrieb:

könnte jemand mal den Namen oder nen Link zu den Algorithmen geben, welcher schneller als Karatsuba sind, heißt nicht das ich sie verstehe würde mir aber gerne einmal das Prinzip ansehen.

Wikipedia "Karastuba" verlinkt auf die schnelleren.
Und sagt, daß bei Deinen 4096 Bits die schnelleren noch gar nicht so greifen.

Okay vielen Dank schöne Pfingsten noch.
MfG SpR

nachtfeuer

Ich würde empfehlen, erstmal zu gucken, und zu verstehen, wie die Karatsuba-Multiplikation arbeitet, denn sie läßt sich überraschend angenehm und leicht in Assembler umsetzen.

Darüber hinaus ist es sehr hilfreich, Binärmultiplikation aufzuschreiben und zu verstehen. Also 1*1 ist ja noch einfach. Aber wieviel ist 11 * 11 (binär)? oder
wieviel wäre dann 1111111111111111 * 1111111111111111?

volkard

nachtfeuer schrieb:

Darüber hinaus ist es sehr hilfreich, Binärmultiplikation aufzuschreiben und zu verstehen.

Binär? Wollten wir nicht übungshalber das 256-er-System nehmen und dann das 4294967296-er benutzen?

nachtfeuer

volkard schrieb:

Binär? Wollten wir nicht übungshalber das 256-er-System nehmen und dann das 4294967296-er benutzen?

4294967296-er? Ist das nicht etwas veraltet? - Bei 32Bit Compi mit SSE ginge zumindest ein 2^128-er Zahlensystem.
Aber Avx...256er Zahlensystem als Einstiegsdroge und dann weiter machen mit 2^256er.

Man fragt sich gelegentlich, was eigentlich so ein funktionaler Sprachinterpreter wie z. B. der ghci Haskellinterpreter intern so treibt, wenn man mit so Begriffen wie 2*(2²⁵⁶⁾3 + 4 * (2²⁵⁶⁾2 + 2* (2^256) usw. herumrechnet oder z.B. (2²⁵⁶⁾64 -1.

(aber 256er? da hat wohl einer mit 8 Bit Compies angefangen was?)

Christoph

nachtfeuer schrieb:

Man fragt sich gelegentlich, was eigentlich so ein funktionaler Sprachinterpreter wie z. B. der ghci Haskellinterpreter intern so treibt, wenn man mit so Begriffen wie 2*(2²⁵⁶⁾3 + 4 * (2²⁵⁶⁾2 + 2* (2^256) usw. herumrechnet oder z.B. (2²⁵⁶⁾64 -1.

Im Zweifelsfall steckt oft libgmp dahinter, weil man es selber nur schwer schneller hinbekommt. Bei ghc ist das so.

nachtfeuer

Mir fällt noch ein, es gab Anfang letzten Jahres einen ganz ähnlichen Thread,
( http://www.c-plusplus.net/forum/280689 )
und ein Problem des TOs war glaub ich, neben der fehlenden Erfahrung mit ein bißchen Herumprobieren, dass die Bitmultiplikation/Addition nicht gleich verstanden wurde, dabei ist das Prinzip eigentlich sehr einfach, wenn man es mal auf Papier gemacht hat. Wenn man die Binärmultiplikation/Addition verstanden hat, und das entgegenkommen der Prozessorbefehle, versteht sich auch der Asm-Code gleich viel besser.

Wenn man von großen Zahlen wie (2²⁵⁶⁾64 lediglich eins abziehen will, und das ist nicht ganz trivial, kann man die Zahl von hinten nach vorn nach der nächstbesten binären 1 absuchen. Bei Zahlen wie 2^x steht die 1 natürlich ganz vorne bzw. ganz links. Bei der Zahl (2²⁵⁶⁾64 - 1 steht dagegen die erste 1 ganz hinten bzw. ganz rechts.

Man sucht also von rechts die nächstbeste 1, macht aus dieser (theoretisch) eine Null und aus dem Rest ( den folgenden Nullen) per

not reg

eine 1.
Direkt auf der Registerebene (der Abschnitt mit der 1) macht man (praktisch) einfach ein

dec reg

(dec reg ginge natürlich auch bei den restlichen, hinteren (nuller) Abschnitten, da ja 00 - 1 FF sein muß)
(und ginge so ganz ohne Carryflag (aber mit CX also Cindy, der kleinen Zähltochter vom Graf Zahl...) )