Wollte eben wie angekündigt meinen obigen Ursprungs-Code ersetzen, aber es kam die Meldung mit dem Limit von 60.000 Zeichen. Deshalb hier der neue optimierte Code in 2 Teilen.
Ich hatte mir selbst das Ziel gesetzt unter 2 Minuten für M49 zu bleiben, aber bei 125 Sekunden war Schluß. Einige Versuche endeten im Stack-Nirwana und toll angedachte Ideen waren Sackgassen. So bleibt es also erstmal bei dieser 1 Thread-Variante. Wollte auch noch den Code optisch aufarbeiten (Einrückungen), aber naja, auch so ist der Code prima übersichtlich und verständlich
.
Code: Alles auswählen
;- Schnelle Integer-Multiplikation mit FFT für große Zahlen, hier als Beispiel die Berechnung von 2-er Potenzen
;- Bis auf die max.Länge keine Einschränkungen bei den Faktoren
;- 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
;- Zur besseren Orientierung im Code ist der Großteil des PB-Codes noch auskommentiert enthalten. Kann natürlich entfernt werden
;- Windows7/64, PureBasic 5.41 LTS (x64)
;- "Helle" Klaus Helbing, 01.02.2016
Procedure.s FFT_Pot2(Expo, Korr)
EnableASM
;EnableExplicit ;war für Übersicht mal nötig
Protected.q BufferCos, BufferCosA, BufferDez, BufferDezA, BufferErg, BufferErgA, BufferFX, BufferFXA, BufferFaktor1, BufferFaktor1A, BufferFaktor2, BufferFaktor2A
Protected.q BufferProd, BufferProdA, BufferProdukt, BufferProduktA, BufferSin, BufferSinA, N
Protected.q Is1Bit, L1, L2, Limit, Stellen, SZ, Ergebnis, Potenz
Protected.s sAusgabe
!mov rcx,[p.v_Expo] ;zur Sicherheit so, sonst rcx direkt (als erstes natürlich)
!mov [Potenz],rcx
!mov rdx,[p.v_Korr]
!mov [KorrWert],rdx
;die Ziffern-Stellenzahl des Ergebnisses ermitteln
;Stellen = (Potenz * Log10(2)) + 32 ;math. +1, hier mehr für Alignment und Freiraum nach hinten (ZeroByte(s))
;Stellen = Stellen + (16 - (Stellen & $1F))
!fldlg2
!fild qword[Potenz]
!fmulp st1,st0
!fistp qword[Stellen]
!mov rax,[Stellen]
!add rax,32
!mov rdx,rax
!and rdx,1fh
!add rax,16
!sub rax,rdx
!mov [Stellen],rax
!push [Stellen]
POP Stellen
!bsr rax,[Potenz] ;höchstes gesetztes Bit
!jnz @f
!mov rax,1 ;Exponent=0
!add rax,[KorrWert]
!push rax
POP Ergebnis
sAusgabe = Str(Ergebnis) ;2^0=1
ProcedureReturn sAusgabe
!@@:
!or rax,rax
!jnz @f
!mov rax,2 ;Exponent=1
!add rax,[KorrWert]
!push rax
POP Ergebnis
sAusgabe = Str(Ergebnis) ;2^1=2
ProcedureReturn sAusgabe
!@@:
;Test auf 2-er Potenz
!bsf rdx,[Potenz]
!cmp rax,rdx
!jne @f
!mov [Is1Bit],1
!@@:
!push rax
POP Limit
!push [Is1Bit]
POP Is1Bit
BufferProd = AllocateMemory(Stellen + 64)
If BufferProd
BufferProdA = BufferProd + (64 - (BufferProd & $3F)) ;64-er Alignment
Else
sAusgabe = #LFCR$ + "Fehler: Nicht genügend Speicher für BufferProd!"
ProcedureReturn sAusgabe
EndIf
PUSH BufferProdA
!pop [BufferProdA]
BufferErg = AllocateMemory(Stellen + 64)
If BufferErg
BufferErgA = BufferErg + (64 - (BufferErg & $3F)) ;64-er Alignment
Else
sAusgabe = #LFCR$ + "Fehler: Nicht genügend Speicher für BufferErg!"
ProcedureReturn sAusgabe
EndIf
PUSH BufferErgA
!pop [BufferErgA]
;Startwerte
!mov r8,[BufferProdA]
!mov byte[r8],2
!mov qword[ProdLen],1
!mov rax,[BufferErgA]
!mov byte[rax],1
!mov [ErgLen],1
!shr [Potenz],1 ;Bit0 testen
!jnc @f ;Bit0 ist 0
!mov byte[rax],2
!@@:
BufferDez = AllocateMemory(32 + 64) ;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
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
;==================================================================================================
;==================================================================================================
For SZ = 1 To Limit
;L1=PeekQ(BufferLenA); + (SZ - 1) * 8)
!mov rax,[ProdLen]
!mov [L1],rax ;L1 wegen Kompatibilität erstmal lassen
;L1 = LFak1
!mov rax,[L1]
;LNull11 = L1 ;echte Länge, auch hier Null-Test möglich
!mov [LNull11],rax ;echte Länge, auch hier Null-Test möglich
;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,[L1] ;niedrigstes gesetztes Bit
!mov [Exponent],rcx
!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,4 ;Produkt-Länge muss mindestens 16 sein, sonst an anderen Stellen Klimmzüge
!sub rcx,[Exponent]
!js @f ;Exponent ist > 4 (Hauptfall)
!jz @f ;Exponent ist = 4
!shl [LenF],cl
!add [Exponent],rcx
!@@:
!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
;L1hinten = LenF >> 1 ;LenF ist 2-er Potenz
!mov rax,[LenF]
!shr rax,1 ;LenF ist 2-er Potenz
!mov [L1hinten],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
;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]
;======================================================================================
;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,[BufferProdA]
!add rax,r10
!sub rax,[LNull1]
!movzx r9d,byte[rax+rdx] ;rdx=j
!dec [L1hinten]
!@@:
!mov rax,[BufferSinA]
!mov dword[rax+r8],r9d
!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,[BufferProdA]
!add rax,r10
!sub rax,[LNull1]
!movzx r9d,byte[rax]
!dec [L1vorn]
!@@:
!mov rax,[BufferSinA]
!mov dword[rax+r8],r9d
;PokeL(BufferCosA + RevPos, WL)
!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)
;j + 40
!add r14,40
!add r13,4
!sub r15,4
;Next
!jnz @b
;========================================================
;Winkel erst hier, weil BufferSinA (und BufferCosA) oben "missbraucht" werden
;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
;========================================================
;Faktor1
;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
;========================================================
;Und jetzt die beiden (gleichen) transformierten Vektoren miteinander multiplizieren
!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
;P3 = BufferFaktor1A + 10 + P0
!mov r8,rcx
!add r8,10 ;r8=P3
;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]
!fmul st0,st0
!fld tword[r8]
!fmul st0,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[r8]
!fmulp st1,st0
!fadd st0,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
;========================================================
;Die Produkt-Werte wieder revers "verteilen"
!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
;==========================================================
FreeMemory(BufferCos)
FreeMemory(BufferSin)
;die Doubles von FFT in Integer konvertieren und Koeffizienten aufaddieren mit Zehner-Versatz
BufferProdukt = AllocateMemory(N + 128)
If BufferProdukt
BufferProduktA = BufferProdukt + (64 - (BufferProdukt & $3F)) ;64-er Alignment
Else
sAusgabe + #LFCR$ + "Fehler: Nicht genügend Speicher für BufferProdukt!"
ProcedureReturn sAusgabe
EndIf
PUSH BufferProduktA
!pop [BufferProduktA]
!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 rdx,[BufferProduktA]
!mov rax,11 ;11 reicht, evtl. sogar anpassen je nach Exponent
;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
;------------------------
;jetzt suchen wir mal von vorn in BufferProduktA die erste von Null verschiedene Ziffer
!pxor xmm1,xmm1 ;xmm1 null
!mov rcx,rdx ;rdx=BufferProduktA
!@@:
!movdqu xmm0,[rdx] ;16 Bytes einlesen
!add rdx,16
!pcmpeqb xmm0,xmm1 ;alle 16 Bytes auf Null (xmm1) testen, Byte=255 wenn Null und Byte=0 wenn nicht
!pmovmskb eax,xmm0 ;das jeweilige MSB des Bytes nach eax (ax) kopieren
!cmp eax,0FFFFh ;noch alles Nullen?
!je @b ;ja
!sub rdx,17 ;zurück, um erste Ziffer ungleich Null zu finden
!@@:
!inc rdx
!shr eax,1
!jc @b
!sub rdx,rcx ;rdx ist Pointer auf 1.Ziffer <> 0 in BufferProduktA
!mov rcx,[N]
!sub rcx,rdx ;Anzahl Bytes
;PokeQ(BufferLenA + (SZ * 8), N - X)
!mov [ProdLen],rcx
;PokeS(BufferPotA + (SZ * Stellen), sAusgabe)
!mov r9,[BufferProduktA]
!add r9,rdx ;Beginn Bytes
;!mov r10,[BufferPotA]
!mov rax,[BufferProdA]
!@@:
!mov r8b,[r9]
!mov [rax],r8b
!inc r9
!inc rax
!dec rcx
!jz .WarKlein
!test rcx,7
!jnz @b
!@@:
!mov r8,[r9]
!mov [rax],r8
!add r9,8
!add rax,8
!sub rcx,8
!jnz @b
!.WarKlein:
FreeMemory(BufferFaktor2)
FreeMemory(BufferProdukt)
;======================================================================================
;======================================================================================
!push [Potenz]
POP Potenz
!shr [Potenz],1
If Is1Bit = 0 And (Potenz & 1) ;um ein Sprung-Adressen-Problem zu umgehen
!mov rax,[ProdLen] ;ProdLen ist immer mindestens genauso groß wie das bisherige ErgLen
!mov [L1],rax
!mov [LNull11],rax ;echte Länge
!mov [LenF],rax ;LenF erstmal setzen, könnte ja schon 2-er Potenz sein
;L1=PeekQ(BufferLenA + ((Limit - 1) * 8)) ;feste Stelle
!mov rdx,[ErgLen]
!mov [L2],rdx ;L1 und L2 bleiben erstmal!
!mov [LNull22],rdx
!bsr rcx,[L1] ;höchstes gesetztes Bit
!bsf rdx,[L1] ;niedrigstes gesetztes Bit
!mov [Exponent],rcx
!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,4 ;Produkt-Länge muss mindestens 16 sein, sonst an anderen Stellen Klimmzüge
!sub rcx,[Exponent]
!js @f ;Exponent ist > 4 (Hauptfall)
!jz @f ;Exponent ist = 4
!shl [LenF],cl
!add [Exponent],rcx
!@@:
!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
;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]
;======================================================================================
;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
!.RevPosLoop1: ;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,[BufferProdA]
!add rax,r10
!sub rax,[LNull1]
!movzx r9d,byte[rax+rdx] ;rdx=j
!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,[BufferErgA]
!add rax,r10
!sub rax,[LNull2]
!movzx r9d,byte[rax+rdx] ;rdx=j
!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,[BufferProdA]
!add rax,r10
!sub rax,[LNull1]
!movzx r9d,byte[rax]
!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,[BufferErgA]
!add rax,r10
!sub rax,[LNull2]
!movzx r9d,byte[rax]
!dec [L2vorn]
!@@:
!mov rax,[BufferCosA]
!mov dword[rax+r8],r9d
;FakPos + 1
!inc r11
!dec r10 ;i
;Next
!jns .RevPosLoop1
;--------------------------------------------------------
;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
;========================================================
;Winkel erst hier, weil BufferSinA und BufferCosA oben "missbraucht" werden
;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
;========================================================
;Faktor1
;Pointer0 = 2
!mov r12,2
;Pointer1 = 0
!xor rdi,rdi
;Pointer2 = Pointer0
!mov rsi,r12
!mov r13,2
;While Pointer2 < N
!.OuterLoopF11:
;While Pointer2 < N
!.InnerLoopF11:
!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 .InnerLoopF11
!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 .OuterLoopF11
;========================================================
;Faktor2
;Pointer0 = 2
!mov r12,2
;Pointer1 = 0
!xor rdi,rdi
;Pointer2 = Pointer0
!mov rsi,r12
!mov r13,2
;While Pointer2 < N
!.OuterLoopF21:
;While Pointer2 < N
!.InnerLoopF21:
;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 .InnerLoopF21
!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 .OuterLoopF21
;========================================================
;Und jetzt die beiden transformierten Vektoren miteinander multiplizieren
!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
;========================================================
;Die Produkt-Werte wieder revers "verteilen"
!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)
;========================================================
! Habe heute vormittag im Büro (freitags geht das schon mal 
Hat 24 Mio Stellen und ist damit größer als die bis jetzt gefunde Primzahl. 
Was soll ich jetzt machen?
!