Frage zu Algorithmus zum Multiplizieren großer Zahlen
-
Hallo, ich wollte eine eigene UINT Klasse für große Zahlen schreiben. Zur Zeit bin ich beim Multiplizieren.
Vorab ein paar Festlegungen, um die Sache einfacher zu gestalten (den code portabler schreiben kann man nachher immernoch):
CHAR_BIT = 8
sizeof(unsigned int) = 4
Endianess = Little EndianDie Klasse sieht so aus:
class Zahl { unsigned int number[ /*beliebige Größe*/ ]; public: static const unsigned int array_size = /*die array size*/; zahl& operator*=(unsigned int a); };Ich habe schon einen Multiplikationsalgorithmus implementiert (shift and add) und wollte jetzt mal die "native" Multiplikation nutzen. Die Grundidee ist, jedes unsigned int in 2 teile zu spalten, um den Überlauf behandeln zu können. Allerdings kann ich keine funktionierende Version implementieren. Die Grundidee ist diese:
zahl& operator*=(unsigned int a) { const unsigned int first_half = 0xFFFFu; const unsigned int second_half = 0xFFFF0000u; for(unsigned int i = 0; i != array_size; ++i) { unsigned int tmp1 = (number[i] & first_half) * (a & first_half); unsigned int tmp2 = (number[i] & second_half) * (a & second_half); number[i] = tmp1; //... } return *this; }Aber irgendwie fehlt mir hier der Ansatz, wie es weitergehen soll. Habt ihr eine Idee?
-
diese frage ist zu schwierig für das c++ forum.
-
GroßeZahlen schrieb:
Die Grundidee ist diese:
zahl& operator*=(unsigned int a) { const unsigned int first_half = 0xFFFFu; const unsigned int second_half = 0xFFFF0000u; for(unsigned int i = 0; i != array_size; ++i) { unsigned int tmp1 = (number[i] & first_half) * (a & first_half); unsigned int tmp2 = (number[i] & second_half) * (a & second_half); number[i] = tmp1; //... } return *this; }Aber irgendwie fehlt mir hier der Ansatz, wie es weitergehen soll. Habt ihr eine Idee?
Du zerlegst also zwei Zahlen in jeweils 16-Bit-Hälften:
a = ahi + alo mit ahi = 0x????0000 und alo = 0x0000???? b = bhi + blo mit bhi = 0x????0000 und blo = 0x0000????Allerdings gilt ahi*bhi mod 2^32 = 0. Und das ist das, was Du mit unsigned int (bei 32bit) wirklich ausrechnest. Du solltest vielleicht das Verschieben nicht vergessen:
Alternativ: ahi = a >> 16; alo = a & 0xFFFF; bhi = b >> 16; blo = b & 0xFFFF; ==> a = ahi * 0x10000 + alo mit ahi = 0x0000???? und alo = 0x0000???? b = bhi * 0x10000 + blo mit bhi = 0x0000???? und blo = 0x0000????Dann ist gilt
a*b = (ahi*0x10000+alo) * (bhi*0x1000+blo) = ahi*alo*0x100000000 + (ahi*blo+alo*bhi)*0x10000 + alo*blo ^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^ (1) (2) (3) 0x????????00000000 (1) + 0x0000????????0000 (2) + 0x00000000???????? (3) ----------------------------- 0x???????????????? (a*b)Vielleicht kommst Du ja so weiter...
Auch noch interessant: http://en.wikipedia.org/wiki/Karatsuba_algorithm
-
Den Karatsuba Algorithmus zur Multiplikation von grossen Zahlen habe ich mal implementiert.
Das ist das ganze Projekt, falls sich jemand dafür interessiert:
http://dl.dropbox.com/u/1716912/karatsuba.rarEs ist eine Implementierung einer Big Integer Klasse.
Fokus war allerding schon in der Implementierung des Multiplikationsalgorithmuses und daher ist der Rest nur schnell hingeklatscht und unvollständig.
Ein paar Verbesserungsvorschläge habe ich damals allerdings noch vermerkt.
-
Ich habe den Ansatz von krümelkacker mal verfolgt. Diese Funktion soll A * B rechnen und das Ergebnis in einem unsigned int Array speichern:
void IO_MULT () { //möglichst viel overflow provozieren :) unsigned int A = 0xFFFFFFFFu; unsigned int B = 0xFFFFFFFFu; unsigned int part1; unsigned int part2; unsigned int ergebnis[2]; // ------------------------------- // Vorbereitung: // ------------------------------- unsigned int ahi = A >> 16; unsigned int alo = A & 0xFFFF; unsigned int bhi = B >> 16; unsigned int blo = B & 0xFFFF; // ------------------------------- // Schritt 1 und 3: // ------------------------------- unsigned int part1 = alo * blo; // = 0xFFFE0001 unsigned int part2 = ahi * alo; // = 0xFFFE0001 // ------------------------------- // Schritt 2: // ------------------------------- unsigned int middle1 = ahi*blo + alo*bhi; unsigned int middle2 = 0; if(middle1 < (ahi*blo)) //overflow bei addition? middle2 = 1; middle2 = (middle2 << 16) + (middle1 >> 16); middle1 = middle1 << 16; // ------------------------------- // middle mit part1 und part2 verrechnen // ------------------------------- part1 += middle1; if(part1 < middle1) //overflow bei addition? ++part2; part2 += middle2; // ------------------------------- // Ergebnis ausgeben // ------------------------------- ergebnis[0] = part1; ergebnis[1] = part2; std::cout<< std::hex << ergebnis[0] << ergebnis[2]; }Abgesehen davon, dass der Stil natürlich grausam ist (Deklaration am Anfang der Funktion, unnötige temporäre Variablen usw) erhalte ich hier das falsche Ergebnis. Die Ausgabe ist:
0x1FFFFFFFe <-- Ausgabe 0xFFFFFFFE00000001 <-- richtiges ErgebnisSieht jemand, wo mein Fehler liegt
-
zum einen grefist du bei der ausgabe auf ergebnis[2] zu. das kann nicht funktionieren.
btw müsste es in zeile 24 nicht "ahi * bhi" heißen, statt "ahi * alo" ?
-
fragezeichen? schrieb:
zum einen grefist du bei der ausgabe auf ergebnis[2] zu. das kann nicht funktionieren.
Ja da hast du recht. Dummer Flüchtigkeitsfehler. Lustigerweise ändert das nichts an der Ausgabe Oo
fragezeichen? schrieb:
btw müsste es in zeile 24 nicht "ahi * bhi" heißen, statt "ahi * alo" ?
Ich habs mal geändert (war nur copy&paste von krümelkacker). Am Ergebnis ändert es nichts, da ja in diesem Spezialsfall alo = ahi = blo = bhi = 0xFFFF .
Hier die geänderte Version:Der Fehler lag an einer Falschen Annahme:
Interne Repräsentation: Little Endian
also ergebnis[0] low ordered
und ergebnis[1] high ordered.Bei der ausgabe muss das allerdings umgedreht sein:
std::cout<< ergebnis[1] << ' ' <<ergebnis[0]. Dann stimmts !
-
Nachdem ich jetzt (dank krümelkacker und fragezeichen?) zwei Zahlen multiplizieren und den Overflow abfangen kann, möchte ich das jetzt in meinem Programm implementieren. Nochmal schnell das Grundkonzept:
Ein internes unsigned int array, was eine Zahl repräsentiert soll mit einem unsigned int multipliziert werden. Das ganze sähe dann so aus:zahl& operator*= (unsigned int b) { for(unsigned int i = 0; i != array_size; ++i ) { // HIER :krümelkackers algorithmus, wobei // A = number[i] und B = b ist. number[i] = part1; number[i+1] += part2; if(i < array_size-1) //noch nicht beim letzten array element angekommen if(number[i+1] < part2) //overflow bei addition? for(int j = i+2; j != array_size; ++j) if( ++number[j] ) break; //overflow handeln } return *this; }Das stimmt so aber nicht (wenn ich das so implementiere, dann bekomme ich viel zu große Zahlen als Ergebnis). Ist denn die Grundidee so richtig? Also das Array durchgehen und das aktuelle Element = part1 zu setzen und das nächste element mit part2 zu addieren?
-
@ drakon:
Ich versuche mich gerade am karatsuba algorithmus und versuche dein Projekt auf meine Anwendung umzuschreiben. Dein code sah so aus:my_integer& my_integer::operator *= ( const my_integer& other ) { int digits = get_digits (); if ( get_digits () <= 1 ) //... else { my_integer d1 ( *this, 0, digits/2 ); my_integer d2 ( *this, digits/2,digits ); my_integer d3 ( other, 0, digits/2 ); my_integer d4 ( other, digits/2,digits ); my_integer x = d1*d3; my_integer y = d2*d4; my_integer z = (d1-d2)*(d4-d3); z.shift ( digits/2 ); my_integer x1 = x; x1.shift ( digits ); my_integer x2 = x; x2.shift ( digits/2 ); my_integer y1 = y; y1.shift ( digits/2 ); my_integer y2 = y; (*this) = x1 + x2 + y1 + y2 + z; } //... }Allerdings hänge ich an dieser Stelle:
my_integer x1 = x; x1.shift ( digits ); my_integer x2 = x; x2.shift ( digits/2 ); my_integer y1 = y; y1.shift ( digits/2 ); my_integer y2 = y; (*this) = x1 + x2 + y1 + y2 + z;1. Frage: Müsste x1 nicht auch um (digits/2) geshiftet werden und nicht um (digits) (Zeile 2)?
2. Frage: ist diese Zeile(*this) = x1 + x2 + y1 + y2 + z;nicht gleichbedeutend mit dieser hier?
(*this) = x.shift(digits/2+1) + y.shift(digits/2+1) + z;da (pseudocode)
geg.: x1 = x2 = x x1<<(digits/2) + x2<<(digits/2) = 2*(x<<(digits/2)) = (x<<(digits/2))<<1 = x<<(digits/2+1)?
-
- (schau dir besser mal die Antwort auf deine zweite Frage an, dann wird sich das hier erübrigen..;))
Nein. Ich denke mein Code ist schon korrekt.
Mal ein kleines Beispiel, wie wir 2 2-stellige Zahlen multiplizieren:
(Unterstriche sind reingeshiftete 0en)
98 *25 ------- x1: 18__ // 9*2 x2: 18_ y1: 40_ // 8*5 y2: 40 z: 3_ // (9-8)(5-2) ------- 2450Ich denke hier liegt das Missverständnis. Es ist kein binäres Shift. Alles, was shift bei mir hier macht ist 0en von rechts einzufügen. Ich denke mal, dass dann auch die erste Frage klar ist.
Ich habe da wohl einen ungünstigen Namen gewählt. Für mich war das gleich das, was ich machen wollte, nämlich einfach die Zahlen nach links shiften. Dachte halt nicht, dass jemand anders den Code sehen wird.
Aber zu optimieren gibt es auf jeden Fall noch Sachen (habe ich ja auch vermerkt). Der Code sollte halt einfach die Methode ziemlich direkt implementieren, so dass man ihn gut nach vollziehen kann.
- (schau dir besser mal die Antwort auf deine zweite Frage an, dann wird sich das hier erübrigen..;))
-
@ drakon:
Keine Ahnung was mich bei meinem letzten Post geritten hat (man bechte dir Uhrzeit
)Vielen Dank nochmal, dass du mir deine Implementation des Karatsuba Algorithmus hochgeladen hast! Ich habs jetzt endlich geschafft, das Ganze bei mir zu implementieren, mit sehr guten Ergebnissen:
Integer mit 2^20 (~1 Mio) Bits: Schulmultiplikation: 6,234 sec Karatsuba mit Schulmultiplikation ab 512 Bits: 0,718 sec
-
Hehe. Kein Problem. War wahrscheinlich auch wirklich ungeschickt benannt von mir.
Wie sind denn da die genauen Testumgebeungen gewesen? Also Compiler, Einstellungen usw.
Und wie ist der Code? - Hast du da noch optimiert oder ihn so in etwa übernommen, wie ich ihn hatte?
Auf jeden Fall ist das ein nettes Ergebnis.
-
Compiler ist Visual Studio 2008 Express. Alle Optimierungen eingeschaltet. Code ist hier zu finden:
http://codepad.org/a34smVCX
Ein bisschen abgewichen von deinem Code bin ich schon. Zb. erstelle ich keine Objekte d1-d4, sondern reiche nur Zeiger weiter.Ich denke vor Allem in der school_mult Funktion könnte man noch ein bisschen Speed rausholen...
-
Warum malloc'st du? Nimm halt new, dann sparst du dir auch Errorchecking
-
Habs grad mal gemessen:
Threshold 2, Bits 2^21Mit malloc: 7,4 bis 7,6 sec
Mit new : 7,8 bis 8,0 secGrad laufen auf meinem Computer aber mehrere Prozesse, die wahrscheinlich stören. Werde später nochmal nachmessen.
-
edit: ERROR_HANDLING ist natürlich nicht aktiviert gewesen.
-
So gerade nochmal eine vernünftige Messung gemacht mit diesen Einstellungen:
Microsoft Visual Studio 2008 Express, Komplette Optimierung. Länge der zu multiplizierenden Zahlen: 2^21 Bits ERROR_HANDLING, also manuelles prüfen auf den Rückgabewert von malloc auf NULL nicht, aktiviert. Grenze für Schulmultiplikation einmal auf 1 (möglichst viele Allokationen) und einmal auf 16 (beste Laufzeit): malloc: 18484 (Grenze 1) 2187 (Grenze 16) new: 19483 (Grenze 1) 2218 (Grenze 16) new (std::nothrow): 20469 (Grenze 1) 2219 (Grenze 16)
-
Was ich nicht verstehe, ist, warum überhaupt viele Allokationen gemacht werden. Weiß man denn nicht vorher schon, wieviel Zwischenspeicher gebraucht wird?
-
static inline void mult_kara(uint_type length, const uint_type * a, const uint_type* b, uint_type* res) { uint_type * tmp = malloc/new/egal(4*length+4711*sin(humidityOfMoon); //evtl if(!tmp) throw mult_kara(length, a, b, res, tmp); free/delete tmp; } //im eigentlichen algo ist dann alles viel billiger. static inline void mult_kara(uint_type length, const uint_type * a, const uint_type* b, uint_type* res, uint_type* tmp) { if(length<=kara_threshold) { std::memset(res, 0u, (length<<1u)*sizeof(uint_type)); mult_school(length, a, b, res); return; } uint_type half_l ( length>>1u ); uint_type * y = tmp; tmp += length;//klappt immer mult_kara(half_l, a, b, y, tmp); uint_type * x = tmp; tmp += length;//klappt immer std::memset(x, 0u, sizeof(uint_type)*length); mult_kara(half_l, a+half_l, b+half_l, x, tmp); ... //die free-Orgie entfällt einfach }Kein Witz.
-
Stimmt. Die Idee ist sehr gut. Ich habs gerade mal so umgesetzt. Ergebnisse:
Grenze 16: 2000 ms (vorher 2187) Grenze 1: 4656 ms (vorher 18484)Allerdings weis ich nicht, wie genau ich die Buffersize in Abhängigkeit von der Abbruchgrenze ermitteln soll (um möglichst wenig Speicher zu belegen). Hat da jemand eine Idee?