Code: Alles auswählen
;- Schnelle Integer-Multiplikation mit FFT für große Zahlen
;- Beide Faktoren müssen gleich groß sein und die Länge muss eine Zweier-Potenz sein. Padding bei Bedarf mit führenden Nullen
;- Berechnungs-Schema läuft ab nach Schmetterlingsgraph
;- Keine Fremd-Bibliotheken!
;- Man beachte: Kein ASM-Anteil :-) !
;- Windows7/64, PureBasic 5.31 (x64)
;- "Helle" Klaus Helbing, 02.10.2015
;Als Quelle (Faktoren) dienen hier 2 Strings
Faktor1$ = "6582018229284824168619876730229402019930943462534319453394436096" ;64-Bytes lang für Test
Faktor2$ = "5486721320215405789450123459785103469751240504554109782315460125"
For i = 1 To 10 ;für Test. Faktoren-Länge von 2MB ist noch o.K.
Faktor1$ + Faktor1$
Faktor2$ + Faktor2$
Next
;Hier könnte Test auf Länge Faktor1$/Faktor2$ usw. hin
N.q = Len(Faktor1$) << 1
Exponent.q
Mantisse.q
PF1.q
PF2.q
EW.q
Pointer0.q
Pointer1.q
Pointer2.q
Pointer_Produkt.q = N - 1
Pointer_Produkt1.q
RevPos.q
FakPos.q
QuadWert.q
QuadWert1.q
Fehler.q ;nur Flag für Fehler-Test
Si.d
Co.d
Pi2.d = #PI * 2
Rad.d
Rad1.d
WR.d
WI.d
ZR.d
ZI.d
ND.d
WL.l
Ziffer.b
If N > 4194304 ;=2^22 als Produkt-Größe = 2^21 als Faktor sind noch o.K.
If MessageRequester("Warnung!", "Die Faktoren sind zu groß! Das Ergebnis wird mit hoher Wahrscheinlichkeit falsch sein!" + #CRLF$ + "Trotzdem weiter?", #PB_MessageRequester_YesNo) = #PB_MessageRequester_No
End
EndIf
EndIf
;Exponent ermitteln; 2^Exponent=N
ND = N
QuadWert = PeekQ(@ND) & $7FFFFFFFFFFFFFFF ;Vorzeichen-Bit ausblenden (Bit 63)
Exponent = QuadWert >> 52 ;Exponent = Bits 52 bis 62 (11 Bits)
Mantisse = QuadWert & $FFFFFFFFFFFFF ;Mantisse = Bits 0 bis 51 (52 Bits)
If Mantisse = $FFFFFFFFFFFFF ;also fast 1
Mantisse = 0 ;kann zu Null aufgerundet werden; ist ja Aufaddierung
QuadWert1 = 1 ;Exponent dafür um eins erhöhen
EndIf
If QuadWert = 0 Or Mantisse <> 0
Exponent = 0 ;Null
Else
Exponent = Exponent + QuadWert1 - 1023 ;1023=Bias
EndIf
If Exponent < 3
MessageRequester("Fehler!", "Möglicherweise Faktoren zu klein!")
End
EndIf
;Speicher-Reservierung; kann sicher noch verbessert werden. Vllt. mal alles zusammenfassen und dann zuteilen
BufferSin = AllocateMemory((N * 8) + 128)
If BufferSin
BufferSinA = BufferSin + (64 - (BufferSin & $3F)) ;64-er Alignment, klotzen, nicht kleckern!
Else
MessageRequester("Fehler!", "Nicht genügend Speicher für BufferSin!") ;kann auch noch ein "End" hinter
EndIf
BufferCos = AllocateMemory((N * 8) + 128)
If BufferCos
BufferCosA = BufferCos + (64 - (BufferCos & $3F)) ;64-er Alignmentr Alignment
Else
MessageRequester("Fehler!", "Nicht genügend Speicher für BufferCos!")
EndIf
BufferFaktor1 = AllocateMemory((2 * N * 8) + 128)
If BufferFaktor1
BufferFaktor1A = BufferFaktor1 + (64 - (BufferFaktor1 & $3F)) ;64-er Alignment
Else
MessageRequester("Fehler!", "Nicht genügend Speicher für BufferFaktor1!")
EndIf
BufferFaktor2 = AllocateMemory((2 * N * 8) + 128)
If BufferFaktor2
BufferFaktor2A = BufferFaktor2 + (64 - (BufferFaktor2 & $3F)) ;64-er Alignment
Else
MessageRequester("Fehler!", "Nicht genügend Speicher für BufferFaktor2!")
EndIf
BufferProdukt = AllocateMemory((N * 8) + 128)
If BufferProdukt
BufferProduktA = BufferProdukt + (64 - (BufferProdukt & $3F)) + 64 ;64-er Alignment, +64 um Unterlauf zu vermeiden
Else
MessageRequester("Fehler!", "Nicht genügend Speicher für BufferProdukt!")
EndIf
BufferDez.q = AllocateMemory(16 + 128)
If BufferDez
BufferDezA = BufferDez + (64 - (BufferDez & $3F)) ;64-er Alignment
Else
MessageRequester("Fehler!", "Nicht genügend Speicher für BufferDez!")
EndIf
;========================================================================================
TA_Gesamt = ElapsedMilliseconds()
PF1 = @Faktor1$
PF2 = @Faktor2$
TA_Faktoren = ElapsedMilliseconds()
;Faktoren bearbeiten
;Zuerst die Ziffern gemäß reverse Verteilung neu setzen, dabei gleich in Long konvertieren
FakPos = 0
For i = (N / 2) - 1 To 0 Step -1
RevPos = FakPos
RevPos = ((RevPos >> 1) & $55555555) | ((RevPos & $55555555) << 1)
RevPos = ((RevPos >> 2) & $33333333) | ((RevPos & $33333333) << 2)
RevPos = ((RevPos >> 4) & $0F0F0F0F) | ((RevPos & $0F0F0F0F) << 4)
RevPos = ((RevPos >> 8) & $00FF00FF) | ((RevPos & $00FF00FF) << 8)
RevPos = ( RevPos >> 16 ) | ( RevPos << 16)
RevPos = ( RevPos >> (32 - Exponent)) & $FFFFFFFF
WL = PeekB(PF1 + i) - $30
PokeL(BufferSinA + 2 * RevPos, WL)
WL = PeekB(PF2 + i) - $30
PokeL(BufferCosA + 2 * RevPos, WL)
FakPos + 1
Next
j = 0
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)
WL = PeekL(BufferCosA + i)
WR = WL
PokeD(BufferFaktor2A + j, WR)
PokeD(BufferFaktor2A + j + 16, WR)
j + 32
Next
TE_Faktoren = ElapsedMilliseconds() - TA_Faktoren
;Winkel erst hier, weil BufferSinA und BufferCosA oben "missbraucht" werden
TA_Winkel = ElapsedMilliseconds()
Rad1 = Pi2 / N ;Konstante
;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
For k = 0 To (N >> 3) ;hier nicht -1 !!!
Rad = Rad1 * k
Si = Sin(Rad)
Co = Cos(Rad)
PokeD(BufferSinA + k << 3, Si) ;0-45°
PokeD(BufferCosA + k << 3, Co) ;0-45°
PokeD(BufferCosA + ((N >> 2) - k) << 3, Si) ;45-90°
PokeD(BufferCosA + ((N >> 2) + k) << 3, -Si) ;90-135°
PokeD(BufferCosA + ((N >> 1) - k) << 3, -Co) ;135-180°
PokeD(BufferSinA + ((N >> 2) - k) << 3, Co) ;45-90°
PokeD(BufferSinA + ((N >> 2) + k) << 3, Co) ;90-135°
PokeD(BufferSinA + ((N >> 1) - k) << 3, Si) ;135-180°
Next
TE_Winkel = ElapsedMilliseconds() - TA_Winkel
;Faktor1
TA_FFT = ElapsedMilliseconds()
Pointer0 = 2
Pointer1 = 0
Pointer2 = Pointer0
While Pointer2 < N
While Pointer2 < N
For k = 1 To Pointer0
EW = (N / Pointer0 << 1) * Pointer1
WR = (PeekD(BufferFaktor1A + 2*Pointer2 * 8) * PeekD(BufferCosA + EW * 8) - PeekD(BufferFaktor1A+8 + 2*Pointer2 * 8) * PeekD(BufferSinA + EW * 8))
WI = (PeekD(BufferFaktor1A+8 + 2*Pointer2 * 8) * PeekD(BufferCosA + EW * 8) + PeekD(BufferFaktor1A + 2*Pointer2 * 8) * PeekD(BufferSinA + EW * 8))
ZR = PeekD(BufferFaktor1A + 2*(Pointer2 - Pointer0) * 8)
ZI = PeekD(BufferFaktor1A+8 + 2*(Pointer2 - Pointer0) * 8)
PokeD(BufferFaktor1A + 2*(Pointer2 - Pointer0) * 8, ZR + WR)
PokeD(BufferFaktor1A+8 + 2*(Pointer2 - Pointer0) * 8, ZI + WI)
PokeD(BufferFaktor1A + 2*Pointer2 * 8, ZR - WR)
PokeD(BufferFaktor1A+8 + 2*Pointer2 * 8, ZI - WI)
Pointer1 + 1
Pointer2 + 1
Next
Pointer1 = 0
Pointer2 + Pointer0
Wend
Pointer0 << 1
Pointer1 = 0
Pointer2 = Pointer0
Wend
Debug "--------------"
;Faktor2
Pointer0 = 2
Pointer1 = 0
Pointer2 = Pointer0
While Pointer2 < N
While Pointer2 < N
For k = 1 To Pointer0
EW = (N / Pointer0 << 1) * Pointer1
WR = (PeekD(BufferFaktor2A + 2*Pointer2 * 8) * PeekD(BufferCosA + EW * 8) - PeekD(BufferFaktor2A+8 + 2*Pointer2 * 8) * PeekD(BufferSinA + EW * 8))
WI = (PeekD(BufferFaktor2A+8 + 2*Pointer2 * 8) * PeekD(BufferCosA + EW * 8) + PeekD(BufferFaktor2A + 2*Pointer2 * 8) * PeekD(BufferSinA + EW * 8))
ZR = PeekD(BufferFaktor2A + 2*(Pointer2 - Pointer0) * 8)
ZI = PeekD(BufferFaktor2A+8 + 2*(Pointer2 - Pointer0) * 8)
PokeD(BufferFaktor2A + 2*(Pointer2 - Pointer0) * 8, ZR + WR)
PokeD(BufferFaktor2A+8 + 2*(Pointer2 - Pointer0) * 8, ZI + WI)
PokeD(BufferFaktor2A + 2*Pointer2 * 8, ZR - WR)
PokeD(BufferFaktor2A+8 + 2*Pointer2 * 8, ZI - WI)
Pointer1 + 1
Pointer2 + 1
Next
Pointer1 = 0
Pointer2 + Pointer0
Wend
Pointer0 << 1
Pointer1 = 0
Pointer2 = Pointer0
Wend
TE_FFT = ElapsedMilliseconds() - TA_FFT
;Und jetzt die beiden transformierten Vektoren miteinander multiplizieren
TA_Multi = ElapsedMilliseconds()
For k = 0 To N - 1
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
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(BufferFaktor1A + 2*k * 8, WR) ;weil vorher BufferFaktor1RA noch benötigt wird
Next
TE_Multi = ElapsedMilliseconds() - TA_Multi
;Die Produkt-Werte wieder revers "verteilen"
TA_Rueck = ElapsedMilliseconds()
FakPos = 0
For i = 0 To N - 1
RevPos = FakPos
RevPos = ((RevPos >> 1) & $55555555) | ((RevPos & $55555555) << 1)
RevPos = ((RevPos >> 2) & $33333333) | ((RevPos & $33333333) << 2)
RevPos = ((RevPos >> 4) & $0F0F0F0F) | ((RevPos & $0F0F0F0F) << 4)
RevPos = ((RevPos >> 8) & $00FF00FF) | ((RevPos & $00FF00FF) << 8)
RevPos = ( RevPos >> 16 ) | ( RevPos << 16)
RevPos = ( RevPos >> (32 - Exponent)) & $FFFFFFFF
PokeD(BufferFaktor2A + 2*8 * RevPos, PeekD(BufferFaktor1A + 2*i * 8))
PokeD(BufferFaktor2A+8 + 2*8 * RevPos, PeekD(BufferFaktor1A+8 + 2*i * 8))
FakPos + 1
Next
;Anfangswerte verteilen für FFT
For i = 0 To N - 2 Step 2
WR = PeekD(BufferFaktor2A + 2*i * 8) + PeekD(BufferFaktor2A + 2*(i + 1) * 8)
WI = PeekD(BufferFaktor2A+8 + 2*i * 8) + PeekD(BufferFaktor2A+8 + 2*(i + 1) * 8)
ZR = PeekD(BufferFaktor2A + 2*i * 8) - PeekD(BufferFaktor2A + 2*(i + 1) * 8)
ZI = PeekD(BufferFaktor2A+8 + 2*i * 8) - PeekD(BufferFaktor2A+8 + 2*(i + 1) * 8)
PokeD(BufferFaktor2A + 2*i * 8, WR)
PokeD(BufferFaktor2A+8 + 2*i * 8, WI)
PokeD(BufferFaktor2A + 2*(i + 1) * 8, ZR)
PokeD(BufferFaktor2A+8 + 2*(i + 1) * 8, ZI)
Next
;nochmal FFT, jetzt aber invers (mit neg.Sinus)
Pointer0 = 2
Pointer1 = 0
Pointer2 = Pointer0
While Pointer2 < N
While Pointer2 < N
For k = 1 To Pointer0
EW = (N / Pointer0 << 1) * Pointer1
WR = (PeekD(BufferFaktor2A + 2*Pointer2 * 8) * PeekD(BufferCosA + EW * 8) + PeekD(BufferFaktor2A+8 + 2*Pointer2 * 8) * PeekD(BufferSinA + EW * 8)) ;neg.Sin, Operand geändert
WI = (PeekD(BufferFaktor2A+8 + 2*Pointer2 * 8) * PeekD(BufferCosA + EW * 8) - PeekD(BufferFaktor2A + 2*Pointer2 * 8) * PeekD(BufferSinA + EW * 8)) ;neg.Sin, Operand geändert
ZR = PeekD(BufferFaktor2A + 2*(Pointer2 - Pointer0) * 8)
ZI = PeekD(BufferFaktor2A+8 + 2*(Pointer2 - Pointer0) * 8)
PokeD(BufferFaktor2A + 2*(Pointer2 - Pointer0) * 8, ZR + WR)
PokeD(BufferFaktor2A+8 + 2*(Pointer2 - Pointer0) * 8, ZI + WI)
PokeD(BufferFaktor2A + 2*Pointer2 * 8, ZR - WR)
PokeD(BufferFaktor2A+8 + 2*Pointer2 * 8, ZI - WI)
Pointer1 + 1
Pointer2 + 1
Next
Pointer1 = 0
Pointer2 + Pointer0
Wend
Pointer0 << 1
Pointer1 = 0
Pointer2 = Pointer0
Wend
TE_Rueck = ElapsedMilliseconds() - TA_Rueck
TA_Ergebnis = ElapsedMilliseconds()
;die Doubles von FFT in Integer konvertieren und Koeffizienten aufaddieren mit Zehner-Versatz
Restore Dezi ;sicher ist sicher!
Pointer0 = 0
For i = 0 To N - 1 ;vom Produkt alle Werte durch
;Konvertierung Double in Integer
QuadWert = Round(PeekD(BufferFaktor2A + Pointer0), #PB_Round_Nearest) ;Double zu Quad
QuadWert >> Exponent ;=div N
; ;Fehler-Abschätzung ;dann Zeile drüber auskommentieren
; !xor rdx,rdx ;sorry, ASM hier der Einfachheit halber
; !mov rcx,[v_N]
; !mov rax,[v_QuadWert]
; !div rcx
; !mov [v_QuadWert],rax
; !or rdx,rdx ;Remainder?
; !jz @f
; !mov [v_Fehler],1 ;Ergebnis anrüchig, Faktoren zu groß
; !@@:
;Hex2Dez
QuadWert1 = QuadWert
Ziffer = 0
Pointer1 = 0
Pointer2 = 0
While PeekQ(?Dezi + Pointer1) <> 0
While (QuadWert - PeekQ(?Dezi + Pointer1)) >= 0
Ziffer + 1
QuadWert - PeekQ(?Dezi + Pointer1)
QuadWert1 = QuadWert
Wend
PokeB(BufferDezA + Pointer2, Ziffer)
Pointer1 + 8
Pointer2 + 1
Ziffer = 0
QuadWert = QuadWert1
Wend
PokeB(BufferDezA + Pointer2, QuadWert1) ;letzte Dezimal-Ziffer abspeichern
Pointer_Produkt1 = Pointer_Produkt
For k = 0 To 15 ;eigentlich hier überdimensioniert
Ziffer = PeekB(BufferProduktA + Pointer_Produkt1)
Ziffer + PeekB(BufferDezA + Pointer2)
PokeB(BufferProduktA + Pointer_Produkt1, Ziffer)
If Ziffer > 9 ;Übertrag
Ziffer - 10
PokeB(BufferProduktA + Pointer_Produkt1, Ziffer)
Ziffer = PeekB(BufferProduktA + Pointer_Produkt1 - 1) + 1
PokeB(BufferProduktA + Pointer_Produkt1 - 1, Ziffer)
EndIf
Pointer_Produkt1 - 1
Pointer2 - 1
Next
Pointer_Produkt - 1
Pointer0 + 16 ;Zeiger in BufferFaktor2A
Next
;------------------------
;Ergebnis (Produkt) in ASCII-Ziffern abspeichern
For i = 0 To N Step 8
PokeQ(BufferProduktA + i, PeekQ(BufferProduktA + i) | $3030303030303030)
Next
PokeB(BufferProduktA + i - 8, 0) ;Zero-Byte zur Sicherheit für String-Auslesen
;Hier könnte man auch sinnvoller Weise vorher Speicher wieder freigeben
Produkt$ = LTrim(PeekS(BufferProduktA), "0") ;evtl. führende Null(en) weg
SetClipboardText(Produkt$) ;oder auskommentieren
TE_Ergebnis = ElapsedMilliseconds() - TA_Ergebnis
TE_Gesamt = ElapsedMilliseconds() - TA_Gesamt
Faktor$ = "Faktoren-Laenge: " + Str(Len(Faktor1$)) + " Bytes" + #CRLF$
Exponent$ = "Exponent 2^ für Faktoren-Laenge: " + Str(Exponent - 1) + #CRLF$
Winkel$ = "Zeit für Sinus/Cosinus-Berechnung: " + Str(TE_Winkel) + " ms" + #CRLF$
Faktoren$ = "Zeit für Faktoren-Aufbereitung: " + Str(TE_Faktoren) + " ms" + #CRLF$
FFT$ = "Zeit für Fast Fourier Transformation: " + Str(TE_FFT) + " ms" + #CRLF$
Multi$ = "Zeit für Multiplikation: " + Str(TE_Multi) + " ms" + #CRLF$
Rueck$ = "Zeit für Rueck-Transformation: " + Str(TE_Rueck) + " ms" + #CRLF$
Ergebnis$ = "Zeit für Ergebnis-Aufbereitung als String: " + Str(TE_Ergebnis) + " ms" + #CRLF$
Gesamt$ = "Zeit gesamt: " + Str(TE_Gesamt) + " ms" + #CRLF$
Fehler$ = "" + #LFCR$ ;oder nachfolgendes
; If Fehler ;wenn Fehler-Test
; Fehler$ = "Ergebnis (möglicherweise) fehlerhaft!" + #LFCR$
; Else
; Fehler$ = "Ergebnis (mit hoher Wahrscheinlichkeit) o.K.!" + #LFCR$
; EndIf
Info$ = Faktor$ + Exponent$ + Winkel$ + Faktoren$ + FFT$ + Multi$ + Rueck$ + Ergebnis$ + Gesamt$ + Fehler$ + #LFCR$
If Len(Produkt$) > 65536 * 2 ;je nach Gusto und was Windows noch sicher anzeigt
MessageRequester("Helles Multiplikation-FFT-Test", Info$)
Else
OpenWindow(0, 0, 0, 1000, 800, "Helles Multiplikation-FFT-Test", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
StringGadget(0, 10, 10, 980, 780, Info$ + Produkt$, #PB_String_ReadOnly | #ES_MULTILINE | #ESB_DISABLE_LEFT | #ESB_DISABLE_RIGHT | #WS_VSCROLL)
Repeat : Until WaitWindowEvent() = #PB_Event_CloseWindow
CloseWindow(0)
EndIf
Z:
DataSection
Dezi:
Data.q 1000000000000000 ;begrenze ich mal auf 16 Stellen; ist damit schon überdimensioniert
Data.q 100000000000000
Data.q 10000000000000
Data.q 1000000000000
Data.q 100000000000
Data.q 10000000000
Data.q 1000000000
Data.q 100000000
Data.q 10000000
Data.q 1000000
Data.q 100000
Data.q 10000
Data.q 1000
Data.q 100
Data.q 10
Data.q 0
EndDataSection
!Helle
P.S.: Natürlich existiert auch eine ASM-Version, aber die ist speziell auf Haswell/Skylake zugeschnitten und teilweise noch in Arbeit.
Edit 3.10.2015: Kleine Optimierungen
Edit 29.12.2015: Speicher-Anforderung optimiert
