Meine 5GHz-CPU spuckt die 23.249.425 Ziffern von Mersenne M50 übrigens nach 44s aus.
Hier der (noch verbesserungsbedürftige) Code:
Code: Alles auswählen
;- Schnelle Unsigned Integer-Multiplication mit FFT
;- see https://en.wikipedia.org/wiki/Butterfly_diagram (radix-2 decimation-in-time FFT algorithm)
;- Windows 64-Bit, getestet mit PureBasic 5.46 LTS (x64) / PureBasic 5.61 (x64)
;- Ascii und Unicode
;- Benötigt CPU mit AVX2-Unterstützung
;- "Helle" Klaus Helbing, 05.01.2018: Mersenne M50 hinzugefügt
;- Um in den Genuss der Info-Anzeige zu kommen, sollten Programme wie der Taskmanager nicht gestartet sein!
;- Version vom 11.01.2018
Global.q Asc_Uni, Bits_M, CalcPot, FakA, GadgetZeile, Potenz, ZeitA, ZeitE
Global.s Basis$, Result$
Declare.s FFT_Mul(Para1, Para2)
Declare Calc()
Declare RangeCharSSE42(Para1, Para2)
Declare FakInfo()
;Test auf AVX2
!xor ecx,ecx
!mov eax,7 ;Test auf 7 geschenkt
!cpuid
!test ebx,20h ;Bit5 für AVX2
!jnz @f ;Yeah!
MessageRequester("Schade!", "Diese CPU unterstützt kein AVX2!" + #LFCR$ + "Ende")
End
!@@:
FensterTitel$ = "Helles FFT-Calculator AVX2-Version"
OpenWindow(0, 0, 0, 800, 600, FensterTitel$, #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
CompilerIf #PB_Compiler_Unicode
Asc_Uni = 2 : #Asc_Uni = 2 ;mit Konstante für Step usw.
CompilerElse
Asc_Uni = 1 : #Asc_Uni = 1
CompilerEndIf
FontHigh = Int(9.0 / (GetDeviceCaps_(GetDC_(WindowID(0)), #LOGPIXELSY) / 96.0)) ;GetDeviceCaps kann bei der Bildschirmanpassung nützlich sein
;LoadFont(0, "Arial", FontHigh)
LoadFont(0, "Trebuchet MS Fett", FontHigh) ;in Windows 7 verfügbar; evtl.testen
Repeat
TextGadget(1, 20, 20, 300, 15, "Was soll berechnet werden:")
OptionGadget(2, 20, 45, 780, 15, "Multiplikation zweier positiver Integer-Faktoren")
OptionGadget(3, 20, 65, 780, 15, "Berechnung einer positiven ganzzahligen Potenz")
OptionGadget(4, 20, 85, 780, 15, "Berechnung einer Mersenne-Prim-Zahl (M25...M50)")
OptionGadget(5, 20, 105, 780, 15, "Fakultät")
OptionGadget(6, 20, 125, 780, 15, "Neuner-Test")
ButtonGadget(7, 300, 560, 200, 30, "W e i t e r")
SetGadgetState(2, 1) : FFT = 2 ;Voreinstellung Multiplikation, falls nichts ausgewählt wird
For i = 1 To 7
SetGadgetFont(i, FontID(0))
Next
Repeat
Event = WaitWindowEvent()
If Event = #PB_Event_CloseWindow
CloseWindow(0)
End
EndIf
If Event = #PB_Event_Gadget
Select EventGadget()
Case 7
Break
Default
FFT = EventGadget()
EndSelect
EndIf
ForEver
For i = 1 To 7
FreeGadget(i)
Next
Select FFT
;=========================================================================
Case 2 ;Multiplikation
Info0$ = "Multiplikation "
SetWindowTitle(0, Info0$)
CharRange$ = Chr(1) + Chr(47) + Chr(58) + Chr(255) ;für Suche nach Chars, die keine Dezimal-Ziffer sind. Chr(0) geht nicht!
TextGadget(1, 50, 0, 300, 20, "Faktor1 eintragen (z.B. einen String mit Strg+V)")
TextGadget(2, 450, 0, 300, 20, "Faktor2 eintragen (z.B. einen String mit Strg+V)")
ButtonGadget(3, 325, 560, 150, 30, "Berechnung starten")
ButtonGadget(4, 125, 560, 150, 30, "Faktor1 überprüfen")
ButtonGadget(5, 525, 560, 150, 30, "Faktor2 überprüfen")
For i = 1 To 5
SetGadgetFont(i, FontID(0))
Next
EditorGadget(6, 10, 20, 385, 530)
EditorGadget(7, 405, 20, 385, 530)
SetActiveGadget(6)
Repeat
Event = WaitWindowEvent()
If Event = #PB_Event_CloseWindow
CloseWindow(0)
End
EndIf
If Event = #PB_Event_Gadget
Gadget = EventGadget()
Select Gadget
Case 3
Break
Case 4, 5
If Gadget = 4
Faktor$ = GetGadgetText(6)
Else
Faktor$ = GetGadgetText(7)
EndIf
RangeChar = RangeCharSSE42(@Faktor$, @CharRange$)
If RangeChar
Ascii = Asc(Mid(Faktor$, RangeChar, 1)) : Char$ = Mid(Faktor$, RangeChar, 1)
MessageRequester("Faktor" + Str(Gadget - 3) + ": Fehler!", "Digit an Stelle " + Str(RangeChar) + " ist keine Dezimal-Ziffer! (ASCII=" + Ascii + "), Char=" + Char$, #MB_ICONERROR)
Else
MessageRequester("Faktor" + Str(Gadget - 3) + ": O.K.!", "Alle Digits sind Dezimal-Ziffern!")
EndIf
EndSelect
EndIf
ForEver
Faktor1$ = GetGadgetText(6)
Faktor2$ = GetGadgetText(7)
For i = 1 To 7
FreeGadget(i)
Next
ZeitA = ElapsedMilliseconds()
Result$ = FFT_Mul(@Faktor1$, @Faktor2$)
ZeitE = ElapsedMilliseconds() - ZeitA
;=========================================================================
Case 3 ;Potenzen
Info0$ = "Berechnung von X^Y "
SetWindowTitle(0, Info0$)
TextGadget(1, 50, 10, 300, 20, "Basis X eintragen für X^Y")
TextGadget(2, 400, 10, 300, 20, "Exponent Y eintragen für X^Y")
StringGadget(3, 50, 50, 300, 20, "2", #PB_String_Numeric)
StringGadget(4, 400, 50, 300, 20, "1", #PB_String_Numeric)
ButtonGadget(5, 300, 560, 200, 30, "Berechnung starten")
SetActiveGadget(2)
For i = 1 To 5 ;das StringGadget einfach mal drin lassen
SetGadgetFont(i, FontID(0))
Next
Repeat
Repeat
Event = WaitWindowEvent()
If Event = #PB_Event_CloseWindow
CloseWindow(0)
End
EndIf
Until EventGadget() = 5
Potenz = Val(GetGadgetText(4)) ;Exponent
If (Len(GetGadgetText(4)) > 17) Or (Potenz > 1E16) ;willkürlich gewählte Werte, sichere Seite
MessageRequester("Abbruch!", "Exponent zu groß! Bitte neue Eingabe!")
Else
Break
EndIf
ForEver
Basis$ = GetGadgetText(3)
Basis = 0
If Len(Basis$) < 2
Basis = Len(Basis$)
EndIf
For i = 1 To 5
FreeGadget(i)
Next
Info0$ = Basis$ + "^" + Str(Potenz) + ": "
SetWindowTitle(0, Info0$)
If Len(Basis$) > 0 Or Basis > 0
If Potenz > 0
Result$ = Basis$
Calc()
Else
Result$ = "1"
EndIf
Else
Result$ = "0" ;0^Y
EndIf
;=========================================================================
Case 4 ;Mersenne
Info0$ = "Berechnung von Mersenne "
SetWindowTitle(0, Info0$ + "Mxx")
;Werte für Mersenne M25-M49, die Strings sind nur zur Info
M25.q = 21701 : M25$ = "(2^21701)-1" : Z25$ = "44867916611904333479...57410828353511882751"
M26.q = 23209 : M26$ = "(2^23209)-1" : Z26$ = "40287411577898877818...36743355523779264511"
M27.q = 44497 : M27$ = "(2^44497)-1" : Z27$ = "85450982430363380319...44867686961011228671"
M28.q = 86243 : M28$ = "(2^86243)-1" : Z28$ = "53692799550275632152...99857021709433438207"
M29.q = 110503 : M29$ = "(2^110503)-1" : Z29$ = "52192831334175505976...69951621083465515007"
M30.q = 132049 : M30$ = "(2^132049)-1" : Z30$ = "51274027626932072381...52138578455730061311"
M31.q = 216091 : M31$ = "(2^216091)-1" : Z31$ = "74609310306466134368...91336204103815528447"
M32.q = 756839 : M32$ = "(2^756839)-1" : Z32$ = "17413590682008709732...02603793328544677887"
M33.q = 859433 : M33$ = "(2^859433)-1" : Z33$ = "12949812560420764966...02414267243500142591"
M34.q = 1257787 : M34$ = "(2^1257787)-1" : Z34$ = "41224577362142867472...31257188976089366527"
M35.q = 1398269 : M35$ = "(2^1398269)-1" : Z35$ = "81471756441257307514...85532025868451315711"
M36.q = 2976221 : M36$ = "(2^2976221)-1" : Z36$ = "62334007624857864988...76506256743729201151"
M37.q = 3021377 : M37$ = "(2^3021377)-1" : Z37$ = "12741168303009336743...25422631973024694271"
M38.q = 6972593 : M38$ = "(2^6972593)-1" : Z38$ = "43707574412708137883...35366526142924193791"
M39.q = 13466917 : M39$ = "(2^13466917)-1" : Z39$ = "92494773800670132224...30073855470256259071"
M40.q = 20996011 : M40$ = "(2^20996011)-1" : Z40$ = "12597689545033010502...94714065762855682047"
M41.q = 24036583 : M41$ = "(2^24036583)-1" : Z41$ = "29941042940415717208...67436921882733969407"
M42.q = 25964951 : M42$ = "(2^25964951)-1" : Z42$ = "12216463006127794810...98933257280577077247"
M43.q = 30402457 : M43$ = "(2^30402457)-1" : Z43$ = "31541647561884608093...11134297411652943871"
M44.q = 32582657 : M44$ = "(2^32582657)-1" : Z44$ = "12457502601536945540...11752880154053967871"
M45.q = 37156667 : M45$ = "(2^37156667)-1" : Z45$ = "20225440689097733553...21340265022308220927"
M46.q = 42643801 : M46$ = "(2^42643801)-1" : Z46$ = "16987351645274162247...84101954765562314751"
M47.q = 43112609 : M47$ = "(2^43112609)-1" : Z47$ = "31647026933025592314...80022181166697152511"
M48.q = 57885161 : M48$ = "(2^57885161)-1" : Z48$ = "58188726623224644217...46141988071724285951"
M49.q = 74207281 : M49$ = "(2^74207281)-1" : Z49$ = "30037641808460618205...87010073391086436351"
M50.q = 77232917 : M50$ = "(2^77232917)-1" : Z50$ = "46733318335923109998...82730618069762179071"
TextGadget(1, 20, 0, 300, 15, "Wähle eine Mersenne-Zahl:")
OptionGadget(2, 20, 20, 780, 15, "Mersenne M25, " + M25$ + ", 6533 Dezimal-Ziffern, " + Z25$)
OptionGadget(3, 20, 40, 780, 15, "Mersenne M26, " + M26$ + ", 6987 Dezimal-Ziffern, " + Z26$)
OptionGadget(4, 20, 60, 780, 15, "Mersenne M27, " + M27$ + ", 13395 Dezimal-Ziffern, " + Z27$)
OptionGadget(5, 20, 80, 780, 15, "Mersenne M28, " + M28$ + ", 25962 Dezimal-Ziffern, " + Z28$)
OptionGadget(6, 20, 100, 780, 15, "Mersenne M29, " + M29$ + ", 33265 Dezimal-Ziffern, " + Z29$)
OptionGadget(7, 20, 120, 780, 15, "Mersenne M30, " + M30$ + ", 39751 Dezimal-Ziffern, " + Z30$)
OptionGadget(8, 20, 140, 780, 15, "Mersenne M31, " + M31$ + ", 65050 Dezimal-Ziffern, " + Z31$)
OptionGadget(9, 20, 160, 780, 15, "Mersenne M32, " + M32$ + ", 227832 Dezimal-Ziffern, " + Z32$)
OptionGadget(10, 20, 180, 780, 15, "Mersenne M33, " + M33$ + ", 258716 Dezimal-Ziffern, " + Z33$)
OptionGadget(11, 20, 200, 780, 15, "Mersenne M34, " + M34$ + ", 378632 Dezimal-Ziffern, " + Z34$)
OptionGadget(12, 20, 220, 780, 15, "Mersenne M35, " + M35$ + ", 420921 Dezimal-Ziffern, " + Z35$)
OptionGadget(13, 20, 240, 780, 15, "Mersenne M36, " + M36$ + ", 895832 Dezimal-Ziffern, " + Z36$)
OptionGadget(14, 20, 260, 780, 15, "Mersenne M37, " + M37$ + ", 909526 Dezimal-Ziffern, " + Z37$)
OptionGadget(15, 20, 280, 780, 15, "Mersenne M38, " + M38$ + ", 2098960 Dezimal-Ziffern, " + Z38$)
OptionGadget(16, 20, 300, 780, 15, "Mersenne M39, " + M39$ + ", 4053946 Dezimal-Ziffern, " + Z39$)
OptionGadget(17, 20, 320, 780, 15, "Mersenne M40, " + M40$ + ", 6320430 Dezimal-Ziffern, " + Z40$)
OptionGadget(18, 20, 340, 780, 15, "Mersenne M41, " + M41$ + ", 7235733 Dezimal-Ziffern, " + Z41$)
OptionGadget(19, 20, 360, 780, 15, "Mersenne M42, " + M42$ + ", 7816230 Dezimal-Ziffern, " + Z42$)
OptionGadget(20, 20, 380, 780, 15, "Mersenne M43, " + M43$ + ", 9152052 Dezimal-Ziffern, " + Z43$)
OptionGadget(21, 20, 400, 780, 15, "Mersenne M44, " + M44$ + ", 9808358 Dezimal-Ziffern, " + Z44$)
OptionGadget(22, 20, 420, 780, 15, "Mersenne M45, " + M45$ + ", 11185272 Dezimal-Ziffern, " + Z45$)
OptionGadget(23, 20, 440, 780, 15, "Mersenne M46, " + M46$ + ", 12837064 Dezimal-Ziffern, " + Z46$)
OptionGadget(24, 20, 460, 780, 15, "Mersenne M47, " + M47$ + ", 12978189 Dezimal-Ziffern, " + Z47$)
OptionGadget(25, 20, 480, 780, 15, "Mersenne M48, " + M48$ + ", 17425170 Dezimal-Ziffern, " + Z48$)
OptionGadget(26, 20, 500, 780, 15, "Mersenne M49, " + M49$ + ", 22338618 Dezimal-Ziffern, " + Z49$)
OptionGadget(27, 20, 520, 780, 15, "Mersenne M50, " + M50$ + ", 23249425 Dezimal-Ziffern, " + Z50$)
SetGadgetState(2, 1)
ButtonGadget(28, 300, 560, 200, 30, "Berechnung starten")
For i = 1 To 28
SetGadgetFont(i, FontID(0))
Next
Basis$ = "2"
Result$ = "2"
Potenz = M25 : M$ = "M25=" + M25$ : Z$ = Z25$ ;Voreinstellung Mersenne, falls nichts ausgewählt wird
Repeat
Event = WaitWindowEvent()
If Event = #PB_Event_CloseWindow
CloseWindow(0)
End
EndIf
If Event = #PB_Event_Gadget
Select EventGadget()
Case 2
Potenz = M25 : M$ = "M25=" + M25$ : Z$ = Z25$
Case 3
Potenz = M26 : M$ = "M26=" + M26$ : Z$ = Z26$
Case 4
Potenz = M27 : M$ = "M27=" + M27$ : Z$ = Z27$
Case 5
Potenz = M28 : M$ = "M28=" + M28$ : Z$ = Z28$
Case 6
Potenz = M29 : M$ = "M29=" + M29$ : Z$ = Z29$
Case 7
Potenz = M30 : M$ = "M30=" + M30$ : Z$ = Z30$
Case 8
Potenz = M31 : M$ = "M31=" + M31$ : Z$ = Z31$
Case 9
Potenz = M32 : M$ = "M32=" + M32$ : Z$ = Z32$
Case 10
Potenz = M33 : M$ = "M33=" + M33$ : Z$ = Z33$
Case 11
Potenz = M34 : M$ = "M34=" + M34$ : Z$ = Z34$
Case 12
Potenz = M35 : M$ = "M35=" + M35$ : Z$ = Z35$
Case 13
Potenz = M36 : M$ = "M36=" + M36$ : Z$ = Z36$
Case 14
Potenz = M37 : M$ = "M37=" + M37$ : Z$ = Z37$
Case 15
Potenz = M38 : M$ = "M38=" + M38$ : Z$ = Z38$
Case 16
Potenz = M39 : M$ = "M39=" + M39$ : Z$ = Z39$
Case 17
Potenz = M40 : M$ = "M40=" + M40$ : Z$ = Z40$
Case 18
Potenz = M41 : M$ = "M41=" + M41$ : Z$ = Z41$
Case 19
Potenz = M42 : M$ = "M42=" + M42$ : Z$ = Z42$
Case 20
Potenz = M43 : M$ = "M43=" + M43$ : Z$ = Z43$
Case 21
Potenz = M44 : M$ = "M44=" + M44$ : Z$ = Z44$
Case 22
Potenz = M45 : M$ = "M45=" + M45$ : Z$ = Z45$
Case 23
Potenz = M46 : M$ = "M46=" + M46$ : Z$ = Z46$
Case 24
Potenz = M47 : M$ = "M47=" + M47$ : Z$ = Z47$
Case 25
Potenz = M48 : M$ = "M48=" + M48$ : Z$ = Z48$
Case 26
Potenz = M49 : M$ = "M49=" + M49$ : Z$ = Z49$
Case 27
Potenz = M50 : M$ = "M50=" + M50$ : Z$ = Z50$
Case 28
Info0$ + M$ + ": "
SetWindowTitle(0, Info0$)
Break
EndSelect
EndIf
ForEver
For i = 1 To 28
FreeGadget(i)
Next
Calc()
ResLen.q = StringByteLength(Result$)
PokeB(@Result$ + ResLen - Asc_Uni, PeekB(@Result$ + ResLen - Asc_Uni) - 1) ;-1 für Mersenne
;=========================================================================
Case 5 ;Fakultät
Info0$ = "Fakultät: "
SetWindowTitle(0, Info0$)
TextGadget(1, 50, 10, 300, 20, "Fakultät:")
StringGadget(2, 50, 50, 300, 20, "2", #PB_String_Numeric)
ButtonGadget(3, 300, 560, 200, 30, "Berechnung starten")
SetActiveGadget(2)
For i = 1 To 3 ;das StringGadget einfach mal drin lassen
SetGadgetFont(i, FontID(0))
Next
Repeat
Repeat
Event = WaitWindowEvent()
If Event = #PB_Event_CloseWindow
CloseWindow(0)
End
EndIf
Until EventGadget() = 3
Fak = Val(GetGadgetText(2)) ;hier reicht der Wertebereich von Val() wohl aus ;-)
;besser doch Test
If (Len(GetGadgetText(2)) > 7) Or (Fak > 1E6) ;willkürlich gewählte Werte, sichere Seite
MessageRequester("Abbruch!", "Wert zu groß! Bitte neue Eingabe!")
Else
Break
EndIf
ForEver
Info0$ + Str(Fak) + "! "
SetWindowTitle(0, Info0$)
Result$ = "1"
For i = 1 To 3
FreeGadget(i)
Next
ListViewGadget(2, 10, 10, 250, 580) ;Info-Fenster
AddWindowTimer(0, 0, 5000) ;für Anzeige im Info-Fenster
BindEvent(#PB_Event_Timer, @FakInfo())
GadgetZeile = 0
AddGadgetItem (2, -1, "Berechnung gestartet... ")
SetGadgetState(2, GadgetZeile)
ZeitA = ElapsedMilliseconds()
For FakA = 1 To Fak
Faktor1$ = Str(FakA)
WindowEvent() ;nur für Anzeige im Info-Fenster
Result$ = FFT_Mul(@Faktor1$, @Result$)
Next
ZeitE = ElapsedMilliseconds() - ZeitA
RemoveWindowTimer(0, 0)
UnbindEvent(#PB_Event_Timer, @FakInfo()) ;nötig!
FreeGadget(2)
;=========================================================================
Case 6 ;Neuner-Test
Info0$ = "Neuner-Test: "
SetWindowTitle(0, Info0$)
Test9 = 540000 ;540000 ist o.K. für Double
Buffer9 = AllocateMemory(32 * Test9 * #Asc_Uni + 64)
Buffer9A = Buffer9 + (64 - (Buffer9 & $3F))
!mov rax,[v_Buffer9A]
!mov byte[rax],9
!mov rcx,[v_Test9]
!cmp [v_Asc_Uni],1
!je .ASCII
!vpbroadcastw ymm0,[rax] ;Unicode Word
!shl rcx,1
!jmp @f
!.ASCII:
!vpbroadcastb ymm0,[rax] ;ASCII Byte
!@@:
!vmovdqa [rax],ymm0
!add rax,32
!dec rcx
!jnz @b
Info0$ + "Berechne " + Str(32 * Test9) + "*" + Str(32 * Test9) + " Neuner-Bytes (Digits): "
SetWindowTitle(0, Info0$)
ZeitA = ElapsedMilliseconds()
Result$ = FFT_Mul(Buffer9A, Buffer9A)
ZeitE = ElapsedMilliseconds() - ZeitA
Info9$ = ""
L = StringByteLength(Result$)
Buffer1 = @Result$
If PeekB(Buffer1 + L - Asc_Uni) <> $31 ;hier könnte auch die Konstante genommen werden
Info9$ = "Fehler 1! "
EndIf
If PeekB(Buffer1 + (L / 2) - Asc_Uni) <> $38
Info9$ + "Fehler 8! "
EndIf
For i = 0 To (L / 2) - 2 * Asc_Uni Step #Asc_Uni
If PeekB(Buffer1 + i) <> $39
Info9$ + "Fehler 9! "
Break
EndIf
Next
For i = L / 2 To L - 2 * Asc_Uni Step #Asc_Uni
If PeekB(Buffer1 + i) <> $30
Info9$ + "Fehler 0! "
Break
EndIf
Next
If Info9$ = ""
Info9$ = "Alles o.K.! "
EndIf
FreeMemory(Buffer9)
Info0$ + Info9$
SetWindowTitle(0, Info0$)
;=========================================================================
EndSelect
;=========================================================================
EditorGadget(1, 10, 10, 780, 540)
Info1$ = Info0$ + "Benötigte Zeit: " + Str(ZeitE) + "ms " + Z$
SetWindowTitle(0, Info1$)
Digits$ = ""
If Asc(Mid(Result$, 1, 1)) < 58 ;also Ziffer (hoffentlich :-) !) und nicht z.B.Fehlermeldung
Digits$ = "Digits:" + Str(Len(Result$)) + #LFCR$
EndIf
SetGadgetText(1, Digits$ + Result$)
ButtonGadget(2, 100, 560, 250, 30, "In die Zwischenablage kopieren")
SetGadgetFont(2, FontID(0))
ButtonGadget(3, 450, 560, 250, 30, "Zurück zur Auswahl")
SetGadgetFont(3, FontID(0))
Repeat
Event = WaitWindowEvent()
If Event = #PB_Event_CloseWindow
Break 2
EndIf
If Event = #PB_Event_Gadget
Select EventGadget()
Case 2
SetClipboardText(Result$)
Case 3
SetWindowTitle(0, FensterTitel$)
Info0$ = "" : Info1$ = "" : Z$ = ""
Break
EndSelect
EndIf
ForEver
For i = 1 To 3
FreeGadget(i)
Next
Until WaitWindowEvent() = #PB_Event_CloseWindow
CloseWindow(0)
End
Procedure.s FFT_Mul(PFac1, PFac2)
EnableASM
Protected.q BufferCosSin, BufferCosSinA, BufferDez, BufferDezA, BufferFX, BufferFXA, BufferFaktor1, BufferFaktor1A
Protected.q BufferFaktor2, BufferFaktor2A, BufferProdukt, BufferProduktA, BufferReverse, BufferReverseA, LG, N
Protected.s sAusgabe
;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 = "-1"
ProcedureReturn sAusgabe
EndIf
PUSH BufferFXA
!pop [BufferFXA]
;!Align 16
;Save Registers
!mov rax,[BufferFXA]
!fxsave64 [rax] ;sichert erstmal komplett FPU (mit MXCSR) und XMM (0-15). Lass ich mal so, obwohl FNSAVE auch reichen würde. Dann aber so keine Register-Sicherung
!mov [rax+464],rdi ;sichert callee-save registers, ohne RBX und RBP!
!mov [rax+472],rsi
!mov [rax+480],r12
!mov [rax+488],r13
!mov [rax+496],r14
!mov [rax+504],r15
;übergebende Parameter
!mov r8,[p.v_PFac1]
!mov [PF1],r8
!mov r9,[p.v_PFac2]
!mov [PF2],r9
!mov rax, -16
!mov rdx, -16
!pxor xmm1,xmm1 ;xmm1 null setzen, ist praktisch das Zero-Byte bzw. 2 Zero-Bytes für Unicode
;Länge der Strings ermitteln
CompilerIf #PB_Compiler_Unicode
!@@:
!add rax,16
!vpcmpistri xmm1,dqword[r8+rax],00001001b ;Bit0 = 1 und Bit1 = 0, also String-Chars werden als unsigned Words interpretiert, Bit2 = 0 und Bit3 = 1, also Test auf equal each
!jnz @b
!shr rax,1 ;Pointer in String muss hier halbiert werden
!add rax,rcx
!pxor xmm1,xmm1 ;xmm1 null setzen, ist praktisch das Zero-Byte bzw. 2 Zero-Bytes für Unicode
!@@:
!add rdx,16
!vpcmpistri xmm1,dqword[r9+rdx],00001001b ;Bit0 = 1 und Bit1 = 0, also String-Chars werden als unsigned Words interpretiert, Bit2 = 0 und Bit3 = 1, also Test auf equal each
!jnz @b
!shr rdx,1 ;Pointer in String muss hier halbiert werden
!add rdx,rcx
CompilerElse
!@@:
!add rax,16
!vpcmpistri xmm1,dqword[r8+rax],00001000b ;Bit0 = 0 und Bit1 = 0, also String-Chars werden als unsigned Bytes interpretiert, Bit2 = 0 und Bit3 = 1, also Test auf equal each
!jnz @b
!add rax,rcx
!pxor xmm1,xmm1 ;xmm1 null setzen, ist praktisch das Zero-Byte bzw. 2 Zero-Bytes für Unicode
!@@:
!add rdx,16
!vpcmpistri xmm1,dqword[r9+rdx],00001000b ;Bit0 = 0 und Bit1 = 0, also String-Chars werden als unsigned Bytes interpretiert, Bit2 = 0 und Bit3 = 1, also Test auf equal each
!jnz @b
!add rdx,rcx
CompilerEndIf
!mov [L1],rax
!mov [L2],rdx
!cmp rax,rdx ;L1 soll der größere Faktor sein
!jae @f
!mov [L1],rdx
!mov [L2],rax
!mov [PF1],r9
!mov [PF2],r8
!@@:
!add rax,rdx
!push rax
POP LG ;Gesamtlänge L1+L2 für "Verstringung" und Speicher-Management
;LenF = L1
!mov rdx,[L1] ;L1 ist größter Faktor
!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 = "Fehler! Mindestens einer der Faktoren hat Länge Null! " ;wegen Kompatibilität zu "normaler" Multiplikation gelassen
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,5 ;Produkt-Länge muss mindestens 32 sein, sonst an anderen Stellen Klimmzüge ("Verstringung")
!sub rcx,[Exponent]
!js @f ;Exponent ist > 5 (Hauptfall)
!jz @f ;Exponent ist = 5
!shl [LenF],cl
!add [Exponent],rcx
!@@:
;zur Sicherheit im MXCSR Bit 13 und 14 auf Null setzen (round to nearest)
;MXCSR[6] : Denormals are zeros - 0
;MXCSR[7:12] : Exception masks all 1's (all exceptions masked)
;MXCSR[13:14] : Rounding control - 0 (round to nearest)
;MXCSR[15] : Flush to zero for masked underflow - 0 (off)
;R+ bit 14 Round Positive
;R- bit 13 Round Negative
;RZ bits 13 and 14 Round to Zero
;RN bits 13 and 14 are 0 Round to Nearest
;load MXCSR
!vstmxcsr [MXCSR] ;load Control-and Status-Register in Memory
!mov eax,[MXCSR]
!mov [MXCSR_Sicher],eax ;trotz fxsave64 (s.o.) hier zur Sicherheit
;MXCSR_Sicher = MXCSR
!and [MXCSR],11111111111111111001111111111111b ;zur Sicherheit Bit 13 und 14 auf Null setzen (round to nearest)
!vldmxcsr [MXCSR] ;store MXCSR
;N = LenF << 1
!mov rax,[LenF]
!shl rax,1
!mov [N],rax
!push rax ;wird dann "PB-N"
!vcvtsi2sd xmm1,xmm2,rax ;Quad-N in Double-ND konvertieren für spätere Division
!vmovsd [ND],xmm1
;Pointer_Produkt.q = N - 1
!dec rax
!mov [Pointer_Produkt],rax
POP N ;wird u.a.für AllocateMemory benötigt!
;Speicher-Reservierung
BufferCosSin = AllocateMemory((2 * N * 8) + 64)
If BufferCosSin
BufferCosSinA = BufferCosSin + (64 - (BufferCosSin & $3F)) ;64-er Alignment, klotzen, nicht kleckern!
Else
sAusgabe = "-1"
ProcedureReturn sAusgabe
EndIf
PUSH BufferCosSinA
!pop [BufferCosSinA]
BufferFaktor1 = AllocateMemory((2 * N * 8) + 64)
If BufferFaktor1
BufferFaktor1A = BufferFaktor1 + (64 - (BufferFaktor1 & $3F)) ;64-er Alignment
Else
sAusgabe = "-1"
ProcedureReturn sAusgabe
EndIf
PUSH BufferFaktor1A
!pop [BufferFaktor1A]
BufferFaktor2 = AllocateMemory((2 * N * 8) + 64)
If BufferFaktor2
BufferFaktor2A = BufferFaktor2 + (64 - (BufferFaktor2 & $3F)) ;64-er Alignment
Else
sAusgabe = "-1"
ProcedureReturn sAusgabe
EndIf
PUSH BufferFaktor2A
!pop [BufferFaktor2A]
;========================================================
;Faktoren bearbeiten
!tzcnt rax,[N] ;BMI1, auch LZCNT möglich
!mov [Exponent],rax
!mov rdx,32
!sub rdx,rax
!vmovq xmm9,rdx ;32 - Exponent
!mov rsi,[N]
!mov r8,[PF1]
!mov r10,[PF2]
!mov r11,[BufferCosSinA] ;Temp-Memory
!mov r9,r11
!add r9,rsi
!add r9,rsi
!mov r14,[L1]
!mov rcx,[v_Asc_Uni]
!dec rcx
!shl r14,cl
!sub r14,[v_Asc_Uni]
!mov r15,[L2]
!shl r15,cl
!sub r15,[v_Asc_Uni]
!mov rcx,2
!sub rcx,[v_Asc_Uni]
!shr rsi,cl
!sub rsi,[v_Asc_Uni]
!mov rcx,7
!@@:
!mov [r11+rcx*4],ecx
!dec rcx
!jns @b
!vmovdqa ymm0,[r11] ;Start-Values 0,1,2,3,4,5,6,7
!vpbroadcastd ymm1,[V1] ;AVX2
!vpbroadcastd ymm2,[V2]
!vpbroadcastd ymm3,[V3]
!vpbroadcastd ymm4,[V4]
!vpbroadcastd ymm11,[VAcht]
!mov rdi,[BufferFaktor1A] ;Temp-Memory
!xor r12,r12
;------------------------------------
!.Loop_Rev_Fak:
!vpsrld ymm5,ymm0,1 ;>> 1
!vpand ymm6,ymm5,ymm1 ;& V1
!vpand ymm7,ymm0,ymm1 ;& V1
!vpslld ymm8,ymm7,1 ;<< 1
!vpor ymm10,ymm6,ymm8 ;|
!vpsrld ymm5,ymm10,2 ;>> 2
!vpand ymm6,ymm5,ymm2 ;& V2
!vpand ymm7,ymm10,ymm2 ;& V2
!vpslld ymm8,ymm7,2 ;<< 2
!vpor ymm10,ymm6,ymm8 ;|
!vpsrld ymm5,ymm10,4 ;>> 4
!vpand ymm6,ymm5,ymm3 ;& V3
!vpand ymm7,ymm10,ymm3 ;& V3
!vpslld ymm8,ymm7,4 ;<< 4
!vpor ymm10,ymm6,ymm8 ;|
!vpsrld ymm5,ymm10,8 ;>> 8
!vpand ymm6,ymm5,ymm4 ;& V4
!vpand ymm7,ymm10,ymm4 ;& V4
!vpslld ymm8,ymm7,8 ;<< 8
!vpor ymm10,ymm6,ymm8 ;|
!vpsrld ymm5,ymm10,16 ;>> 16
!vpslld ymm8,ymm10,16 ;<< 16
!vpor ymm10,ymm5,ymm8 ;|
!vpsrld ymm5,ymm10,xmm9 ;>> (32 - Exponent)
!vmovdqa [rdi],ymm5
!xor r13,r13
!.Next_Byte8:
!mov r12d,[rdi+r13*4] ;8 reverse Positionen
!xor rax,rax
!mov rdx,r14
!mov rcx,[v_Asc_Uni]
!dec rcx
!shr rdx,cl
!cmp [L1],rdx
!jb @f
!movzx eax,byte[r8+r14] ;next Byte Faktor1$
!and eax,0fh ;for String (=sub 30h)
!@@:
!mov [r9+r12*2],eax
!xor rax,rax
!mov rdx,r15
!mov rcx,[v_Asc_Uni]
!dec rcx
!shr rdx,cl
!cmp [L2],rdx
!jb @f
!movzx eax,byte[r10+r15] ;next Byte Faktor2$
!and eax,0fh
!@@:
!mov [r11+r12*2],eax
!sub rsi,[v_Asc_Uni]
!sub r14,[v_Asc_Uni]
!sub r15,[v_Asc_Uni]
!inc r13
!cmp r13,8
!jne .Next_Byte8
!vpaddd ymm0,ymm0,ymm11
!or rsi,rsi
!jns .Loop_Rev_Fak
;========================================================
;Integer in Double konvertieren, verdoppeln und abspeichern
!xor rdx,rdx
!xor r8,r8
!mov rcx,[N]
!shr rcx,3
!mov r10,[BufferFaktor1A]
!mov r11,[BufferCosSinA]
!mov r12,[BufferFaktor2A]
!mov r9,r11
!mov rsi,[N]
!add r9,rsi
!add r9,rsi
!@@:
;Faktor1
!vcvtdq2pd ymm0,[r9+r8] ;4 Ziffern (als DWords) von Faktor1 in 4 Doubles konvertieren
!vmovsd [r10+rdx],xmm0
!vmovsd [r10+rdx+16],xmm0 ;+8 ist der I-Anteil, bleibt hier Null!
!vpsrldq ymm1,ymm0,8
!vmovsd [r10+rdx+32],xmm1
!vmovsd [r10+rdx+48],xmm1
!vperm2f128 ymm2,ymm0,ymm0,1
!vmovsd [r10+rdx+64],xmm2
!vmovsd [r10+rdx+80],xmm2
!vpsrldq ymm1,ymm2,8
!vmovsd [r10+rdx+96],xmm1
!vmovsd [r10+rdx+112],xmm1
;Faktor2
!vcvtdq2pd ymm0,[r11+r8]
!vmovsd [r12+rdx],xmm0
!vmovsd [r12+rdx+16],xmm0
!vpsrldq ymm1,ymm0,8
!vmovsd [r12+rdx+32],xmm1
!vmovsd [r12+rdx+48],xmm1
!vperm2f128 ymm2,ymm0,ymm0,1
!vmovsd [r12+rdx+64],xmm2
!vmovsd [r12+rdx+80],xmm2
!vpsrldq ymm1,ymm2,8
!vmovsd [r12+rdx+96],xmm1
!vmovsd [r12+rdx+112],xmm1
!add r8,16
!add rdx,128
!sub rcx,1
!jnz @b
;========================================================
;Winkel erst hier, weil BufferCosSinA oben "missbraucht" wird
!vcvtsi2sd xmm0,xmm0,qword[N]
!vmovsd xmm1,[Pi2]
!vdivsd xmm7,xmm1,xmm0 ;Schrittweite Konstante
!vmovq xmm6,qword[Signum] ;Konstante
!mov rcx,[N]
!shr rcx,3
!xor r8,r8 ;Schritt-Nr. Konstante
!mov r9,[N]
!shr r9,2 ;N/4 für Adressierung unten
!mov r10,[N]
!shr r10,1 ;N/2
!mov r11,[BufferCosSinA]
!inc rcx
!@@:
!vcvtsi2sd xmm15,xmm14,r8
!vmulsd xmm0,xmm7,xmm15
!vbroadcastsd ymm1,xmm0 ;YMM1: 1=x^1 2=x^1 3=x^1 4=x^1 AVX2, da Quelle Register ist
!vmulpd ymm2,ymm1,ymm1 ;YMM2: 1=x^2 2=x^2 3=x^2 4=x^2
!vmulsd xmm3,xmm2,xmm2 ;YMM3: 1=x^4 2=x^2 3=0 4=0
!vmulpd ymm4,ymm2,ymm2 ;YMM4: 1=x^4 2=x^4 3=x^4 4=x^4
!vmulpd ymm2,ymm4,ymm4 ;YMM2: 1=x^8 2=x^8 3=x^8 4=x^8 for the next 4 terms
!vmovdqa xmm8,xmm3 ;save for Cos
!vmulpd xmm9,xmm8,xmm4
!vperm2f128 ymm10,ymm9,ymm8,100000b
!vmulpd xmm14,xmm3,xmm1 ;YMM14: 1=x^5 2=x^3 3=0 4=0
!vmulpd xmm1,xmm14,xmm4 ;YMM1: 1=x^9 2=x^7 3=0 4=0
!vperm2f128 ymm5,ymm1,ymm14,100000b ;YMM5: 1=x^9 2=x^7 3=x^5 4=x^3 YMM5(0-127)=YMM1(0-127), YMM5(128-255)=YMM14(0-127)
!vmulpd ymm1,ymm5,yword[RezFakSin] ;multiply the 4 values in YMM5 with RezFak -1/3! ... 1/9!, result in YMM1
!vmulpd ymm11,ymm10,yword[RezFakCos]
;next 4 terms, without loop
!vmulpd ymm14,ymm5,ymm2 ;YMM14: 1=x^17 2=x^15 3=x^13 4=x^11
!vmulpd ymm10,ymm10,ymm2
!vfmadd231pd ymm1,ymm14,yword[RezFakSin+32]
!vfmadd231pd ymm11,ymm10,yword[RezFakCos+32]
;next 4 terms
!vmulpd ymm5,ymm14,ymm2 ;YMM5: 1=x^25 2=x^23 3=x^21 4=x^19
!vmulpd ymm10,ymm10,ymm2
!vfmadd231pd ymm1,ymm5,yword[RezFakSin+64]
!vfmadd231pd ymm11,ymm10,yword[RezFakCos+64]
;result
!vhaddpd ymm2,ymm1,ymm1 ;YMM2: 1=1+2 of YMM1, 3=3+4 of YMM1
!vhaddpd ymm12,ymm11,ymm11
!vextractf128 xmm1,ymm2,1b ;XMM1: 3+4 of YMM2
!vextractf128 xmm11,ymm12,1b
!vaddsd xmm3,xmm2,xmm1 ;XMM3: 1=sum of iterations
!vaddsd xmm13,xmm12,xmm11
!vaddsd xmm12,xmm0,xmm3 ;Sine XMM12: 1=sum of iterations plus start-value (1.term=x)
!vaddsd xmm14,xmm13,[Eins] ;Cosine
;------------------------------------
;store values 0-180°
!mov rax,r8
!shl rax,1
!vmovsd qword[r11+rax*8+8],xmm12 ;Sin 0-45°
!vorpd xmm1,xmm12,xmm6 ;neg Sin
!vmovsd qword[r11+rax*8],xmm14 ;Cos 0-45°
!mov rax,r9
!sub rax,r8
!shl rax,1
!vmovsd qword[r11+rax*8],xmm12 ;Cos 45-90°
!mov rdx,r9
!add rdx,r8
!shl rdx,1
!vmovsd qword[r11+rdx*8],xmm1 ;Cos 90-135°
!mov r13,r10
!sub r13,r8
!shl r13,1
!vorpd xmm2,xmm14,xmm6
!vmovsd qword[r11+r13*8],xmm2 ;Cos 135-180°
!vmovsd qword[r11+rax*8+8],xmm14 ;Sin 45-90°
!vmovsd qword[r11+rdx*8+8],xmm14 ;Sin 90-135°
!vmovsd qword[r11+r13*8+8],xmm12 ;Sin 135-180°
!inc r8
!dec rcx
!jnz @b
;========================================================
;jetzt die beiden Faktoren transformieren
!mov rsi,[N]
!mov r8,[BufferCosSinA]
!mov rcx,1 ;Exponent_FFT
!mov r12,2 ;Pointer0
!mov r9,2 ;Pointer2
!mov r10,[BufferFaktor1A]
!mov rdi,[BufferFaktor2A]
!.OuterLoopF:
!inc rcx ;Exponent_FFT
!mov rax,rsi ;N
!shr rax,cl
!jz .EndF
!mov r14,rax
!.InnerLoopF:
!xor r11,r11 ;Pointer1
!mov r15,r12 ;Pointer0
;------------------------------------
!@@:
!mov rax,r14
!mul r11 ;Pointer1
!shl rax,3 ;EW*8, Pointer in CosSin
!mov rdx,r9 ;Pointer2
!mov r13,rdx
!shl rdx,4 ;Pointer2*8
!sub r13,r12 ;[v_Pointer0]
!shl r13,4 ;(Pointer2 - Pointer0) * 8
!vperm2f128 ymm0,ymm9,yword[r8+rax],100010b ;2x CosSinA in ymm0, ymm9 is Dummy
!vmovapd xmm7,dqword[rdi+rdx] ;BufferFaktor2A complete (R and I)
!vperm2f128 ymm2,ymm7,yword[r10+rdx],10b
!vmulpd ymm4,ymm2,ymm0
!vshufpd ymm1,ymm0,ymm0,101b ;change Sin and Cos
!vhsubpd ymm3,ymm4,ymm4 ;WR
!vmulpd ymm4,ymm2,ymm1
!vhaddpd ymm5,ymm4,ymm4 ;WI
!vshufpd ymm15,ymm3,ymm5,0b ;WI in ymm15-high
!vmovapd xmm1,dqword[r10+r13]
!vmovapd xmm7,dqword[rdi+r13]
!vperm2f128 ymm2,ymm7,ymm1,10b
!vaddpd ymm14,ymm15,ymm2
!vsubpd ymm13,ymm2,ymm15
!vmovapd [r10+r13],xmm14
!vperm2f128 ymm2,ymm7,ymm14,11b
!vmovapd [rdi+r13],xmm2
!vmovapd [r10+rdx],xmm13
!vperm2f128 ymm2,ymm7,ymm13,11b
!vmovapd [rdi+rdx],xmm2
!add r11,2 ;Pointer1
!inc r9 ;Pointer2
!dec r15
!jnz @b
!add r9,r12 ;Pointer0
!cmp r9,rsi ;N
!jb .InnerLoopF
!shl r12,1
!mov r9,r12
!jmp .OuterLoopF
!.EndF:
;========================================================
;Und jetzt die beiden transformierten Vektoren miteinander multiplizieren
!mov r8,[BufferFaktor1A] ;nimmt auch das Produkt auf
!mov r9,[BufferFaktor2A]
!xor r11,r11 ;Pointer
!mov rcx,[N]
!shr rcx,1
!vmovdqu ymm5,[Null1] ;Null,Signum,Null,Signum
!@@:
!vmovapd ymm0,[r8+r11]
!vmovapd ymm3,[r9+r11]
!vmulpd ymm1,ymm0,ymm3 ;ymm1:R*R,I*I,R*R,I*I
!vshufpd ymm2,ymm3,ymm3,101b ;change R and I
!vxorpd ymm4,ymm1,ymm5 ;Vorzeichen ändern, damit nicht sub sondern add für beide Terme verwendet werden kann
!vmulpd ymm1,ymm0,ymm2 ;R*I, R*I
!vhaddpd ymm6,ymm4,ymm1
!vmovapd [r8+r11],ymm6
!add r11,32
!dec rcx
!jnz @b
;========================================================
;Rück-Transformation
;Produkt bearbeiten
BufferReverse = AllocateMemory(32 + 64)
If BufferReverse
BufferReverseA = BufferReverse + (64 - (BufferReverse & $3F)) ;64-er Alignment
Else
sAusgabe = "-1"
ProcedureReturn sAusgabe
EndIf
PUSH BufferReverseA
!pop [BufferReverseA]
!mov rsi,[N] ;für Abbruch Schleife Loop_Rev_Prod; letzte Pos.steht ja fest (wie erste)
!shl rsi,1
!sub rsi,2
!mov r8,[BufferFaktor1A] ;Ergebnis von Multiplikation oben
!mov r9,[BufferFaktor2A] ;hier neue Reihenfolge rein
!mov rdx,32
!sub rdx,[Exponent]
!vmovq xmm9,rdx
!mov rcx,7
!@@:
!mov [r9+rcx*4],ecx
!dec rcx
!jns @b
!vmovdqa ymm0,[r9] ;Startwerte (8 DWords=32 Bytes) 0,1,2,3,4,5,6,7 r9 hier "missbraucht"
!vpbroadcastd ymm1,[V1] ;AVX2
!vpbroadcastd ymm2,[V2]
!vpbroadcastd ymm3,[V3]
!vpbroadcastd ymm4,[V4]
!vpbroadcastd ymm11,[VAcht]
!mov rdi,[BufferReverseA]
!xor r12,r12
;------------------------------------
!.Loop_Rev_Prod:
!vpsrld ymm5,ymm0,1 ;>> 1
!vpand ymm6,ymm5,ymm1 ;& V1
!vpand ymm7,ymm0,ymm1 ;& V1
!vpslld ymm8,ymm7,1 ;<< 1
!vpor ymm10,ymm6,ymm8 ;|
!vpsrld ymm5,ymm10,2 ;>> 2
!vpand ymm6,ymm5,ymm2 ;& V2
!vpand ymm7,ymm10,ymm2 ;& V2
!vpslld ymm8,ymm7,2 ;<< 2
!vpor ymm10,ymm6,ymm8 ;|
!vpsrld ymm5,ymm10,4 ;>> 4
!vpand ymm6,ymm5,ymm3 ;& V3
!vpand ymm7,ymm10,ymm3 ;& V3
!vpslld ymm8,ymm7,4 ;<< 4
!vpor ymm10,ymm6,ymm8 ;|
!vpsrld ymm5,ymm10,8 ;>> 8
!vpand ymm6,ymm5,ymm4 ;& V4
!vpand ymm7,ymm10,ymm4 ;& V4
!vpslld ymm8,ymm7,8 ;<< 8
!vpor ymm10,ymm6,ymm8 ;|
!vpsrld ymm5,ymm10,16 ;>> 16
!vpslld ymm8,ymm10,16 ;<< 16
!vpor ymm10,ymm5,ymm8 ;|
!vpsrld ymm5,ymm10,xmm9 ;>> (32 - Exponent)
!vmovdqa [rdi],ymm5
!xor rcx,rcx
!@@:
!mov r12d,[rdi+rcx*4]
!shl r12,1
!vmovapd xmm12,[r8] ;16 Bytes (R.d und I.d)
!vmovapd [r9+r12*8],xmm12
!add r8,16
!inc rcx
!cmp rcx,8
!jne @b
!vpaddd ymm0,ymm0,ymm11
!cmp r12,rsi ;Position letzter Wert?
!jne .Loop_Rev_Prod
;========================================================
!mov r8,[BufferFaktor2A]
!mov rcx,[N]
!vpbroadcastq ymm4,qword[Signum] ;Signum kann auch hier für Gather verwendet werden
!vmovdqu ymm6,[Adresse2] ;0,8,32,40 auch add abspeichern
!vmovdqu ymm7,[Adresse3] ;16,24,48,56 auch sub abspeichern
!vmovdqa ymm5,ymm4 ;ymm4 sichern, wird auf Null gesetzt!
!shr rcx,2
!@@:
!vgatherqpd ymm1,[r8+ymm6],ymm4 ;ymm4 s.o.
!vmovdqa ymm4,ymm5
!vgatherqpd ymm2,[r8+ymm7],ymm4
!vaddpd ymm8,ymm1,ymm2
!vmovapd [r8],xmm8
!vextractf128 [r8+32],ymm8,1b
!vsubpd ymm9,ymm1,ymm2
!vmovapd [r8+16],xmm9
!vextractf128 [r8+48],ymm9,1b
!vmovdqa ymm4,ymm5
!add r8,64
!dec rcx
!jnz @b
;------------------------------------
;------------------------------------
!mov rsi,[N]
!mov r8,[BufferCosSinA]
!mov rcx,1 ;Exponent_FFT
!mov r12,2 ;Pointer0
!mov r9,2 ;Pointer2
!mov r10,[BufferFaktor2A]
!.OuterLoopP:
!inc rcx ;Exponent_FFT
!mov rax,rsi ;N
!shr rax,cl
!jz .EndP
!mov r14,rax
!shl r14,3 ;EW*8
!.InnerLoopP:
!xor r11,r11 ;Pointer1
!mov r15,r12
;------------------------------------
!@@:
!mov rax,r14
!mul r11 ;Pointer1 oder imul rax,r11 ;dann wäre rdx frei
!mov rdx,r9 ;Pointer2
!mov r13,rdx
!shl rdx,4
!sub r13,r12 ;[v_Pointer0]
!shl r13,4
!vmovapd xmm7,dqword[r10+r13]
!vmovapd xmm0,dqword[r10+rdx] ;BufferFaktor2A
!vmovapd xmm1,dqword[r8+rax] ;CosSinA
!vshufpd xmm2,xmm0,xmm0,1b ;R und I vertauschen
!vmulpd xmm3,xmm2,xmm1
!vmulpd xmm4,xmm0,xmm1
!vxorpd xmm5,xmm3,[Null1] ;Signum, xmm3 negieren um vhaddpd verwenden zu können (sonst sub) wie xmm4
!vhaddpd xmm6,xmm4,xmm5 ;WI in xmm15-high, WR in Low
!vaddpd xmm14,xmm7,xmm6 ;dqword[r10+r13]
!vsubpd xmm13,xmm7,xmm6
!vmovapd [r10+r13],xmm14
!vmovapd [r10+rdx],xmm13
!add r11,2 ;Pointer1
!inc r9 ;Pointer2
!dec r15
!jnz @b
;---------------------------------------
!add r9,r12 ;Pointer0
!cmp r9,rsi ;N
!jb .InnerLoopP
!shl r12,1
!mov r9,r12
!jmp .OuterLoopP
!.EndP:
;========================================================
FreeMemory(BufferReverse)
FreeMemory(BufferCosSin)
FreeMemory(BufferFaktor1)
BufferProdukt = AllocateMemory(N + 64)
If BufferProdukt
BufferProduktA = BufferProdukt + (64 - (BufferProdukt & $3F)) ;64-er Alignment
Else
sAusgabe = "-1"
ProcedureReturn sAusgabe
EndIf
PUSH BufferProduktA
!pop [BufferProduktA]
BufferDez = AllocateMemory(16 + 64)
If BufferDez
BufferDezA = BufferDez + (64 - (BufferDez & $3F)) ;64-er Alignment
Else
sAusgabe = "-1"
ProcedureReturn sAusgabe
EndIf
PUSH BufferDezA
!pop [BufferDezA]
;========================================================
;Ergebnis für Ausgabe aufbereiten
;Konvertierung Double in Integer, Rundung gemäß Rounding Control Bits im MXCSR-Register, oben auf "nearest" gesetzt
;da BufferFaktor2A noch den I-Anteil enthält, wird hier "komprimiert", d.h. in die I-Bereiche werden die Produkte reinkopiert
;die Double-Werte in Long-Integer konvertieren, Division ist notwendig für höhere Genauigkeit; nicht später shiften!
!vpbroadcastq ymm4,qword[Signum] ;Signum kann auch hier für Gather verwendet werden
!vmovdqu ymm6,[Adresse1] ;Adressen der Real-Teile 0,16,32,48
!vmovdqa ymm5,ymm4 ;ymm4 wird auf Null gesetzt!
!mov r15,[N]
!mov r14,[BufferFaktor2A]
!shr r15,2
!mov rcx,r14
!vbroadcastsd ymm1,[ND] ;"nur" AVX, da Quelle Memory ist
!@@:
!vgatherqpd ymm0,[r14+ymm6],ymm4 ;ymm4 s.o.
!vdivpd ymm2,ymm0,ymm1
!vcvtpd2dq xmm3,ymm2
!vmovdqa [rcx],xmm3
!vmovdqa ymm4,ymm5
!add r14,64
!add rcx,16
!dec r15
!jnz @b
;Hex2Dez und "Verstringung"
;Pointer_Produkt = N - 1
;For i = 0 To N - 1 ;vom Produkt alle Werte durch
!mov r15,[N]
!mov r14,[BufferFaktor2A]
!mov rdx,[BufferProduktA]
!mov rdi,r15
!dec rdi
!lea r10,[Dezi]
!mov r11,[BufferDezA]
!.DeziLoop:
!mov eax,[r14] ;eax=Integer-Wert
;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,10 ;Anzahl Dezi-Werte ohne Null, 9 reicht für M49, ansonsten auch Dezi anpassen an Faktoren(Produkt)-Länge!
!.OuterLoopDez1:
;While (QuadWert - PeekQ(?Dezi + Pointer1)) >= 0
!mov r13d,eax
!.InnerLoopDez1:
!sub r13d,[r10+r8]
!js .InnerLoopDez1End
;Ziffer + 1
!inc r12b
;QuadWert - PeekQ(?Dezi + Pointer1)
;!sub rax,[r10+r8]
!mov eax,r13d ;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,10 ;Anzahl Dezi-Werte ohne Null, 9 reicht für M49, ansonsten auch Dezi anpassen an Faktoren(Produkt)-Länge!
;!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
!add r14,4
!dec r15
;Next
!jnz .DeziLoop
;------------------------
;Ergebnis (Produkt) als ASCII-Ziffern abspeichern
!vpbroadcastq ymm0,[AsciiNull] ;load 32x48 ("Null") AVX2
!mov r9,[BufferProduktA]
!xor r10,r10
!mov rcx,[N]
!shr rcx,5 ;/32, YMM-Length
!@@:
!vpor ymm1,ymm0,[r9+r10] ;AVX2
!vmovdqa [r9+r10],ymm1
!add r10,32
!dec rcx
!jnz @b
!mov byte[r9+r10],0 ;Zero-Byte für String, zur Sicherheit, evtl.word
FreeMemory(BufferFaktor2)
FreeMemory(BufferDez)
sAusgabe = LTrim(PeekS(BufferProduktA + N - LG, LG, #PB_Ascii), "0") ;evtl. führende Null(en) entfernen
If sAusgabe = "" ;mindestens einer der Faktoren war Null, Kompatibilität zu "normaler" Multiplikation
sAusgabe = "0"
EndIf
FreeMemory(BufferProdukt)
;Restore Registers
!mov rax,[BufferFXA]
!fxrstor64 [rax]
!mov rdi,[rax+464] ;restore callee-save registers
!mov rsi,[rax+472]
!mov r12,[rax+480]
!mov r13,[rax+488]
!mov r14,[rax+496]
!mov r15,[rax+504]
FreeMemory(BufferFX)
;Restore MXCSR ;nach fxrstor64!
!vldmxcsr [MXCSR_Sicher] ;schreibt den gesicherten Wert zurück ins MXCSR-Register
!vzeroupper
DisableASM
ProcedureReturn sAusgabe
;========================================================
DataSection
!RezFakSin: ;für Sinus
!dq 2.755731922398589065e-6 ; 1/9! 4.Iteration
!dq -1.984126984126984127e-4 ;-1/7! 3.Iteration
!dq 8.333333333333333333e-3 ; 1/5! 2.Iteration
!dq -1.666666666666666667e-1 ;-1/3! 1.Iteration
!dq 2.811457254345520763e-15 ; 1/17! 8.Iteration
!dq -7.647163731819816476e-13 ;-1/15! 7.Iteration
!dq 1.605904383682161460e-10 ; 1/13! 6.Iteration
!dq -2.505210838544171878e-8 ;-1/11! 5.Iteration
!dq 6.446950284384473396e-26 ; 1/25! 12.Iteration
!dq -3.868170170630684038e-23 ;-1/23! 11.Iteration
!dq 1.957294106339126123e-20 ; 1/21! 10.Iteration
!dq -8.220635246624329717e-18 ;-1/19! 9.Iteration
!RezFakCos: ;für Cosinus
!dq 2.480158730158730159e-5 ; 1/8! 4.Iteration
!dq -1.388888888888888889e-3 ;-1/6! 3.Iteration
!dq 4.166666666666666667e-2 ; 1/4! 2.Iteration
!dq -5.000000000000000000e-1 ;-1/2! 1.Iteration
!dq 4.779477332387385297e-14 ; 1/16! 8.Iteration
!dq -1.147074559772972471e-11 ;-1/14! 7.Iteration
!dq 2.087675698786809898e-9 ; 1/12! 6.Iteration
!dq -2.755731922398589065e-7 ;-1/10! 5.Iteration
!dq 1.611737571096118349e-24 ; 1/24! 12.Iteration
!dq -8.896791392450573287e-22 ;-1/22! 11.Iteration
!dq 4.110317623312164858e-19 ; 1/20! 10.Iteration
!dq -1.561920696858622646e-16 ;-1/18! 9.Iteration
!Eins:
!dq 1.0
!Pi2:
!dq 6.283185307179586476925
!Null1:
!dq 0h
!Signum:
!dq 8000000000000000h ;für neg Sin/Cos und Gather
!dq 0h
!dq 8000000000000000h
!AsciiNull:
!dq 3030303030303030h ;für String
!Dezi:
!dq 10000000000
!dq 1000000000 ;für M49 ist der größte Integer-Wert etwas über 200.000.000
!dq 100000000
!dq 10000000
!dq 1000000
!dq 100000
!dq 10000
!dq 1000
!dq 100
!dq 10
!dq 0 ;blanke Sicherheit
!Adresse1: ;für Gather
!dq 0
!dq 16
!dq 32
!dq 48
!Adresse2: ;für Gather
!dq 0
!dq 8
!dq 32
!dq 40
!Adresse3: ;für Gather
!dq 16
!dq 24
!dq 48
!dq 56
!BufferCosSinA dq ?
!BufferDezA dq ?
!BufferFaktor1A dq ?
!BufferFaktor2A dq ?
!BufferFXA dq ?
!BufferProduktA dq ?
!BufferReverseA dq ?
!Exponent dq ?
!L1 dq ?
!L2 dq ?
!LenF dq ?
!N dq ?
!ND dq ?
!PF1 dq ?
!PF2 dq ?
!Pointer_Produkt dq ?
!MXCSR dd ?
!MXCSR_Sicher dd ?
!V1 dd 55555555h
!V2 dd 33333333h
!V3 dd 0F0F0F0Fh
!V4 dd 00FF00FFh
!VAcht dd 8
EndDataSection
EndProcedure
Procedure Calc() ;Aufruf durch Potenzen (Case 3) und Mersenne (Case 4)
ZeitA = ElapsedMilliseconds()
Stellen = Potenz * Len(Result$) + 1
Stellen = Stellen + (64 - (Stellen & $3F))
!mov rax,[v_Potenz]
!bsr rdx,rax
!mov [v_Bits_M],rdx
CalcPot = Potenz ;CalcPot wird geshiftet!
BufferPot = AllocateMemory((Stellen * 2 * #Asc_Uni) + 64)
If BufferPot
BufferPotA = BufferPot + (64 - (BufferPot & $3F)) ;64-er Align
Else
MessageRequester("Fehler!", "Nicht genügend freier Arbeitsspeicher!", #MB_ICONERROR)
End
EndIf
ListViewGadget(2, 10, 10, 250, 580) ;Info-Fenster
PokeS(BufferPotA + Stellen * #Asc_Uni, Result$) ;hier noch Result$=Basis$, zum besseren Verständnis
PokeS(BufferPotA, "1") ;falls Bit0 von Potenz nicht gesetzt, X^0=1
!shr [v_CalcPot],1
!jnc @f
PokeS(BufferPotA, Basis$)
m = 1
!@@:
p = 1
j = 2
k = Bits_M - 1
z = 0 ;Zeile im ListViewGadget (Info)
For i = 0 To k
AddGadgetItem (2, -1, "Berechne " + Basis$ + "^" + Str(j))
SetGadgetState(2, z)
z + 1
Result$ = FFT_Mul(BufferPotA + Stellen * #Asc_Uni, BufferPotA + Stellen * #Asc_Uni) ;Potenzen berechnen
If Result$ = "-1"
MessageRequester("Fehler!", "Nicht genügend freier Arbeitsspeicher!", #MB_ICONERROR)
End
EndIf
PokeS(BufferPotA + Stellen * #Asc_Uni, Result$)
j * 2
!shr [v_CalcPot],1
!jnc @f
n + p
AddGadgetItem (2, -1, "Berechne " + Basis$ + "^" + Str(n * 2 + m))
SetGadgetState(2, z)
z + 1
Result$ = FFT_Mul(BufferPotA, BufferPotA + Stellen * #Asc_Uni) ;(Zwischen-) Ergebnis berechnen
If Result$ = "-1"
MessageRequester("Fehler!", "Nicht genügend freier Arbeitsspeicher!", #MB_ICONERROR)
End
EndIf
If i < k
PokeS(BufferPotA, Result$)
EndIf
!@@:
p * 2
Next
FreeMemory(BufferPot)
FreeGadget(2)
ZeitE = ElapsedMilliseconds() - ZeitA
EndProcedure
Procedure.i RangeCharSSE42(StrPointer, ChrPointer)
!mov rax,[p.v_ChrPointer]
!movdqu xmm1,[rax]
!mov rdx,[p.v_StrPointer]
!mov rax, -16
CompilerIf #PB_Compiler_Unicode ;Unicode
!align 16 ;kann speed-mäßig was bringen
!@@:
!add rax,16
!pcmpistri xmm1,dqword[rdx+rax],00000101b ;Bit0 = 1 und Bit1 = 0, also String-Chars werden als unsigned Words interpretiert, Bit2 = 1 und Bit3 = 0, also Test auf Ranges
!ja @b
!shr rax,1 ;Pointer in String muss hier halbiert werden (Word)
!add rax,rcx
!shl rcx,1 ;weil sonst nur 8 (Words), unten wird aber mit 16 verglichen
CompilerElse ;kein Unicode
!align 16 ;kann speed-mäßig was bringen
!@@:
!add rax,16
!pcmpistri xmm1,dqword[rdx+rax],00000100b ;Bit0 = 0 und Bit1 = 0, also String-Chars werden als unsigned Bytes interpretiert, Bit2 = 1 und Bit3 = 0, also Test auf Ranges
!ja @b
!add rax,rcx
CompilerEndIf
!cmp rcx,16 ;rcx=16 oder 8, wenn nichts gefunden! 8 schon erweitert
!je @f
!inc rax
ProcedureReturn ;rax = Pos., beginnt bei 1
!@@:
!xor rax,rax
ProcedureReturn ;rax = 0, Char nicht gefunden
EndProcedure
Procedure FakInfo()
GadgetZeile + 1
AddGadgetItem (2, -1, "Berechne aktuell " + Str(FakA) + "!")
SetGadgetState(2, GadgetZeile)
EndProcedure
Helle
Edit 7.1.2018: Kosmetik
Edit 11.1.2018: Kosmetik und minimale Speicherplatz-Reduzierung