Frage zu Algorithmus zum Multiplizieren großer Zahlen

drakon

(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)
   -------
     2450

Ich 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.

@ 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

drakon

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

Icematix

Warum malloc'st du? Nimm halt new, dann sparst du dir auch Errorchecking

Habs grad mal gemessen:
Threshold 2, Bits 2^21

Mit malloc: 7,4 bis 7,6 sec
Mit new : 7,8 bis 8,0 sec

Grad 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)

volkard

Was ich nicht verstehe, ist, warum überhaupt viele Allokationen gemacht werden. Weiß man denn nicht vorher schon, wieviel Zwischenspeicher gebraucht wird?

volkard

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?

volkard

Kannst Du bestimmen, wieviel für sagen wie mal length=1024 und kara_threshold=64 nötig ist?

Ich habs jetzt so gelöst (scheint zu funktionieren):

uint_t tmp (length);
uint_t size(length);
while(tmp>Grenze)
{
    tmp/=2; //simuliert Rekursion :D
    size += tmp;
}
nrs* buffer((nrs*)malloc(sizeof(nrs)*size*3)); //*3 weil 3 mal die Funktion in der Funktion aufgerufen wird.

volkard

GroßeZahlen schrieb:

uint_t tmp (length);
uint_t size(length);
while(tmp>Grenze)
{
    tmp/=2; //simuliert Rekursion :D
    size += tmp;
}

Die Schleife macht also
size = length/2 + length/4 + length/8 + ... + 2Grenze
?
Das wäre wenn ich recht rate gleich zu
size = length - 2Grenze

GroßeZahlen schrieb:

nrs* buffer((nrs*)malloc(sizeof(nrs)*size*3)); //*3 weil 3 mal die Funktion in der Funktion aufgerufen wird.

Jo, klingt schlau.

Ich muss diesen Thread nochmal aus der Versenkung holen. Es geht um den "Standard" Multiplikationsalgorithmus (siehe unten). Der verhält sich bei mir sehr komisch. Nochmal zum Konzept: Das Array number (unsigned int number[]) ist Member der Klasse Zahl und repräsentiert eine Zahl. Wenn ich zwei solcher "Zahlen" multiplizieren will, funktioniert das NUR, wenn die Anzahl der Elemente der arrays (array_size) eine Potenz von zwei ist.

//number und a.number sind gleich lang, result ist doppelt so lang
void mult_school(const Zahl& a, Zahl& result) const
{
	std::memset(result.number, 0, array_size*2*sizeof(unsigned int)); //result auf 0 setzen

	for(unsigned int i=0;i!=array_size;++i)
	{
		unsigned int carry = 0;

		for(unsigned int j=0;j!=array_size;++j)
		{
			unsigned int index = i+j; //position, ab der "result" bearbeitet wird

			//-----------------------------------------
			// Errechne number[i]*a.number[j]
			// Splitte das Ergebnis in part1 und part2
			//-----------------------------------------
			unsigned int ahi = number[i] >> 16;
			unsigned int alo = number[i] & 0xFFFF; 
			unsigned int bhi = a.number[j] >> 16;
			unsigned int blo = a.number[j] & 0xFFFF;

			unsigned int part1 = alo * blo;
			unsigned int part2 = ahi * bhi;

			unsigned int middle1 = ahi * blo;
			unsigned int tmp = alo * bhi;
			middle1 += tmp;
			unsigned int middle2 = middle1 < tmp;
			middle2 <<= 16;
			middle2 |= (middle1>>16);
			middle1 <<= 16;

			part1 += middle1;
			unsigned int cyt = part1 < middle1;
			part2 += cyt;
			part2 += middle2;
			/*
			unsigned long long eins = (unsigned long long)number[i] * a.number[j];
			unsigned long long zwei = part1 + ((unsigned long long)part2<<32);
			if(eins!=zwei) std::cerr<<"FEHLER"<<std::endl;
			*/

			//-----------------------------------------
			// Verrechne part1 und part2 mit result
			//-----------------------------------------
			result.number[index] += part1;
			unsigned int cy1 = result.number[index] < part1; //overflow bei addition

			result.number[index+1]+=cy1; //overflow vom Schritt hiervor mit dem nächsten element verrechnen
			unsigned int cy2 = result.number[index+1] < cy1;

			result.number[index+1]+=part2; //nächstes element mir part2 verrechnen
			unsigned int cy3 = result.number[index+1] < part2;

			result.number[index+1]+=carry; //nächstes element mit übertrag vom letzten schleifendurchlauf verrechnen

			carry = result.number[index+1] < carry; //übertrag für nächsten schleifendurchlauf festsetzen 
			carry += (cy3|cy2);
		}
	}
}

Im Prinzip multipliziere ich das ite Element von this mit dem jten Element von a und speichere das Ergebnis in part1 und part2. Das Funktioniert immer und ist auch immer richtig. Anschließend verrechne ich das mit result.
Der Algorithmus an sich scheint richtig zu sein, da alle Arraylängen, die eine Potenz von 2 sind, das richtige Ergebnis produzieren. Wenn die Arraylänge von this und a keine Potenz von 2 ist, erhalte ich solche Fehler:

Array Size = 3;
Aufgabe:      45675872532487774540918482702 * 62469334509438103030645792939
Richtige Lösung: 2853341360242434495933484427440446062475097909763545241178
Erzeuge Lösung:  2853341360242434495933484427440446062457389175715234054144

Auffällig ist, dass die oberen paar bits immer richtig sind, alle weiteren Bits aber auf 0 stehen (nochmal: mit einer Arraysize, die eine Potenz von 2 ist, ist das Ergebnis immer richtig)

Zur Aufgabe oben:
Richtige Lösung: 11101000101111001000110000101111110111111000110101010001101010111001000001011101011111101110100101000001100011110101011011111011111010111000010000011111100101011000011011101000001011001011010
Erzeugte Lösung: 1110100010111100100011000010111111011111100011010101000110101011100100000101110101111110111010010100000110001111010101101111101[b]0000000000000000000000000000000000000000000000000000000000000000[/b]

Viel extremer ist es bei größeren Zahlen. Dann sind die oberen hundert bits richtig gesetzt und dann kommen 500 nullen.

Kann sich jemand vorstellen, woran sowas liegen kann? Provoziere ich irgendwo undefiniertes Verhalten?