Guter Algorithmus zum multiplizieren großer Zahlen gesucht

[FAlter]

Hallo,

ich bin gerade dabei, darüber nachzudenken, wie ich am gescheitesten zwei 128 Bit große Integerzahlen auf meinem 32 Bit System multiplizieren kann.

Ein Gedanke wäre, die schriftliche Multiplikation aus der Grundschule nachzubilden und die Vorzeichen vorher extra behandeln/hinterher entsprechend wieder korrigieren, sodass man immer positive Zahlen multipliziert. Dabei würde ich 32bittig multiplizieren und die Ergebnisse dann entsprechend addieren.

Nun kommt mir dieser Ansatz aber sehr komplex vor. Gibt es vielleicht einen besseren Ansatz (einfacher, schneller?)

Ein noch schlimmerer und schon längst wieder verworfener Gedanke war es, das ganze so oft mit sich zu addieren wie es der andere Faktor angibt. Aber das ist zu umständlich und wird wohl ewig dauern.

Ideen können z. B. in Worten, Delphi (class operator Multiply(const Left, Right: Int128): Int128 oder Assembler (eax zeigt auf/adressiert Ergebnis, edx auf linken und ecx auf rechten Faktor) formuliert werden.
Int128 ist ein Record, dass zwei Int64 (Lo und Hi) enthält. Addition und Subtraktion, Vergleiche und Negation sind per Operatorüberladung bereits problemlos möglich.

So, jetzt schlaf ich erstmal ne Nacht darüber...

Mfg
FAlter

[Neutral General]

Zitat:

ich bin gerade dabei, darüber nachzudenken, wie ich am gescheitesten zwei 128 Bit große Integerzahlen auf meinem 32 Bit System multiplizieren kann.

Weiter habe ich nicht gelesen. Aber das reicht schon.

Mein Tipp: Schlafen und morgen weiterdenken
Das werde ich jetzt auch tun.

Gute Nacht

[Corpsman]

Ich denke die Variante mit dem Zerlegen in die 4 Int32bit zahlen dürfte das beste sein.

so werden die ja schlieslich intern auch gerechnet ...

Die Schulmethode Shift add wird da wohl zu langsam sein .

[jfheins]

Also ich würde sagen, du kannst die Zahlen einfach ausmultiplizieren

Soll heißen, wenn du die in 2 64bit Zahlen splittest (verschnkt übrigens ein bit, da du 2 Vorzeichenbits hast)

Und dann kannnst du rechnen: A*B = a1*b1*2^128 + a1*b2*2^64 + a2*b1*2^64 + a2*b2

Jetzt musst du nur noch eine Funktion schreiben, die 2 64bit Zahlen multipliziert und das Ergebnis als 128bit Zahl zurückgibt (die 2^64 kannste ja shiften...)

(Mit 32bit wirds natürlich entsprechend länger ...)

Btw.: Auf fast allen modernen CPUs gibt es doch bereits vorgefertigte 128bit Arithmetik, oder nicht?

[himitsu]

Die 64-Bit-Operationen sind aber auch nicht grad schnell (wenn man nicht grad in 64 Bit programmiert)
Dann doch lieber direkt 32 Bit-Operationen(?)
ups, stümmt, da gab's doch noch was (ok, dann halt weiter in ASM und SSE2)

[FAlter]

Hi,

Zitat von jfheins:

Btw.: Auf fast allen modernen CPUs gibt es doch bereits vorgefertigte 128bit Arithmetik, oder nicht?

Und das wäre? Meine kann MMX, SSE und SSE2. Aber das sind alles Themen mit denen ich mich noch nie beschäftigt habe. Nur mal schnell nachgesehen, was das für Anweisungen kennt.

Die MMX-Register sind nur 64 Bit groß. Das ist defginitiv keine 128er Arithmetik. Außerdem, so weit ich es lese, arbeitet das maximal mit 32 bit, und man verarbeitet eben 2 32bittige Zahlen zusammen. Also ohne dass sich der Übertrag der einen z. B. auf die andere auswirkt.

Bei SSE lese ich nur was von Floating Point Operationen, auch da eben mehrere Zahlen glecihzeitig.

Bei SSE2 ist beides möglich, Packed Floats und Packed Integers; und da kann man schonmal endlich mit 64-Bit-Werten rechnen, von denen zwei in dem 128-Bit-Register stehen.

Ich sehe also keine direkte Möglichkeit mit meiner CPU 128-Bittig zu multiplizieren.

Mfg
FAlter

[jfheins]

Naja, kann sein dass es das nicht in einem Befehl gibt ... ich hab nur SSE&128bit Register gelesen

Aber du brauchst ja auch nur 64bit Arithmetik

Denn: A*B = a1*b1*2^128 + a1*b2*2^64 + a2*b1*2^64 + a2*b2

Den ersten Term kannst du vergessen. (Würde deinen 128bit Typen sprengen)

Den 2. kannste ausrechnen mit 64bit und zum dritten dazuaddieren *wuppdi* haste den "hohen-ergbnis-64bit-wert" und der vierte Term liefert dir den niedrigen 64bit wert.

Also keine native 128bit Miltiplikation von Ganzzahlen aber alles was du brauchst um es mit wenigen Assemble Befehlen zu schreiben

[FAlter]

Hi,

aber wie sieht es mit den Überträgen aus, die bei der 64Bit-Multiplikation ggf. entstehen? Ich habe das Gefühl, diese werden wohl verworfen... [edit] Unsinniger Satz... [/edit]

Aber mit den 64 Bit (Quad Words):

http://flatassembler.net/docs.php?article=manual#2.1.16

SSE2 mit ghepackten QWords:

paddq performs the addition of packed quad words, psubq performs the substraction of packed quad words, pmuludq performs an unsigned multiply of low double words from each corresponding quad words and returns the results in packed quad words. These instructions follow the same rules for operands as the general MMX operations described in 2.1.14.

Das heißt die einzige 64-Bit-Multiplikation die es mit SSE2 gibt ist eine 32-Bit-Multiplikation mit 63 64-Bit ergebnissen, ansonsten kann man nur 64-Bit addieren und subtrahieren. Und zwar bei allen drei Anweisungen jeweils zwei Operationen (zwei Zahlen) gleichzeitig.

Also geht das mit der 64-Bit-Multiplikation auch nur über den 32-Bit-Umweg. Da werde ich wohl letztlich es doch per Hand mit 32-Bit stricken. Erstmal ohne SSE und dann wenns klappt und ich auch ne Division (mit Rest!) habe ggf. nochmal zum lernen / als Übung mit SSE.

Mfg
FAlter

[edit] Was mach ich denn alles für Fehler? Ich will meine Rechtschreibprüfung zurück haben! [/edit]

[FAlter]

Hi,

hier mein aktueller Stand, ich habe in Assembler eine Proc zum unsigned multiplizieren geschrieben, die sogar ein 256 bittiges Ergebnis zurückliefert. Und dann In Delphi die Behandlung des Vorzeichens. Also jeweils in der Sprache, wo ich es einfacher fand.

Hier die Unsigned Multiplikation, ist sicherlich noch optimierungsfähig:

Delphi-Quellcode:

procedure MulPInt128U256U(const Result: PInt128; Left, Right: PInt128);
//Result: 256 bit --> zwei PInt128 (Lo, Hi) nacheinander!!!
//Unsigned!
//for 32-bit; nur wenn ss = ds, was bei Delphi immer zu sein scheint
//(kleine Optimierung durch Voraussetzung)
var
  buf1, buf2, buf3: array[0..4] of LongWord;
asm
//mul performs an unsigned multiplication of the operand and the accumulator.
//[...] If the operand is a double word, the processor multiplies it
//by the contents of EAX and returns the 64-bit result in EDX and EAX. mul
//sets CF and OF when the upper half of the result is nonzero, otherwise they
//are cleared. Rules for the operand are the same as for the inc instruction.

  push ebx
  push esi

  push eax
  push ecx

  mov esi, edx

  mov ecx, [ecx]
  mov ebx, eax
  call @@multiply

  mov ebx, [esp] //=original ecx
  mov ecx, [ebx+4]
  lea ebx, [buf1]
  call @@multiply

  mov ebx, [esp]
  mov ecx, [ebx+8]
  lea ebx, [buf2]
  call @@multiply

  mov ebx, [esp]
  mov ecx, [ebx+12]
  lea ebx, [buf3]
  call @@multiply

  //pop ecx
  add esp, 4 //Wert von ecx nicht wichtig, add schneller
  pop eax

  lea ebx, [buf1]
  xor ecx, ecx
  call @@add

  lea ebx, [buf2]
  call @@add

  lea ebx, [buf3]
  call @@add

  jmp @@end

  @@multiply:
    //eine Zeile multiplizieren
    //ecx = dword-faktor
    //abx = buf für ergebnis
    //und esi: Zeigt auf 128-Bit-Faktor!!!

    mov eax, [esi]
    mul ecx
    mov [ebx], eax
    mov [ebx+4], edx

    mov eax, [esi+4]
    mul ecx
    add[ebx+4], eax
    adc edx, 0
    mov [ebx+8], edx
    jc @@m1             //Übertrag speichern
      mov [ebx+12], 0
      jmp @@m2
    @@m1:
      mov [ebx+12], 1
    @@m2:

    mov eax, [esi+8]
    mul ecx
    add [ebx+8], eax
    adc [ebx+12], edx
    jc @@m3             //Übertrag speichern
      mov [ebx+16], 0
      jmp @@m4
    @@m3:
      mov [ebx+16], 1
    @@m4:

    mov eax, [esi+12]
    mul ecx
    add [ebx+12], eax
    adc [ebx+16], edx

  //Wenn jetzt noch ein Bit Übertrag ist, scheiss ich drauf!
  //Oder quäle mich wenn es immer total falsch rechnet :(
  //Meiner Berechnung sollten die 160 Bit (5*32 Bit) aber genau passen für eine
  //Multiplikation 32-Bit mit 128-Bit
  ret

  @@addcarry:
    //eax --> Ergebnis-Stelle
    //ebx ist wieder buf --> Stelle
    //ecx - letzter Übertrag

    add ecx, [ebx]  //Übertrag verrechnen

    jnc @@a1        //neuer Übertrag?
      mov edx, 1    //Speichern
    jmp @@a2
    @@a1:           //kein neuer Übertrag
      xor edx, edx  //Speichern
    @@a2:

    add [eax], ecx  //eigentliche Addition

    jnc @@a3        //neuer Übertrag?
      inc edx       //Speichern/erhöhen
    @@a3:

    mov ecx, edx
  ret

  @@add:
    //EAX -> Ergebnis, letzte Stelle
    //(nicht sie Stelle, die zuerst addiert werden soll, sondern DAVOR)
    //zwei der Zeilen zusammenaddieren
    //ebx ist wieder buf
    //ecx - letzter Übertrag

    add eax, 4  //ne Stelle weiterrutschen

    push eax    //4 Stellen addieren...
      call @@addcarry

      add eax, 4
      add ebx, 4
      call @@addcarry

      add eax, 4
      add ebx, 4
      call @@addcarry

      add eax, 4
      add ebx, 4
      call @@addcarry
    pop eax
  ret

  @@end:

  pop esi
  pop ebx
end;

Und hier die Vorzeichenbehaftete Multiplikation, die darauf zurückgreift:

Delphi-Quellcode:

procedure MulPInt128(const Result: PInt128; Left, Right: PInt128);
var
  Negative: Boolean;
  Help1, Help2: PInt128;
  UseH1, UseH2: Boolean;
  Reslt{$IFOPT Q+}, Reslt2{$ENDIF}: PInt128;
begin
  Help1 := nil; UseH1 := false;
  Help2 := nil; UseH2 := false;

  Negative := false;

  try

    if Left^ < 0 then
    begin
      Negative := true;
      UseH1 := true;
      New(Help1);
      Help1^ := -Left^;
    end
    else
      Help1 := Left;
    if Right^ < 0 then
    begin
      Negative := not Negative;
      UseH2 := true;
      New(Help2);
      Help2^ := -Right^;
    end
    else
      Help2 := Right;

    GetMem(Reslt, SizeOf(Int128) * 2);
    try
      MulPInt128U256U(Reslt, Help1, Help2);

      Result^ := Reslt^;

      {$IFOPT Q+}
      Reslt2 := Reslt;
      inc(Reslt2);

      if (Reslt2^ <> 0) or (Reslt < 0) then
        raise EIntOverflow.Create('Integer Overflow.');
      {$ENDIF}
    finally
      FreeMem(Reslt);
    end;

    if Negative then
      Result^ := -Result^;
  finally
    if UseH2 then Dispose(Help2);
    if UseH1 then Dispose(Help1);
  end;
end;

Mfg
FAlter

[jfheins]

Ich kenne mich jetzt nicht mit Assembler aus, aber ... du addierst doch nicht etwa, oder?

Außerdem kann ich mir nicht vorstellen, dass das effizient ist - das sieht mir einfach zu lang aus

Um nochmal auf meinen Vorschlag zurückzukommen:

128bit - A*B = a1*b2*2^64 + a2*b1*2^64 + a2*b2

Kannst du weiter runterbrechen - d.h. eine Multiplikation von 2 64bit Werten sieht dann so aus:

x*y = x1*y1*2^64 + x1*y2*2^32 + x2*y1*2^32 + x2*y2

Wobei x1 das high-order dword ist und x2 das low-order dword. [s]Außerdem kannst du wieder den ersten Term vernachlässigen (würde wieder zu einem Überlauf führen)[s] Vergiss das ^^

Du musst dann eben Funktionen für 64bit und 32bit Multiplikation sowie zur Multiplikation mit 2er Potenzen zur Verfügung stellen

Du könntest daraus auch etwas rekursives machen nach dem Motto "Teile und Herrsche"