Schnelle Rechenoperationen mit großen zahlen

[masm]

96Bits == 3 DWORDs
Das Carry flag entspricht dem Übertrag (Grundschulmathematik )

ADC x,0 beachtet den möglichen Übertrag aus niederwertigeren Stellen.

[SpR]

Ahh danke tut mir leid du hast recht mit den 3 DWORDS = 96 bits
Irgendwie ist mir da ein Fehler unterlaufen^^. Das Ausgeben der Zahlen als hexadezimale Zahlen klappt auch sehr gut allerdings würde ich noch gerne wissen wie man die Zahlen dezimal ausgibt hat das jemand eine Idee?

[SpR]

SpR schrieb:

Ahh danke tut mir leid du hast recht mit den 3 DWORDS = 96 bits
Irgendwie ist mir da ein Fehler unterlaufen^^. Nur verstehe ich deinen Algorithmus nicht ganz, in welchem du ausgerechnet hast wieviele bits für ich für 30 Dezimal Stellen benötige, ich habe mal versucht deinen mit deinem Algorithmus auf dein Ergebnis zu kommen leider scheine ich das Prinzip noh nicht ganz raus zu haben hier der Code:

unsigned long x=(unsigned long)(ceil(log(1E30)/log(2))+31);
x=x&-32;

Leider bekomme ich bei diesem Programm die Zahl 128 raus,welche nichts mit 96 Bits und auch nichts mit 3 DWORDS zutun hat.
Das Ausgeben der Zahlen als hexadezimale Zahlen klappt auch sehr gut allerdings würde ich noch gerne wissen wie man die Zahlen dezimal ausgibt hat das jemand eine Idee?

[rkhb]

SpR schrieb:

Leider bekomme ich bei diesem Programm die Zahl 128 raus,welche nichts mit 96 Bits und auch nichts mit 3 DWORDS zutun hat.

Ich habe mal auf folgender Seite die Bits gezählt:

http://manderc.manderby.com/concepts/umrechner/index.php

Für 1E30, also für eine 1 mit 30 Nullen, braucht man 100 Bits, also 4 DWords zu je 32 Bits.

Die Formel lautet also richtig:

int dwords = (int) ceil(log(1E30)/log(2)/32);

Das Ausgeben der Zahlen als hexadezimale Zahlen klappt auch sehr gut allerdings würde ich noch gerne wissen wie man die Zahlen dezimal ausgibt hat das jemand eine Idee?

Beherrscht Du die Umwandlung Integer zu Dezimal (http://dcla.rkhb.de/umwandlung/int2dez.html)? Bei der Division durch 10 musst Du auf die Grundlagen der Division zurückgreifen: http://dcla.rkhb.de/div.html. Ein Beispiel für zwei Dword-Register gibts hier: http://dcla.rkhb.de/div64.html. Die Erweiterung auf 4 Dwords überlasse ich Dir.

Ich frage mich, ob man dafür auch die 128-Bit-SSE-Register nutzbar machen könnte. Weiteres weiß ich aber nicht.

viele grüße
ralph

[SpR]

Vielen dank Ralph das reicht vollkommen aus ich werde mich sofort einmal reinlesen.
Mfg SpR

[rkhb]

SpR schrieb:

Vielen dank Ralph

Bitte, gern geschehen. Danke für die Rückmeldung.

das reicht vollkommen aus ich werde mich sofort einmal reinlesen.

Erst spitz machen, und dann "genug" schreien? Nix da :

#include <stdio.h>
#include <string.h>

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

void bigint_add (void* big1, void* big2)
{
	_asm
	{
		mov edi, big1
		mov esi, big2

		mov eax, [esi+0]
		mov edx, [esi+4]
		mov ecx, [esi+8]
		mov ebx, [esi+12]

		add [edi+0], eax
		adc [edi+4], edx
		adc [edi+8], ecx
		adc [edi+12], ebx
	}
}

void bigint_sub (void* big1, void* big2)
{
	_asm
	{
		mov edi, big1
		mov esi, big2

		mov eax, [esi+0]
		mov edx, [esi+4]
		mov ecx, [esi+8]
		mov ebx, [esi+12]

		sub [edi+0], eax
		sbb [edi+4], edx
		sbb [edi+8], ecx
		sbb [edi+12], ebx
	}
}

void bigint_inc (void* big)
{
	_asm
	{
		mov edi, big

		mov eax, 1
		xor ecx, ecx

		add [edi+0], eax
		adc [edi+4], ecx
		adc [edi+8], ecx
		adc [edi+12], ecx
	}
}

void bigint_dec (void* big)
{
	_asm
	{
		mov edi, big

		mov eax, 1
		xor ecx, ecx

		sub [edi+0], eax
		sbb [edi+4], ecx
		sbb [edi+8], ecx
		sbb [edi+12], ecx
	}
}

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
	}
}

__declspec(naked) void bigint_isnull () // Übergabe ESI: Zeiger auf BigInt
{
	_asm
	{
		mov eax, [esi]
		or eax, [esi+4]
		or eax, [esi+8]
		or eax, [esi+12]
		test eax, eax
		ret
	}
}

int bigint_is_null (void* big)
{
	_asm
	{
		mov esi, big
		mov eax, [esi]
		or eax, [esi+4]
		or eax, [esi+8]
		or eax, [esi+12]
		test eax, eax
		setz al
		movzx eax, al
	}
}

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

	_asm								// http://dcla.rkhb.de/umwandlung/int2dez.html
	{
		xor cl, cl
		lea esi, big1

		Schleife_1:
		push esi
		call bigint_div10
		add esp, 4
		push eax
		add cl, 1
		call bigint_isnull
		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;
	char dez[40] = "buf";

	big1.ll.ll = 0x4674edea40000000; big1.ll.big = 0x0000000c9f2c9cd0; // 1e30
	printf ("hex\t0x%016I64X%016I64X\n",big1.ll.big,big1.ll.ll);
	bigint2dez (&big1, dez);
	printf ("dez\t%s\n",dez);

	bigint_inc (&big1);
	bigint2dez (&big1, dez);
	printf ("inc\t%s\n",dez);

	big2.ll.ll = 123456; big2.ll.big=0;
	bigint_add (&big1, &big2);
	bigint2dez (&big1, dez);
	printf ("add\t%s\n",dez);

	bigint_dec (&big1);
	bigint2dez (&big1, dez);
	printf ("dec\t%s\n",dez);

	big2.ll.ll = 123457; big2.ll.big=0;
	bigint_sub (&big1, &big2);
	bigint2dez (&big1, dez);
	printf ("sub\t%s\n",dez);

	bigint_inc (&big1);
	bigint2dez (&big1, dez);
	printf ("inc\t%s\n",dez);
	if (bigint_is_null (&big1)) puts ("Fehler: Falsches Null");

	bigint_sub (&big1, &big1);
	bigint2dez (&big1, dez);
	printf ("sub\t%s\n",dez);
	if (bigint_is_null (&big1)) puts ("Nu is aber Null");

	return 0;
}

An dem Programm gibts noch einiges zu optimieren und zu lernen. Du könntest Dir auch überlegen, wo Du drehen musst, um mit 50-stelligen Zahlen zu arbeiten. Da ich die Routinen auch für mich selbst nutze, würde ich mich freuen, wenn ich auf Fehler aufmerksam gemacht werde.

viele grüße
ralph

[masm]

rkhb schrieb:

Da ich die Routinen auch für mich selbst nutze, würde ich mich freuen, wenn ich auf Fehler aufmerksam gemacht werde.

Sollten die Routinen unter Windows benutz werden, müssen ESI, EDI und EBX gesichert werden (Win/IntelABI)

[rkhb]

masm schrieb:

Sollten die Routinen unter Windows benutz werden, müssen ESI, EDI und EBX gesichert werden (Win/IntelABI)

Danke für den Hinweis.

Netterweise macht das mein Microsoft-Compiler automatisch. Der PellesC-Compiler macht es aber nicht automatisch. Ich habe mir hin- und herüberlegt, ob ich die Register sichern soll oder nicht. Meiner Meinung nach sind diese Register aber nur innerhalb der Windows- und C-Routinen gefährdet. Ich bewege mich hier aber außerhalb davon. Mit einem Beispiel, wo die Veränderung von ESI, EDI und EBX zur Katastrophe führt, wäre mir sehr geholfen.

viele grüße
ralph

[masm]

Problem gibt es nur bei Callback Funktionen wie z.B. Fensterprozedur, Enumeratoren, ...
Was passiert, hängt von der Windows Version ab. Neuere Versionen scheinen wohl etwas toleranter zu sein, trotzallem sollte man sich einfach an die ABI halten da ansonsten die Potentielle Absturzgefahr besteht

[SpR]

Ich werde darauf achten allerdings hätte ich jetzt noch ein Problem ich berechne die Maximale Dezimalzahl von 4096 bits allerdings kann ich schlecht überprüfen ob die dezimalzahl stimmt, kennt jemand eine Seite oder kann seine Ergebnisse mit meinen abgleichen?
Vielen Dank nocheinmal für die gute Hilfe dank euch verstehe ich jetzt nicht nur wie man mit diesen großen Zahlen rechnet sondern auch die Hintergründe.
MfG SpR

[volkard]

SpR schrieb:

Ich werde darauf achten allerdings hätte ich jetzt noch ein Problem ich berechne die Maximale Dezimalzahl von 4096 bits allerdings kann ich schlecht überprüfen ob die dezimalzahl stimmt, kennt jemand eine Seite oder kann seine Ergebnisse mit meinen abgleichen?
Vielen Dank nocheinmal für die gute Hilfe dank euch verstehe ich jetzt nicht nur wie man mit diesen großen Zahlen rechnet sondern auch die Hintergründe.
MfG SpR

Also 2^4096 brauchste?
104438888141315250669175271071662438257996424904738378038423
348328395390797155745684882681193499755834089010671443926283
798757343818579360726323608785136527794595697654370999834036
159013438371831442807001185594622637631883939771274567233468
434458661749680790870580370407128404874011860911446797778359
802900668693897688178778594690563019026094059957945343282346
930302669644305902501597239986771421554169383555988529148631
823791443449673408781187263949647510018904134900841706167509
366833385055103297208826955076998361636941193301521379682583
718809183365675122131849284636812555022599830041234478486259
567449219461702380650591324561082573183538008760862210283427
019769820231316901767800667519548507992163641937028537512478
401490715913545998279051339961155179427110683113409058427288
427979155484978295432353451706522326906139490598769300212296
339568778287894844061600741294567491982305057164237715481632
138063104590291613692670834285644073044789997190178146576347
322385026725305989979599609079946920177462481771844986745565
925017832907047311943316555080756822184657174637329688491281
952031745700244092661691087414838507841192980452298185733897
764810312608590300130241346718972667321649151113160292078173
8033436090243804708340403154190336

[rkhb]

volkard schrieb:

SpR schrieb:

ich berechne die Maximale Dezimalzahl von 4096 bits allerdings kann ich schlecht überprüfen ob die dezimalzahl stimmt, kennt jemand eine Seite oder kann seine Ergebnisse mit meinen abgleichen?

Also 2^4096 brauchste?
1044...336

Das ist die Anzahl möglicher Dezimalzahlen (ab 0). Die maximale Dezimalzahl, die mit 4096 Bits dargestellt werden kann, ist eins niedriger.

1 Bit -> 2^1-1 -> 1
8 Bit -> 2^8-1 -> 255
16 Bit -> 2^16-1 -> 65535
4096 Bit -> 2^4096-1 -> 1044...335

@SpR:

Die bereits gelinkte Seite http://manderc.manderby.com/concepts/umrechner/index.php schluckt beliebig große Eingaben (bis Dein Speicher ausgeht). Gib mal 4096 Einsen ein. Es dauert ein paar Sekunden, bis das Programm reagiert.

viele grüße
ralph

[SpR]

Vielen Dank werde mal anfangen es abzugleichen^^

[SpR]

Volkard du hast vergsessen dein Ergebnis -1 zu nehmen da das erste bit nur 1 wertig ist. Bis auf die letzte Ziffer stimmt mein Programm mit deinem überein. Meine letzte Ziffer ist 5 und deine 6 ich glaube es liegt dadran dass es nicht 2^4096 ist sondern (2^4096)-1.
MfG SpR außerdem auch dir danke Ralph.

[SpR]

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.

[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

[SpR]

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.

[SpR]

Okay vielen Dank schöne Pfingsten noch.
MfG SpR