Schnelle Integer-Multiplikation mit FFT

Hier könnt Ihr gute, von Euch geschriebene Codes posten. Sie müssen auf jeden Fall funktionieren und sollten möglichst effizient, elegant und beispielhaft oder einfach nur cool sein.
Benutzeravatar
Helle
Beiträge: 566
Registriert: 11.11.2004 16:13
Wohnort: Magdeburg

Re: Schnelle Integer-Multiplikation mit FFT

Beitrag von Helle »

Wegen des heutigen Feiertages (Sachsen-Anhalt ist da nicht so :mrgreen: ) war Zeit für die Procedure-Umsetzung. Die eigentliche Berechnung wurde komplett auf ASM (80-Bit-FPU) umgestellt; alles andere wäre Unfug.

Code: Alles auswählen

;- Schnelle Integer-Multiplikation mit FFT für große Zahlen
;- Bis auf die max.Länge keine Einschränkungen mehr bei den Faktoren!
;- Kein Unicode bei Strings!
;- Zur Erhöhung des Wertebereiches wird mit 80-Bit-FPU gerechnet
;- Getestet bis Faktoren-Länge 128MB mit CPU Intel i7-6700K @4.7GHz und RAM 32GB DDR4
;- Berechnungs-Schema läuft ab nach Schmetterlingsgraph
;- Keine Fremd-Bibliotheken, keine globalen Variablen
;- Zur besseren Orientierung im Code ist der Großteil des PB-Codes noch auskommentiert enthalten. Kann natürlich entfernt werden (wie auch die Infos)
;- Windows7/64, PureBasic 5.41 LTS (x64)
;- "Helle" Klaus Helbing, 06.01.2016

Procedure.s FFT_Mul(LFak1, LFak2, PFak1, PFak2)
  EnableASM 
  ;EnableExplicit ;war für Übersicht mal nötig
  Protected.q BufferCos, BufferCosA, BufferDez, BufferDezA, BufferFX, BufferFXA, BufferFaktor1, BufferFaktor1A
  Protected.q BufferFaktor2, BufferFaktor2A, BufferProdukt, BufferProduktA, BufferSin, BufferSinA, N
  Protected.q Exponent, L1, TA_Gesamt, TE_Gesamt, TA_Faktoren, TE_Faktoren, TA_Winkel, TE_Winkel, TA_FFT, TE_FFT
  Protected.q TA_Multi, TE_Multi, TA_Rueck, TE_Rueck, TA_Ergebnis, TE_Ergebnis
  Protected.s sAusgabe, sProdukt, sFaktor, sExpo, sWinkel, sFaktoren, sFFT, sMulti, sRueck, sErgebnis, sGesamt

  If #PB_Compiler_Unicode      ;das jetzt noch für Unicode zu machen schenke ich mir
    sAusgabe = "Achtung! Bei Strings als Faktoren kein Unicode verwenden!" + #LFCR$ + #LFCR$
  EndIf

  ;PF1 = PFak1
  !mov r8,[p.v_PFak1]        ;zur Sicherheit so, sonst r8 direkt (als erstes natürlich) 
  !mov [PF1],r8
  ;PF2 = PFak2
  !mov r9,[p.v_PFak2]        ;zur Sicherheit so, sonst r9 direkt  
  !mov [PF2],r9

  ;L1 = LFak1
  !mov rax,[p.v_LFak1]       ;zur Sicherheit so, sonst rcx direkt  
  !mov [L1],rax
  ;LNull11 = L1               ;echte Länge, auch hier Null-Test möglich
  !mov [LNull11],rax         ;echte Länge, auch hier Null-Test möglich
  ;L2 = LFak2
  !mov rdx,[p.v_LFak2]       ;zur Sicherheit so, sonst rdx direkt  
  !mov [L2],rdx
  ;LNull22 = L2
  !mov [LNull22],rdx

  ;If L1 < L2
  ;  LNull1 = L1              ;hier als Zwischenspeicher
  ;  L1 = L2                  ;größte Länge ermitteln; evtl. hier auf Null testen
  ;  L2 = LNull1
  ;EndIf
  !cmp rax,rdx
  !jae @f
    !mov rcx,rax
    !mov [L1],rdx
    !mov [L2],rcx
  !@@:
  ;LenF = L1
  !mov rdx,[L1]
  !mov [LenF],rdx            ;LenF erstmal setzen, könnte ja schon 2-er Potenz sein 

  !bsr rcx,[L1]              ;höchstes gesetztes Bit
  !bsf rdx,[L2]              ;Test auf Länge Null
  !jz @f
    !mov [Exponent],rcx
    !bsf rdx,[L1]            ;niedrigstes gesetztes Bit
    !jmp .Faktoren_OK        ;local label
  !@@:
  sAusgabe + #LFCR$ + "Fehler: Mindestens einer der Faktoren hat Länge Null!"
 ProcedureReturn sAusgabe

  !.Faktoren_OK:
  !cmp rcx,rdx
  !je @f                     ;also nur 1 Bit (identisch) gesetzt, ist 2-er Potenz
    !inc rcx
    !mov rax,1
    !shl rax,cl
    !mov [LenF],rax          ;LenF ist jetzt 2-er Potenz
    !add [Exponent],1
  !@@:
  !add [Exponent],1

  !mov rcx,3                 ;Produkt-Länge muss mindestens 8 sein, sonst an anderen Stellen Klimmzüge
  !sub rcx,[Exponent]
  !js @f                     ;Exponent ist > 3 (Hauptfall)
  !jz @f                     ;Exponent ist = 3
    !shl [LenF],cl
    !add [Exponent],rcx
  !@@:

  !push [Exponent]           ;die beiden nur für Anzeige, blanke Optik!
  POP Exponent
  !push [L1]
  POP L1  

  !cmp [Exponent],28
  !jbe @f
    sAusgabe + " Hinweis: Diese Größenordnung ist ungetestet! Keine Garantie auf Richtigkeit!"+ #LFCR$  ;Space am Anfang beachten!
  !@@:

  ;N = LenF << 1
  !mov rax,[LenF]
  !shl rax,1
  !mov [N],rax
  !push [N]
  POP N                      ;wird für AllocateMemory benötigt!

  ;LNull1 = LenF - LNull11    ;LNull11=echte Länge. Bei LNull1 fängt hinterstes Byte an 
  !mov rax,[LenF]
  !sub rax,[LNull11]
  !mov [LNull1],rax
  ;LNull2 = LenF - LNull22
  !mov rax,[LenF]
  !sub rax,[LNull22]
  !mov [LNull2],rax
  ;L1hinten = LenF >> 1       ;LenF ist 2-er Potenz
  ;L2hinten = LenF >> 1
  !mov rax,[LenF]
  !shr rax,1
  !mov [L1hinten],rax
  !mov [L2hinten],rax
  ;If LNull11 < L1hinten
    ;L1hinten = LNull11
  ;EndIf
  !cmp rax,[LNull11]
  !jbe @f  
    !mov rax,[LNull11]
    !mov [L1hinten],rax
  !@@:
  ;L1vorn = LNull11 - L1hinten
  !mov rax,[LNull11]
  !sub rax,[L1hinten]
  !mov [L1vorn],rax
  ;If LNull22 < L2hinten
    ;L2hinten = LNull22
  ;EndIf
  !mov rax,[L2hinten]
  !cmp rax,[LNull22]
  !jbe @f  
    !mov rax,[LNull22]
    !mov [L2hinten],rax
  !@@:
  ;L2vorn = LNull22 - L2hinten
  !mov rax,[LNull22]
  !sub rax,[L2hinten]
  !mov [L2vorn],rax

  ;Register-, Status- und Control- Speicher-Reservierung
  BufferFX = AllocateMemory(512 + 64)
  If BufferFX
    BufferFXA = BufferFX + (64 - (BufferFX & $3F))    ;64-er Alignment, klotzen, nicht kleckern!
   Else 
    sAusgabe + #LFCR$ + "Fehler: Nicht genügend Speicher für BufferFX!"
 ProcedureReturn sAusgabe 
  EndIf
  PUSH BufferFXA
  !pop [BufferFXA]

  !mov rax,[BufferFXA]
  !fxsave64 [rax]            ;sichert erstmal komplett FPU und XMM (0-15). Lass ich mal so, obwohl FNSAVE auch reichen würde. Dann aber so keine Register-Sicherung
  !mov [rax+464],rdi         ;callee-save registers sichern, also auf RBX und RBP verzichten
  !mov [rax+472],rsi
  !mov [rax+480],r12
  !mov [rax+488],r13
  !mov [rax+496],r14
  !mov [rax+504],r15

  ;FFT-Speicher-Reservierung; 10=80Bit 
  BufferSin = AllocateMemory((N * 10) + 128)          ;#PB_Memory_NoClear besser nicht
  If BufferSin
    BufferSinA = BufferSin + (64 - (BufferSin & $3F)) ;64-er Alignment, klotzen, nicht kleckern!
   Else 
    sAusgabe + #LFCR$ + "Fehler: Nicht genügend Speicher für BufferSin!"
 ProcedureReturn sAusgabe
  EndIf
  PUSH BufferSinA            ;Protected nach FAsm
  !pop [BufferSinA]

  BufferCos = AllocateMemory((N * 10) + 128)
  If BufferCos
    BufferCosA = BufferCos + (64 - (BufferCos & $3F)) ;64-er Alignment
   Else 
    sAusgabe + #LFCR$ + "Fehler: Nicht genügend Speicher für BufferCos!"
 ProcedureReturn sAusgabe
  EndIf
  PUSH BufferCosA
  !pop [BufferCosA]

  BufferFaktor1 = AllocateMemory((2 * N * 10) + 128)
  If BufferFaktor1
    BufferFaktor1A = BufferFaktor1 + (64 - (BufferFaktor1 & $3F))    ;64-er Alignment
   Else 
    sAusgabe + #LFCR$ + "Fehler: Nicht genügend Speicher für BufferFaktor1!"
 ProcedureReturn sAusgabe
  EndIf
  PUSH BufferFaktor1A
  !pop [BufferFaktor1A]

  BufferFaktor2 = AllocateMemory((2 * N * 10) + 128)
  If BufferFaktor2
    BufferFaktor2A = BufferFaktor2 + (64 - (BufferFaktor2 & $3F))    ;64-er Alignment
   Else 
    sAusgabe + #LFCR$ + "Fehler: Nicht genügend Speicher für BufferFaktor2!"
 ProcedureReturn sAusgabe
  EndIf
  PUSH BufferFaktor2A
  !pop [BufferFaktor2A]

  ;======================================================================================
  TA_Gesamt = ElapsedMilliseconds()
  TA_Faktoren = ElapsedMilliseconds()
  ;Faktoren bearbeiten
  ;Zuerst die Ziffern gemäß reverse Verteilung neu setzen, dabei gleich in Long konvertieren
  !mov r12d,55555555h
  !mov r13d,33333333h
  !mov r14d,0F0F0F0Fh
  !mov r15d,00FF00FFh
  !mov rcx,32
  !sub rcx,[Exponent]

  ;For i = (N / 4) - 1 To 0 Step -1
  !xor r11,r11               ;FakPos
  !mov r10,[N]
  !shr r10,2
  !mov rdx,r10               ;j = N >> 2
  !dec r10                   ;hier nötig

  !.RevPosLoop:              ;local label
    !mov r8,r11              ;müssen für reverse Bits dann 32-Bit-Register sein!
    !shr r8d,1
    !and r8d,r12d
    !mov edi,r11d
    !and edi,r12d
    !shl edi,1
    !or r8d,edi

    !mov edi,r8d             ;sichern
    !shr r8d,2
    !and r8d,r13d
    !and edi,r13d
    !shl edi,2
    !or r8d,edi

    !mov edi,r8d             ;sichern
    !shr r8d,4
    !and r8d,r14d
    !and edi,r14d
    !shl edi,4
    !or r8d,edi

    !mov edi,r8d             ;sichern
    !shr r8d,8
    !and r8d,r15d
    !and edi,r15d
    !shl edi,8
    !or r8d,edi

    !mov edi,r8d             ;sichern
    !shr r8d,16
    !shl edi,16
    !or r8d,edi              ;r8=RevPos

    !shr r8d,cl              ;cl=32-Exponent
    !shl r8d,1               ;lieber doch so

    ;If L1hinten <= 0
      ;WL = 0
     ;Else
      ;WL = PeekB(PF1 + i + j - LNull1) & $F     ;statt -30h; sicherer bei anderer Quelle als String
      ;L1hinten - 1
    ;EndIf 
    ;PokeL(BufferSinA + RevPos, WL)
    !xor r9,r9               ;WL=r9d
    !cmp [L1hinten],0
    !jbe @f
      !mov rax,[PF1]
      !add rax,r10
      !sub rax,[LNull1]
      !movzx r9d,byte[rax+rdx]    ;rdx=j
      !and r9d,0fh
      !dec [L1hinten]
    !@@:
    !mov rax,[BufferSinA]
    !mov dword[rax+r8],r9d

    ;If L2hinten <= 0
      ;WL = 0
     ;Else
      ;WL = PeekB(PF2 + i + j - LNull2) & $F
      ;L2hinten - 1
    ;EndIf
    ;PokeL(BufferCosA + RevPos, WL)
    !xor r9,r9               ;WL=r9d
    !cmp [L2hinten],0
    !jbe @f
      !mov rax,[PF2]
      !add rax,r10
      !sub rax,[LNull2]
      !movzx r9d,byte[rax+rdx]    ;rdx=j 
      !and r9d,0fh
      !dec [L2hinten]
    !@@:
    !mov rax,[BufferCosA]
    !mov dword[rax+r8],r9d

    ;RevPos + 4
    !add r8,4

    ;If L1vorn <= 0
      ;WL = 0
     ;Else
      ;WL = PeekB(PF1 + i - LNull1) & $F
      ;L1vorn - 1
    ;EndIf
    ;PokeL(BufferSinA + RevPos, WL)
    !xor r9,r9               ;WL=r9d
    !cmp [L1vorn],0
    !jbe @f
      !mov rax,[PF1]
      !add rax,r10
      !sub rax,[LNull1]
      !movzx r9d,byte[rax]
      !and r9d,0fh
      !dec [L1vorn]
    !@@:
    !mov rax,[BufferSinA]
    !mov dword[rax+r8],r9d

    ;If L2vorn <= 0
      ;WL = 0
     ;Else
      ;WL = PeekB(PF2 + i - LNull2) & $F
      ;L2vorn - 1
    ;EndIf
    ;PokeL(BufferCosA + RevPos, WL)
    !xor r9,r9               ;WL=r9d
    !cmp [L2vorn],0
    !jbe @f
      !mov rax,[PF2]
      !add rax,r10
      !sub rax,[LNull2]
      !movzx r9d,byte[rax]
      !and r9d,0fh
      !dec [L2vorn]
    !@@:
    !mov rax,[BufferCosA]
    !mov dword[rax+r8],r9d
    ;FakPos + 1
    !inc r11
    !dec r10                 ;i
  ;Next
  !jns .RevPosLoop  

  ;--------------------------------------------------------
  ;die Integer-DWords in 80-Bit-Floats konvertieren und abspeichern
  ;j = 0
  !fninit                              ;FPU wurde gesichert mit fxsave64

  !xor r13,r13               ;i
  !xor r14,r14               ;j
  !mov r15,[N]
  !shl r15,1
  ;For i = 0 To 2 * N Step 4            ;Step 4 = Long
  !@@:
    ;WL = PeekL(BufferSinA + i)
    ;WR = WL
    ;PokeD(BufferFaktor1A + j, WR)
    ;PokeD(BufferFaktor1A + j + 16, WR)
    ;P1 = BufferSinA + i
    !mov rcx,[BufferSinA]
    !add rcx,r13             ;P1=rcx
    ;P2 = BufferFaktor1A + j
    !mov rdx,[BufferFaktor1A]
    !add rdx,r14             ;P2=rdx
    !fild dword[rcx]
    !fld st0                 ;st0 und st1
    !fstp tword[rdx]         ;Real-Wert
    !fstp tword[rdx+20]      ;gedoppelter Real-Wert. Die dazwischen liegenden Imaginär-Werte sind Null (von AllocateMemory, deshalb dort nicht #PB_Memory_NoClear)
    ;WL = PeekL(BufferCosA + i)
    ;WR = WL
    ;PokeD(BufferFaktor2A + j, WR)
    ;PokeD(BufferFaktor2A + j + 16, WR)
    ;P1 = BufferCosA + i
    !mov rcx,[BufferCosA]
    !add rcx,r13             ;P1=rcx
    ;P2 = BufferFaktor2A + j
    !mov rdx,[BufferFaktor2A]
    !add rdx,r14             ;P2=rdx
    !fild dword[rcx]
    !fld st0                 ;st0 und st1
    !fstp tword[rdx]
    !fstp tword[rdx+20]
    ;j + 40
    !add r14,40

    !add r13,4
    !sub r15,4
  ;Next
  !jnz @b
  TE_Faktoren = ElapsedMilliseconds() - TA_Faktoren

  ;========================================================
  ;Winkel erst hier, weil BufferSinA und BufferCosA oben "missbraucht" werden
  TA_Winkel = ElapsedMilliseconds()
  ;Rad1 = Pi / LenF
  !fldpi
  !fild qword[LenF]
  !fdivp st1,st0
  !fstp tword[Rad1]

  ;SinCos; da dies lahm ist und sowieso beide Werte benötigt werden, wird nur bis Pi/4 (45°) ermittelt und der Rest (bis 180°) nur umkopiert
  !xor r14,r14                         ;P00 
  !mov r8,[BufferSinA]
  !mov r9,[BufferCosA]
  !mov r10,10                          ;80Bit=10 Bytes
  !lea r11,[Sin]
  !lea r12,[Cos]
  !mov rcx,[N]
  !shr rcx,3
  ;For k = 0 To (N >> 3)
  !@@:
    ;Rad = Rad1 * k                     ;nicht aufaddieren! Zu ungenau    
    !fld tword[Rad1]
    !push r14
    !fild qword[rsp]                   ;geht nur mit Mem
    !pop r14
    !fmulp st1,st0
    ;Si = Sin(Rad)
    ;Co = Cos(Rad)
    !fsincos
    !fstp tword[Cos]                   ;nebenbei, FST kann tword nicht
    !fstp tword[Sin]
    ;PokeD(BufferSinA + k * 8, Si)     ;0-45°
    !mov rax,r14                       ;[v_P00]
    !mul r10
    !mov r13,rax                       ;k * 8
    !add rax,r8
    !add r13,r9
    !fld tword[r11]
    !fstp tword[rax]                   ;Sin 0-45°
    ;PokeD(BufferCosA + k * 8, Co)     ;0-45°
    !fld tword[r12]
    !fstp tword[r13]                   ;Cos 0-45°
    ;PokeD(BufferCosA + ((N >> 2) - k) * 8, Si)  ;45-90°
    !mov rax,[LenF]
    !shr rax,1
    !sub rax,r14                       ;[v_P00]
    !mul r10
    !mov r13,rax
    !add rax,r9
    !add r13,r8
    !fld tword[r11]
    !fstp tword[rax]                   ;Cos 45-90°
    ;PokeD(BufferSinA + ((N >> 2) - k) * 8, Co)  ;45-90°
    !fld tword[r12]
    !fstp tword[r13]                   ;Sin 45-90°
    ;PokeD(BufferCosA + ((N >> 2) + k) * 8, -Si) ;90-135°
    !mov rax,[LenF]
    !shr rax,1
    !add rax,r14                       ;[v_P00]
    !mul r10
    !mov r13,rax
    !add rax,r9
    !add r13,r8
    !fld tword[r11]
    !fchs
    !fstp tword[rax]                   ;Cos 90-135°  
    ;PokeD(BufferSinA + ((N >> 2) + k) * 8, Co)  ;90-135°
    !fld tword[r12]
    !fstp tword[r13]                   ;Sin 90-135°
    ;PokeD(BufferCosA + (LenF - k) * 8, -Co)     ;135-180°
    !mov rax,[LenF]
    !sub rax,r14                       ;[v_P00]
    !mul r10
    !mov r13,rax
    !add rax,r9
    !add r13,r8
    !fld tword[r12]
    !fchs
    !fstp tword[rax]                   ;Cos 135-180°
    ;PokeD(BufferSinA + (LenF - k) * 8, Si) ;135-180°
    !fld tword[r11]
    !fstp tword[r13]                   ;Sin 135-180°

    !inc r14                           ;[v_P00]
    !dec rcx
  ;Next
  !jns @b
  TE_Winkel = ElapsedMilliseconds() - TA_Winkel

  ;========================================================
  ;Faktor1
  TA_FFT = ElapsedMilliseconds()
  ;Pointer0 = 2
  !mov r12,2
  ;Pointer1 = 0
  !xor rdi,rdi
  ;Pointer2 = Pointer0
  !mov rsi,r12
  !mov r13,2
  ;While Pointer2 < N
  !.OuterLoopF1:
    ;While Pointer2 < N
    !.InnerLoopF1:
      !mov r14,r12
      ;For k = 1 To Pointer0
      !@@:
        ;EW = (N / Pointer0 << 1) * Pointer1 * 10
        !mov r8,[N]
        !mov rax,10
        !mul rdi
        !mov rcx,r13
        !shr r8,cl
        !mul r8
        !mov r15,rax              ;EW=r15     

        ;P0 = Pointer2 * 20
        !mov rax,20
        !mul rsi                  ;P0=rax
        ;P1 = BufferFaktor1A + P0
        !mov rcx,[BufferFaktor1A]
        !add rcx,rax              ;P1=rcx
        ;P5 = (Pointer2 - Pointer0) * 20  
        ;P6 = BufferFaktor1A + P5
        !mov rax,rsi
        !sub rax,r12
        !mov r10,20               ;imul will ich nicht
        !mul r10                  ;P5=rax   an dieser Stelle wegen rdx
        !add rax,[BufferFaktor1A] ;=P6
        ;P2 = BufferCosA + EW
        !mov rdx,[BufferCosA]
        !add rdx,r15
        ;P3 = BufferFaktor1A + 10 + P0
        !mov r8,rcx
        !add r8,10 
        ;P4 = BufferSinA + EW
        !mov r9,[BufferSinA]
        !add r9,r15               ;P4=r9
        ;WR = (PeekD(BufferFaktor1A + 2 * Pointer2 * 8) * PeekD(BufferCosA + EW) - PeekD(BufferFaktor1A + 8 + 2*Pointer2 * 8) * PeekD(BufferSinA + EW))
        ;->WR = PeekD(P1) * PeekD(P2) - PeekD(P3) * PeekD(P4)
        !fld tword[rcx]
        !fld tword[rdx]
        !fmulp st1,st0
        !fld tword[r8]
        !fld tword[r9]
        !fmulp st1,st0
        !fsubp st1,st0
        !fstp tword[WR]
        ;WI = (PeekD(BufferFaktor1A + 8 + 2 * Pointer2 * 8) * PeekD(BufferCosA + EW) + PeekD(BufferFaktor1A + 2 * Pointer2 * 8) * PeekD(BufferSinA + EW))
        ;->WI = PeekD(P3) * PeekD(P2) + PeekD(P1) * PeekD(P4)
        !fld tword[r8]
        !fld tword[rdx]
        !fmulp st1,st0
        !fld tword[rcx]
        !fld tword[r9]
        !fmulp st1,st0
        !faddp st1,st0
        !fstp tword[WI]
        ;ZR = PeekD(BufferFaktor1A + 2 * (Pointer2 - Pointer0) * 8)
        ;->ZR = PeekD(P6)
        ;PokeD(BufferFaktor1A + 2 * (Pointer2 - Pointer0) * 8, ZR + WR)
        ;->PokeD(P6, ZR + WR)
        !fld tword[rax]
        !fld st0
        !fstp tword[ZR]
        !fld tword[WR]
        !faddp st1,st0
        !fstp tword[rax]
        ;P7 = P6 + 10             ;P6=rax
        ;ZI = PeekD(BufferFaktor1A + 8 + 2 * (Pointer2 - Pointer0) * 8)
        ;->ZI = PeekD(P7)
        ;PokeD(BufferFaktor1A + 8 + 2 * (Pointer2 - Pointer0) * 8, ZI + WI)
        ;->PokeD(P7, ZI + WI)
        !fld tword[rax+10]
        !fld st0
        !fstp tword[ZI]
        !fld tword[WI]
        !faddp st1,st0
        !fstp tword[rax+10]
        ;PokeD(BufferFaktor1A + 2 * Pointer2 * 8, ZR - WR)
        ;->PokeD(P1, ZR - WR)
        !fld tword[ZR]
        !fld tword[WR]
        !fsubp st1,st0
        !fstp tword[rcx]
        ;PokeD(BufferFaktor1A + 8 + 2 * Pointer2 * 8, ZI - WI)
        ;->PokeD(P3, ZI - WI)
        !fld tword[ZI]
        !fld tword[WI]
        !fsubp st1,st0
        !fstp tword[r8]
        ;Pointer1 + 1
        !inc rdi;[v_Pointer1]
        ;Pointer2 + 1
        !inc rsi;[v_Pointer2]
        !dec r14
      ;Next
      !jnz @b

      ;Pointer1 = 0
      !xor rdi,rdi
     ;Pointer2 + Pointer0
      !add rsi,r12
    ;Wend     While Pointer2 < N
      !cmp rsi,[N]
    !jb .InnerLoopF1

    !inc r13
    ;Pointer0 << 1
    !shl r12,1
    ;Pointer1 = 0
    !xor rdi,rdi
    ;Pointer2 = Pointer0
    !mov rsi,r12
  ;Wend     While Pointer2 < N
    !cmp rsi,[N]
  !jb .OuterLoopF1

  ;========================================================
  ;Faktor2
  ;Pointer0 = 2
  !mov r12,2
  ;Pointer1 = 0
  !xor rdi,rdi
  ;Pointer2 = Pointer0
  !mov rsi,r12
  !mov r13,2
  ;While Pointer2 < N
  !.OuterLoopF2:
    ;While Pointer2 < N
    !.InnerLoopF2:
      ;For k = 1 To Pointer0
        !mov r14,r12
      !@@:
        ;EW = (N / Pointer0 << 1) * Pointer1 * 10
        !mov r8,[N]
        !mov rax,10
        !mul rdi
        !mov rcx,r13
        !shr r8,cl
        !mul r8
        !mov r15,rax              ;EW=r15     

        ;P0 = Pointer2 * 20
        !mov rax,20
        !mul rsi                  ;P0=rax
        ;P1 = BufferFaktor2A + P0
        !mov rcx,[BufferFaktor2A]
        !add rcx,rax   ;P1=rcx
        ;P5 = (Pointer2 - Pointer0) * 20  
        ;P6 = BufferFaktor2A + P5
        !mov rax,rsi
        !sub rax,r12
        !mov r10,20               ;imul will ich nicht
        !mul r10                  ;P5=rax   an dieser Stelle wegen rdx
        !add rax,[BufferFaktor2A] ;=P6
        ;P2 = BufferCosA + EW
        !mov rdx,[BufferCosA]
        !add rdx,r15;[v_EW];P2=rdx
        ;P3 = BufferFaktor2A + 10 + P0
        !mov r8,rcx
        !add r8,10 
        ;P4 = BufferSinA + EW
        !mov r9,[BufferSinA]
        !add r9,r15   ;[v_EW];P4=r9
        ;WR = (PeekD(BufferFaktor2A + 2 * Pointer2 * 8) * PeekD(BufferCosA + EW) - PeekD(BufferFaktor2A + 8 + 2*Pointer2 * 8) * PeekD(BufferSinA + EW))
        ;->WR = PeekD(P1) * PeekD(P2) - PeekD(P3) * PeekD(P4)
        !fld tword[rcx]
        !fld tword[rdx]
        !fmulp st1,st0
        !fld tword[r8]
        !fld tword[r9]
        !fmulp st1,st0
        !fsubp st1,st0
        !fstp tword[WR]
        ;WI = (PeekD(BufferFaktor2A + 8 + 2 * Pointer2 * 8) * PeekD(BufferCosA + EW) + PeekD(BufferFaktor2A + 2 * Pointer2 * 8) * PeekD(BufferSinA + EW))
        ;->WI = PeekD(P3) * PeekD(P2) + PeekD(P1) * PeekD(P4)
        !fld tword[r8]
        !fld tword[rdx]
        !fmulp st1,st0
        !fld tword[rcx]
        !fld tword[r9]
        !fmulp st1,st0
        !faddp st1,st0
        !fstp tword[WI]
        ;ZR = PeekD(BufferFaktor2A + 2 * (Pointer2 - Pointer0) * 8)
        ;->ZR = PeekD(P6)
        ;PokeD(BufferFaktor2A + 2 * (Pointer2 - Pointer0) * 8, ZR + WR)
        ;->PokeD(P6, ZR + WR)
        !fld tword[rax]
        !fld st0
        !fstp tword[ZR]
        !fld tword[WR]
        !faddp st1,st0
        !fstp tword[rax]
        ;P7 = P6 + 10             ;P6=rax
        ;ZI = PeekD(BufferFaktor2A + 8 + 2 * (Pointer2 - Pointer0) * 8)
        ;->ZI = PeekD(P7)
        ;PokeD(BufferFaktor2A + 8 + 2 * (Pointer2 - Pointer0) * 8, ZI + WI)
        ;->PokeD(P7, ZI + WI)
        !fld tword[rax+10]
        !fld st0
        !fstp tword[ZI]
        !fld tword[WI]
        !faddp st1,st0
        !fstp tword[rax+10]
        ;PokeD(BufferFaktor2A + 2 * Pointer2 * 8, ZR - WR)
        ;->PokeD(P1, ZR - WR)
        !fld tword[ZR]
        !fld tword[WR]
        !fsubp st1,st0
        !fstp tword[rcx]
        ;PokeD(BufferFaktor2A + 8 + 2 * Pointer2 * 8, ZI - WI)
        ;->PokeD(P3, ZI - WI)
        !fld tword[ZI]
        !fld tword[WI]
        !fsubp st1,st0
        !fstp tword[r8]
        ;Pointer1 + 1
        !inc rdi
        ;Pointer2 + 1
        !inc rsi
        !dec r14
      ;Next
      !jnz @b

      ;Pointer1 = 0
      !xor rdi,rdi
      ;Pointer2 + Pointer0
      !add rsi,r12
    ;Wend     While Pointer2 < N
      !cmp rsi,[N]
    !jb .InnerLoopF2

    !inc r13
    ;Pointer0 << 1
    !shl r12,1
    ;Pointer1 = 0
    !xor rdi,rdi
    ;Pointer2 = Pointer0
    !mov rsi,r12
  ;Wend     While Pointer2 < N
    !cmp rsi,[N]
  !jb .OuterLoopF2
  TE_FFT = ElapsedMilliseconds() - TA_FFT 

  ;========================================================
  ;Und jetzt die beiden transformierten Vektoren miteinander multiplizieren
  TA_Multi = ElapsedMilliseconds()
  !mov r10,[N]               ;ohne -1, dafür unten jnz und nicht jns
  !xor rax,rax               ;rax=P0
  ;For k = 0 To N - 1
  !@@:
    ;P1 = BufferFaktor1A + P0
    !mov rcx,[BufferFaktor1A]
    !add rcx,rax             ;rcx=P1
    ;P2 = BufferFaktor2A + P0
    !mov rdx,[BufferFaktor2A]
    !add rdx,rax
    ;P3 = BufferFaktor1A + 10 + P0
    !mov r8,rcx
    !add r8,10               ;r8=P3 
    ;P4 = BufferFaktor2A + 10 + P0
    !mov r9,rdx
    !add r9,10               ;r9=P4 
    ;WR = (PeekD(BufferFaktor1A + 2 * k * 8) * PeekD(BufferFaktor2A + 2 * k * 8)) - (PeekD(BufferFaktor1A + 8 + 2 * k * 8) * PeekD(BufferFaktor2A + 8 + 2 * k * 8))    ;hier kann nicht direkt in BufferFaktor1RA geschrieben werden
    ;->WR = PeekD(P1) * PeekD(P2) - PeekD(P3) * PeekD(P4)
    !fld tword[rcx]
    !fld tword[rdx]
    !fmulp st1,st0
    !fld tword[r8]
    !fld tword[r9]
    !fmulp st1,st0
    !fsubp st1,st0
    !fstp tword[WR]

    ;PokeD(BufferFaktor1A + 8 + 2 * k * 8, (PeekD(BufferFaktor1A + 2 * k * 8) * PeekD(BufferFaktor2A + 8 + 2 * k * 8)) + (PeekD(BufferFaktor1A + 8 + 2 * k * 8) * PeekD(BufferFaktor2A + 2 * k * 8)))
    ;->PokeD(P3, PeekD(P1) * PeekD(P4) + PeekD(P3) * PeekD(P2))
    !fld tword[rcx]
    !fld tword[r9]
    !fmulp st1,st0
    !fld tword[r8]
    !fld tword[rdx]
    !fmulp st1,st0
    !faddp st1,st0
    !fstp tword[r8]

    ;PokeD(BufferFaktor1A + 2 * k * 8, WR)  ;weil vorher BufferFaktor1RA noch benötigt wird
    ;->PokeD(P1, WR)
    !fld tword[WR]
    !fstp tword[rcx]

    !add rax,20              ;rax=P0=k*20
    !dec r10
  ;Next
  !jnz @b
  TE_Multi = ElapsedMilliseconds() - TA_Multi 

  ;========================================================
  ;Die Produkt-Werte wieder revers "verteilen"
  TA_Rueck = ElapsedMilliseconds()
  !mov r12d,55555555h
  !mov r13d,33333333h
  !mov r14d,0F0F0F0Fh
  !mov r15d,00FF00FFh
  !mov rcx,32
  !sub rcx,[Exponent]

  !xor r11,r11               ;Zähler und FakPos zugleich
  !mov r10,[N]
  !shr r10,1
  ;For i = 0 To (N / 2) - 1 
  !@@:                       ;erste und letzte Ziffern könnten übergangen werden
    !mov r8,r11              ;müssen für reverse Bits 32-Bit-Register sein!!!
    !shr r8d,1
    !and r8d,r12d
    !mov edi,r11d
    !and edi,r12d
    !shl edi,1
    !or r8d,edi

    !mov edi,r8d             ;sichern
    !shr r8d,2
    !and r8d,r13d
    !and edi,r13d
    !shl edi,2
    !or r8d,edi

    !mov edi,r8d             ;sichern
    !shr r8d,4
    !and r8d,r14d
    !and edi,r14d
    !shl edi,4
    !or r8d,edi

    !mov edi,r8d             ;sichern
    !shr r8d,8
    !and r8d,r15d
    !and edi,r15d
    !shl edi,8
    !or r8d,edi

    !mov edi,r8d             ;sichern
    !shr r8d,16
    !shl edi,16
    !or r8d,edi              ;r8=RevPos

    !shr r8d,cl              ;cl=32-Exponent

    ;PokeD(BufferFaktor2A + 2 * 8 * RevPos, PeekD(BufferFaktor1A + 2 * i * 8))
    ;P1 = BufferFaktor2A + 2 * 10 * RevPos
    !mov rax,20
    !mul r8                  ;[v_RevPos]
    !add rax,[BufferFaktor2A]
    !mov rsi,rax
    ;P2 = BufferFaktor1A + 2 * i * 10
    !mov rax,20
    !mul r11
    !add rax,[BufferFaktor1A]
    !mov r9,rax

    !fld tword[rax]
    !fstp tword[rsi]
    ;PokeD(BufferFaktor2A + 8 + 2 * 8 * RevPos, PeekD(BufferFaktor1A + 8 + 2 * i * 8))
    ;P1 = P1 + 10
    ;P2 = P2 + 10
    !fld tword[rax+10]
    !fstp tword[rsi+10]
    ;PokeD(BufferFaktor2A + 2 * 8 * (RevPos + 1), PeekD(BufferFaktor1A + (2 * i * 8) + N * 8))
    ;P1 = BufferFaktor2A + 2 * 10 * (RevPos + 1)
    !inc r8
    !mov rax,20
    !mul r8
    !add rax,[BufferFaktor2A]
    !mov rsi,rax
    ;P2 = BufferFaktor1A + (2 * i * 10) + N * 10
    ;P2 + N * 10
    !mov rax,10
    !mul [N]
    !add rax,r9
    !fld tword[rax]
    !fstp tword[rsi]
    ;PokeD(BufferFaktor2A + 8 + 2 * 8 * (RevPos + 1), PeekD(BufferFaktor1A + 8 + (2 * i * 8) + N * 8))
    ;P1 = P1 + 10
    ;P2 = P2 + 10
    !fld tword[rax+10]
    !fstp tword[rsi+10]

    !inc r11                 ;ist auch FakPos
    !dec r10
  ;Next
  !jnz @b
  FreeMemory(BufferFaktor1)

  ;========================================================
  ;Anfangswerte verteilen für FFT
  !xor rax,rax               ;rax=P0
  !mov r10,[N]               ;ohne -2, dafür unten jnz und nicht jns
  ;For i = 0 To N - 2 Step 2
  !@@:
    ;P00 = (i + 1) * 20
    ;->P00 = P0 + 20
    ;P1 = BufferFaktor2A + P0
    !mov rcx,[BufferFaktor2A]
    !add rcx,rax             ;rcx=P1
    ;P2 = BufferFaktor2A + P00
    !mov rdx,rcx
    !add rdx,20              ;rdx=P2
    ;P3 = BufferFaktor2A + 10 + P0
    !mov r8,rcx
    !add r8,10               ;r8=P3
    ;P4 = BufferFaktor2A + 10 + P00
    !mov r9,r8
    !add r9,20               ;r9=P4

    ;WR = PeekD(BufferFaktor2A + 2 * i * 8) + PeekD(BufferFaktor2A + 2 * (i + 1) * 8)
    ;->WR = PeekD(P1) + PeekD(P2)
    !fld tword[rcx]
    !fld tword[rdx]
    !faddp st1,st0
    !fstp tword[WR]
    ;WI = PeekD(BufferFaktor2A + 8 + 2 * i * 8) + PeekD(BufferFaktor2A + 8 + 2 * (i + 1) * 8)
    ;->WI = PeekD(P3) + PeekD(P4)
    !fld tword[r8]
    !fld tword[r9]
    !faddp st1,st0
    !fstp tword[WI]
    ;ZR = PeekD(BufferFaktor2A + 2 * i * 8) - PeekD(BufferFaktor2A + 2 * (i + 1) * 8)
    ;->ZR = PeekD(P1) - PeekD(P2)
    !fld tword[rcx]
    !fld tword[rdx]
    !fsubp st1,st0
    ;PokeD(BufferFaktor2A + 2 * (i + 1) * 8, ZR)
    !fstp tword[rdx]
    ;ZI = PeekD(BufferFaktor2A + 8 + 2 * i * 8) - PeekD(BufferFaktor2A + 8 + 2 * (i + 1) * 8)
    ;->ZI = PeekD(P3) - PeekD(P4)
    !fld tword[r8]
    !fld tword[r9]
    !fsubp st1,st0
    ;PokeD(BufferFaktor2A + 8 + 2 * (i + 1) * 8, ZI)
    !fstp tword[r9]
    ;PokeD(BufferFaktor2A + 2 * i * 8, WR)
    !fld tword[WR]
    !fstp tword[rcx]
    ;PokeD(BufferFaktor2A + 8 + 2 * i * 8, WI)
    !fld tword[WI]
    !fstp tword[r8]

    !add rax,40              ;P0 = i * 20
    !sub r10,2
  ;Next
  !jnz @b

  ;==========================================================
  ;nochmal FFT, jetzt aber invers (mit neg.Sinus)
  ;Pointer0 = 2
  !mov r12,2
  ;Pointer1 = 0
  !xor rdi,rdi
  ;Pointer2 = Pointer0
  !mov rsi,r12

  !mov r13,2
  ;While Pointer2 < N               ;da erstmal wahr hier fußgesteuert
  !.OuterLoopP:
    ;While Pointer2 < N
    !.InnerLoopP:
      !mov r14,r12
      ;For k = 1 To Pointer0
      !@@:
        ;EW = (N / Pointer0 << 1) * Pointer1 * 10
        !mov r8,[N]
        !mov rax,10
        !mul rdi;[v_Pointer1]
        !mov rcx,r13
        !shr r8,cl
        !mul r8
        !mov r15,rax              ;EW=r15     

        ;P0 = Pointer2 * 20
        !mov rax,20
        !mul rsi                  ;P0=rax
        ;P1 = BufferFaktor2A + P0
        !mov rcx,[BufferFaktor2A]
        !add rcx,rax              ;P1=rcx
        ;P5 = (Pointer2 - Pointer0) * 20  
        ;P6 = BufferFaktor2A + P5
        !mov rax,rsi 
        !sub rax,r12 
        !mov r10,20               ;imul will ich nicht
        !mul r10                  ;P5=rax   an dieser Stelle wegen rdx
        !add rax,[BufferFaktor2A] ;=P6
        ;P2 = BufferCosA + EW
        !mov rdx,[BufferCosA]
        !add rdx,r15
        ;P3 = BufferFaktor2A + 10 + P0
        !mov r8,rcx
        !add r8,10 
        ;P4 = BufferSinA + EW
        !mov r9,[BufferSinA]
        !add r9,r15               ;P4=r9
        ;WR = (PeekD(BufferFaktor2A + 2 * Pointer2 * 8) * PeekD(BufferCosA + EW) + PeekD(BufferFaktor2A + 8 + 2*Pointer2 * 8) * PeekD(BufferSinA + EW))
        ;->WR = PeekD(P1) * PeekD(P2) - PeekD(P3) * PeekD(P4)
        !fld tword[rcx]
        !fld tword[rdx]
        !fmulp st1,st0
        !fld tword[r8]
        !fld tword[r9]
        !fmulp st1,st0
        !faddp st1,st0
        !fstp tword[WR]
        ;WI = (PeekD(BufferFaktor2A + 8 + 2 * Pointer2 * 8) * PeekD(BufferCosA + EW) - PeekD(BufferFaktor2A + 2 * Pointer2 * 8) * PeekD(BufferSinA + EW))
        ;->WI = PeekD(P3) * PeekD(P2) + PeekD(P1) * PeekD(P4)
        !fld tword[r8]
        !fld tword[rdx]
        !fmulp st1,st0
        !fld tword[rcx]
        !fld tword[r9]
        !fmulp st1,st0
        !fsubp st1,st0
        !fstp tword[WI]
        ;ZR = PeekD(BufferFaktor2A + 2 * (Pointer2 - Pointer0) * 8)
        ;->ZR = PeekD(P6)
        ;PokeD(BufferFaktor2A + 2 * (Pointer2 - Pointer0) * 8, ZR + WR)
        ;->PokeD(P6, ZR + WR)
        !fld tword[rax]
        !fld st0
        !fstp tword[ZR]
        !fld tword[WR]
        !faddp st1,st0
        !fstp tword[rax]
        ;P7 = P6 + 10             ;P6=rax
        ;ZI = PeekD(BufferFaktor2A + 8 + 2 * (Pointer2 - Pointer0) * 8)
        ;->ZI = PeekD(P7)
        ;PokeD(BufferFaktor2A + 8 + 2 * (Pointer2 - Pointer0) * 8, ZI + WI)
        ;->PokeD(P7, ZI + WI)
        !fld tword[rax+10]
        !fld st0
        !fstp tword[ZI]
        !fld tword[WI]
        !faddp st1,st0
        !fstp tword[rax+10]
        ;PokeD(BufferFaktor2A + 2 * Pointer2 * 8, ZR - WR)
        ;->PokeD(P1, ZR - WR)
        !fld tword[ZR]
        !fld tword[WR]
        !fsubp st1,st0
        !fstp tword[rcx]
        ;PokeD(BufferFaktor2A + 8 + 2 * Pointer2 * 8, ZI - WI)
        ;->PokeD(P3, ZI - WI)
        !fld tword[ZI]
        !fld tword[WI]
        !fsubp st1,st0
        !fstp tword[r8]
        ;Pointer1 + 1
        !inc rdi
        ;Pointer2 + 1
        !inc rsi
        !dec r14
      ;Next
      !jnz @b

      ;Pointer1 = 0
      !xor rdi,rdi
      ;Pointer2 + Pointer0
      !add rsi,r12
    ;Wend     While Pointer2 < N
      !cmp rsi,[N]
    !jb .InnerLoopP

    !inc r13
    ;Pointer0 << 1
    !shl r12,1
    ;Pointer1 = 0
    !xor rdi,rdi
    ;Pointer2 = Pointer0
    !mov rsi,r12

  ;Wend     While Pointer2 < N
    !cmp rsi,[N]
  !jb .OuterLoopP
  TE_Rueck = ElapsedMilliseconds() - TA_Rueck

  ;==========================================================
  FreeMemory(BufferCos)
  FreeMemory(BufferSin)

  TA_Ergebnis = ElapsedMilliseconds()
  ;die Doubles von FFT in Integer konvertieren und Koeffizienten aufaddieren mit Zehner-Versatz
  BufferProdukt = AllocateMemory((N * 8) + 128)
  If BufferProdukt
    BufferProduktA = BufferProdukt + (64 - (BufferProdukt & $3F)); + 64 ;64-er Alignment, +64 zur Vermeidung Unterlauf (unsauber, aber schneller!)
   Else 
    sAusgabe + #LFCR$ + "Fehler: Nicht genügend Speicher für BufferProdukt!"
 ProcedureReturn sAusgabe
  EndIf
  PUSH BufferProduktA
  !pop [BufferProduktA]

  BufferDez = AllocateMemory(32 + 128)   ;BufferDez ist Protected
  If BufferDez
    BufferDezA = BufferDez + (64 - (BufferDez & $3F)) ;64-er Alignment
   Else 
    sAusgabe + #LFCR$ + "Fehler: Nicht genügend Speicher für BufferDez!"
 ProcedureReturn sAusgabe
  EndIf
  PUSH BufferDezA
  !pop [BufferDezA]          ;BufferDezA jetzt für FAsm 

  !mov r15,[N] 
  !xor r14,r14
  !lea rcx,[WR]
  !lea r10,[Dezi]
  !mov r11,[BufferDezA]
  ;Pointer_Produkt = N - 1
  !mov rdi,r15
  !dec rdi
  ;For i = 0 To N - 1         ;vom Produkt alle Werte durch
  !fninit                    ;zur Sicherheit, setzt auch Rundung auf round to nearest
  !.DeziLoop:
    ;Division der Produktwerte durch N und Konvertierung Double in Integer
    ;QuadWert = Round(PeekD(BufferFaktor2A + i * 16), #PB_Round_Nearest)  ;Double zu Quad
    !mov rax,20
    !mul r14
    !mov r9,[BufferFaktor2A]
    !fld tword[r9+rax]       ;Reihenfolge getauscht, FNINIT oben reicht zur Sicherheit
    !fild [N] 
    !fdivp st1,st0           ;P(OP) jetzt nötig
    ;!fabs                    ;wohl nicht nötig  
    !fistp qword[rcx]        ;jetzt als Quad-Integer abspeichern. SSE2 (cvtsd2si) bringt jetzt hier nichts  
    !mov rax,[rcx]           ;rax=QuadWert

    ;Hex2Dez
    ;Ziffer = 0
    !xor r12b,r12b 
    ;Pointer1 = 0
    !xor r8,r8
    ;Pointer2 = 0
    !xor r9,r9
 
    ;While PeekQ(?Dezi + Pointer1) <> 0 oder Zähler wie jetzt
    !mov rsi,11              ;Anzahl Dezi-Werte ohne Null
    !.OuterLoopDez1:
      ;While (QuadWert - PeekQ(?Dezi + Pointer1)) >= 0
      !mov r13,rax
      !.InnerLoopDez1:
        !sub r13,[r10+r8]
        !js .InnerLoopDez1End
        ;Ziffer + 1
        !inc r12b
        ;QuadWert - PeekQ(?Dezi + Pointer1)
        ;!sub rax,[r10+r8]
        !mov rax,r13         ;AL ist unten letzte Ziffer
      ;Wend     (QuadWert - PeekQ(?Dezi + Pointer1)) >= 0
      !jmp .InnerLoopDez1

      !.InnerLoopDez1End:
      ;PokeB(BufferDezA + Pointer2, Ziffer)
      !mov [r11+r9],r12b
      ;Pointer1 + 8
      !add r8,8
      ;Pointer2 + 1
      !inc r9
      ;Ziffer = 0
      !xor r12b,r12b 
      !dec rsi
    ;Wend      PeekQ(?Dezi + Pointer1) <> 0 oder Zähler wie jetzt
    !jnz .OuterLoopDez1

    ;PokeB(BufferDezA + Pointer2, QuadWert1)   ;letzte Dezimal-Ziffer abspeichern
    !mov [r11+r9],al         ;AL ist letzte Ziffer von oben
    ;Pointer_Produkt1 = Pointer_Produkt 
    !mov rsi,rdi
    !mov rax,11              ;11 reicht, evtl. sogar anpassen je nach Exponent
    !mov rdx,[BufferProduktA]
    ;For k = 0 To 15
    !@@:
      ;Ziffer = PeekB(BufferProduktA + Pointer_Produkt1)
      !mov r12b,[rdx+rsi]
      ;Ziffer + PeekB(BufferDezA + Pointer2)
      !add r12b,[r11+r9]
      ;PokeB(BufferProduktA + Pointer_Produkt1, Ziffer)
      !mov [rdx+rsi],r12b
      ;If Ziffer > 9                           ;Übertrag
      !cmp r12b,9              ;Übertrag?
      !jbe .ZifferOK
        ;Ziffer - 10
        !sub r12b,10
        ;PokeB(BufferProduktA + Pointer_Produkt1, Ziffer)
        !mov [rdx+rsi],r12b
        ;Ziffer = PeekB(BufferProduktA + Pointer_Produkt1 - 1) + 1
        !mov r12b,[rdx+rsi-1]
        !inc r12b
        ;PokeB(BufferProduktA + Pointer_Produkt1 - 1, Ziffer)
        !mov [rdx+rsi-1],r12b
      ;EndIf
      !.ZifferOK:
      ;Pointer_Produkt1 - 1
      !dec rsi
      !js .Reicht
      ;Pointer2 - 1
      !dec r9
      !dec rax
    ;Next
    !jnz @b

    !.Reicht:
    ;Pointer_Produkt - 1
    !dec rdi
    !inc r14
    !dec r15
  ;Next
  !jnz .DeziLoop

  ;------------------------
  ;Ergebnis (Produkt) in ASCII-Ziffern abspeichern
  !mov rdx,3030303030303030h
  !mov rcx,[N]
  !mov rax,[BufferProduktA]
  !mov byte[rax+rcx],0       ;Zero-Byte zur Sicherheit für String-Auslesen
  !sub rcx,8
  !@@:
    !or [rax+rcx],rdx
    !sub rcx,8
  !jns @b

  FreeMemory(BufferFaktor2)
  FreeMemory(BufferDez)

  sProdukt = LTrim(PeekS(BufferProduktA), "0")   ;evtl. führende Null(en) weg
  If sProdukt = ""
    sProdukt = "0"
  EndIf

  FreeMemory(BufferProdukt)

  TE_Ergebnis = ElapsedMilliseconds() - TA_Ergebnis

  TE_Gesamt = ElapsedMilliseconds() - TA_Gesamt

  sFaktor = "Faktor-Laenge (max.): " + Str(L1) + " Bytes (=Ziffern)" + #CRLF$
  sExpo = "Exponent 2^ für angepasste Faktoren-Laenge: " + Str(Exponent - 1) + #CRLF$
  sWinkel = "Zeit für Sinus/Cosinus-Berechnung: " + Str(TE_Winkel) + " ms" + #CRLF$
  sFaktoren = "Zeit für Faktoren-Aufbereitung: " + Str(TE_Faktoren) + " ms" + #CRLF$
  sFFT = "Zeit für Fast Fourier Transformation: " + Str(TE_FFT) + " ms" + #CRLF$
  sMulti = "Zeit für Multiplikation: " + Str(TE_Multi) + " ms" + #CRLF$
  sRueck = "Zeit für Rueck-Transformation: " + Str(TE_Rueck) + " ms" + #CRLF$
  sErgebnis = "Zeit für Ergebnis-Aufbereitung als String: " + Str(TE_Ergebnis) + " ms" + #CRLF$
  sGesamt = "Zeit gesamt: " + Str(TE_Gesamt) + " ms"
  sAusgabe + sFaktor + sExpo + sWinkel + sFaktoren + sFFT + sMulti + sRueck + sErgebnis + sGesamt + "XX" + sProdukt    ;"XX" dient als Trenner zwischen Infos und Produkt

  !mov rax,[BufferFXA]
  !fxrstor64 [rax]
  !mov rdi,[rax+464]         ;callee-save registers zurücksichern, RBX und RBP wurden nicht gesichert (und von mir nicht benutzt)
  !mov rsi,[rax+472]
  !mov r12,[rax+480]
  !mov r13,[rax+488]
  !mov r14,[rax+496]
  !mov r15,[rax+504]
  FreeMemory(BufferFX)       ;Register-mäßig interessant :-)

  ;DisableASM                 ;geschenkt
 ProcedureReturn sAusgabe

  ;--------------------------------------------------------
  DataSection
    !Rad1            dt ?    ;10-Byte-Variablen
    !                dp ?    ;6 Bytes Dummy wegen Alignment
    !Sin             dt ?
    !                dp ?
    !Cos             dt ?
    !                dp ?
    !WR              dt ?
    !                dp ?
    !WI              dt ?
    !                dp ?
    !ZR              dt ?
    !                dp ?
    !ZI              dt ?
    !                dp ?

    !BufferCosA      dq ?
    !BufferDezA      dq ?
    !BufferFaktor1A  dq ?
    !BufferFaktor2A  dq ?
    !BufferFXA       dq ?
    !BufferProduktA  dq ?
    !BufferSinA      dq ?
    !Exponent        dq ?
    !L1              dq ?
    !L2              dq ?
    !LenF            dq ?
    !LNull1          dq ?
    !LNull11         dq ?
    !L1hinten        dq ?
    !L1vorn          dq ?
    !LNull2          dq ?
    !LNull22         dq ?
    !L2hinten        dq ?
    !L2vorn          dq ?
    !N               dq ?
    !PF1             dq ?
    !PF2             dq ?

    !Dezi            dq     100000000000    ;einhundert Milliarden reichen bis 1GB Faktoren-Länge!
    !                dq      10000000000
    !                dq       1000000000
    !                dq        100000000
    !                dq         10000000
    !                dq          1000000
    !                dq           100000
    !                dq            10000
    !                dq             1000
    !                dq              100
    !                dq               10
    !                dq                0    ;bleibt zur Sicherheit
  EndDataSection
 
EndProcedure


;ein Auswahl-Menü schenke ich mir hier. Einen der Tests auswählen, Goto oder nicht Goto, das ist hier die Frage... ;-)
;Quellen-Auswahl: (Faktoren erstellen, kann jeder selbst generieren)
;Goto NixLongString
Faktor1$ = "6582018229284824168619876730229402019930943462534319453394436096"  ;64-Bytes
Faktor2$ = "5486721320215405789450123459785103469751240504554109782315460125"

FL = 6                  ;Faktoren-Länge: FL=6=4KB(2ms), wird im MessageRequester noch mit Produkt angezeigt, Werte von 0 bis 21 (je nach RAM-Ausstattung) möglich, 22=Out of Memory bei mir  
For i = 1 To FL         ;für Test. Faktoren-Länge: FL=14=1MB(1s); 15=2MB(3s); 16=4MB(8s); 17=8MB(17s); 18=16MB(37s); 19=32MB(80s); 20=64MB(173s); 21=128MB(369s); 22=Out of Memory (32GB!)
  Faktor1$ + Faktor1$   ;um auf "Länge" zu kommen :-)
  Faktor2$ + Faktor2$  
Next
String = 1
NixLongString:
;----------------------------------------------------------
Goto NixLongStringKrumm
;oder
Faktor1$ = "6730229402019930943462534319453394436096"
Faktor2$ = "9751240504554109782315460125"
For i = 1 To 10
  Faktor1$ + Faktor1$
  Faktor2$ + Faktor2$  
Next
String = 1
NixLongStringKrumm:
;----------------------------------------------------------
Goto NixShortString
;oder
Faktor1$ = "10"
Faktor2$ = "15"
String = 1
NixShortString:
;----------------------------------------------------------
Goto NixMem
;oder, wenn die Faktoren irgendwo im Speicher stehen, z.B.:
LFMem1 = 6                        ;z.B.6 Ziffern
PFMem1 = AllocateMemory(LFMem1)
For i = 0 To LFMem1 - 1
  PokeB(PFMem1 + i, i + 1)        ;sollen die 6 Ziffern sein. Bei höher 9 natürlich andere Wertezuweisung (Bereich 0-9)
Next
LFMem2 = 5                        ;z.B.5 Ziffern
PFMem2 = AllocateMemory(LFMem2)
For i = 0 To LFMem2 - 1
  PokeB(PFMem2 + i, i + 1)        ;sollen die 5 Ziffern sein. Bei höher 9 natürlich andere Wertezuweisung (Bereich 0-9)
Next
;String = 0
Nixmem:
;----------------------------------------------------------

;die FFT_Mul-Procedure aufrufen. Um die Unicode-Problematik zu umgehen empfehle ich natürlich die Memory-Variante, was in einer Anwendung sowieso so sein sollte  
If String
  Result$ = FFT_Mul(Len(Faktor1$), Len(Faktor2$), @Faktor1$, @Faktor2$)
 Else
  Result$ = FFT_Mul(LFMem1, LFMem2, PFMem1, PFMem2)
EndIf

;die Rückgabe der Procedure auswerten
Pos_XX = FindString(Result$, "XX", 1, #PB_String_CaseSensitive)
If Pos_XX                              ;Pos_XX ist Null bei Fehler  
  L_Result = Len(Result$)
  L_Produkt = L_Result - Pos_XX - 2    ;2=XX
  Produkt$ = Mid(Result$, Pos_XX + 2)
  Info$ = Mid(Result$, 1, L_Result - L_Produkt - 3)

  SetClipboardText(Produkt$)           ;oder eigene Methode

  If L_Result < 10000                  ;damit der MessageRequester nicht schlapp macht. Direkt angezeigt wird bis Faktor-Länge 4KB
    MessageRequester("Helles Integer-Multiplikation-FFT-Test", Info$ + #LFCR$ + "Produkt:" + #LFCR$ + Produkt$ + #LFCR$ + #LFCR$ + "Das Produkt ist auch in der Zwischenablage zu finden!")
   Else
    MessageRequester("Helles Integer-Multiplikation-FFT-Test", Info$ + #LFCR$ + #LFCR$ + "Das Produkt ist für die Anzeige zu groß und deshalb nur in der Zwischenablage zu finden!")
  EndIf

 Else
  MessageRequester("Helles Integer-Multiplikation-FFT-Test", Result$)     ;irgendeine Fehlermeldung wurde generiert
EndIf
Viel Spaß!
Helle

Edit 09.01.2016: Kleine Optimierung im Bereich DeziLoop
Zuletzt geändert von Helle am 09.01.2016 11:17, insgesamt 1-mal geändert.
C h r i s :)
Beiträge: 56
Registriert: 12.11.2015 17:41

Re: Schnelle Integer-Multiplikation mit FFT

Beitrag von C h r i s :) »

so viel asm codes :shock: :o
ich versteh nur bahnhof :freak:
ist das überhaupt pflegbar?
Benutzeravatar
NicTheQuick
Ein Admin
Beiträge: 8679
Registriert: 29.08.2004 20:20
Computerausstattung: Ryzen 7 5800X, 32 GB DDR4-3200
Ubuntu 22.04.3 LTS
GeForce RTX 3080 Ti
Wohnort: Saarbrücken
Kontaktdaten:

Re: Schnelle Integer-Multiplikation mit FFT

Beitrag von NicTheQuick »

Mein Rechner ist deutlich langsamer als deiner, aber bis FL = 21 kam ich auf dem Laptop auch. :mrgreen:
Ich hab 304 Sekunden für FL = 20 gebraucht. FL = 21 rechnet noch.

Ansonsten top! :allright:

Schon mal überlegt z.B. unter Linux die fertigen FFT-Libraries zu nutzen? Ich weiß gar nicht, ob ich meinen Import dafür schon hier im Forum gepostet habe. Das muss ich mal prüfen.
Bild
Benutzeravatar
7x7
Beiträge: 591
Registriert: 14.08.2007 15:41
Computerausstattung: ganz toll
Wohnort: Lelbach

Re: Schnelle Integer-Multiplikation mit FFT

Beitrag von 7x7 »

ASM...nur fliegen ist schöner! :mrgreen:
- alles was ich hier im Forum sage/schreibe ist lediglich meine Meinung und keine Tatsachenbehauptung
- unkommentierter Quellcode = unqualifizierter Müll
Benutzeravatar
Sicro
Beiträge: 955
Registriert: 11.08.2005 19:08
Kontaktdaten:

Re: Schnelle Integer-Multiplikation mit FFT

Beitrag von Sicro »

Wird in CodeArchiv unter Maths/Fast_integer_multiplication_with_FFT.pbi aufgenommen.
Bild
Warum OpenSource eine Lizenz haben sollte :: PB-CodeArchiv-Rebirth :: Pleasant-Dark (Syntax-Farbschema) :: RegEx-Engine (kompiliert RegExes zu NFA/DFA)
Manjaro Xfce x64 (Hauptsystem) :: Windows 10 Home (VirtualBox) :: Neueste PureBasic-Version
Antworten