Neue größte bekannte Primzahl
Verfasst: 20.01.2016 23:22
Auf Heise war ja heute zu lesen, dass die bisher größte bekannte Primzahl (M49?) als solche bestätigt wurde. Wer die "Zahl" (Ziffern) haben will, kann unten aufgeführten Code benutzen. Es wird auf die Schnelle die Procedure aus meinem FFT-Thread verwendet und ich habe überlegt, ob ich wirklich einen neuen Thread aufmachen soll (darf). Da ich aber die Procedure hier noch für diesen speziellen Zweck anpassen will und diese Proceduren dann doch sehr unterschiedlich sein werden, habe ich diesen Thread aufgemacht. Dies als Entschuldigung an die Mods!
Also wie eben schon angedeutet, das ganze hier ist absolut unoptimiert, verschwendet Speicher usw.; werde ich noch ändern. Die Heise-Meldung hat mich doch sehr angestachelt und ich wollte schnell sein...
Die Berechnung dauert bei mir (i7-6700K@4.7GHz) ca.145s.
Viel Spaß!
Helle
Also wie eben schon angedeutet, das ganze hier ist absolut unoptimiert, verschwendet Speicher usw.; werde ich noch ändern. Die Heise-Meldung hat mich doch sehr angestachelt und ich wollte schnell sein...
Die Berechnung dauert bei mir (i7-6700K@4.7GHz) ca.145s.
Code: Alles auswählen
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
ZeitA.q=ElapsedMilliseconds()
Global Is1Bit.q
Global Limit.q
Potenz.q = 74207281 ;M49=(2^74207281)-1 Potenz = 2 ist Minimum!
Loggi.d = Log10(2)
Stellen.q = (Potenz * Loggi) + 1 ;+1 ist nicht immer nötig, aber natürlich so lassen
!bsr rax,[v_Potenz] ;höchstes gesetztes Bit
!jnz @f
;Fehler, Null und Ende; 2 ist Minimum!
!@@:
!or rax,rax
!jnz @f
;Fehler, 1 und Ende
!@@:
;Test auf 2-er Potenz
!bsf rdx,[v_Potenz]
!cmp rax,rdx
!jne @f
!mov [v_Is1Bit],1
!@@:
!add rax,2 ;für Bit0 und Zwischenergebnis
!mov [v_Limit],rax
BufferPot = AllocateMemory((Stellen * Limit) + 128)
If BufferPot
BufferPotA = BufferPot + (64 - (BufferPot & $3F)) ;64-er Align
Else
MessageRequester("Hallo", "Fehler!")
EndIf
BufferLen = AllocateMemory(8 * Limit + 128)
If BufferLen
BufferLenA = BufferLen + (64 - (BufferLen & $3F)) ;64-er Align
Else
MessageRequester("Hallo", "Fehler!")
EndIf
PokeS(BufferPotA, "1")
PokeS(BufferPotA + (Limit - 1) * Stellen, "1")
PokeQ(BufferLenA, 1)
PokeQ(BufferLenA + (Limit - 1) * 8, 1)
!shr [v_Potenz],1 ;Bit0 testen
!jnc @f ;Bit0 ist 0
PokeS(BufferPotA, "2")
PokeS(BufferPotA + (Limit - 1) * Stellen, "2")
!@@:
Result1$ = "2"
For i = 1 To Limit - 2
Result$ = FFT_Mul(Len(Result1$), Len(Result1$), @Result1$, @Result1$)
Pos_XX = FindString(Result$, "XX", 1, #PB_String_CaseSensitive)
Result1$ = Mid(Result$, Pos_XX + 2)
PokeS(BufferPotA + (i * Stellen), Result1$)
PokeQ(BufferLenA + (i * 8), Len(Result1$))
Next
If Is1Bit = 0
For i = 1 To Limit - 2
!shr [v_Potenz],1
!jnc @f
Result$ = FFT_Mul(PeekQ(BufferLenA + ((Limit - 1) * 8)), PeekQ(BufferLenA + (i * 8)), BufferPotA + ((Limit - 1) * Stellen), BufferPotA + (i * Stellen))
Pos_XX = FindString(Result$, "XX", 1, #PB_String_CaseSensitive)
Result1$ = Mid(Result$, Pos_XX + 2)
PokeS(BufferPotA + ((Limit - 1) * Stellen), Result1$)
PokeQ(BufferLenA + ((Limit - 1) * 8), Len(Result1$))
!@@:
Next
EndIf
LM = Len(Result1$)
LetzteZiffern$ = Mid(Result1$, LM - 1) ;mit Kanonen auf Spatzen...
NeueZiffern = Val(LetzteZiffern$) - 1
M$ = InsertString(Left(Result1$, LM - 2), Str(NeueZiffern), LM - 1)
ZeitE.q = ElapsedMilliseconds() - ZeitA
SetClipboardText(M$)
MessageRequester("Helles Mersenne-M49-Berechnung mittels FFT", Str(ZeitE) + #LFCR$ + "Das Ergebnis ist für die Anzeige zu groß und deshalb nur in der Zwischenablage zu finden!")
Helle
!