PureBoard

Das deutsche PureBasic-Forum
https://www.purebasic.fr/german/

Schnelle Cosinus-Berechnung

https://www.purebasic.fr/german/viewtopic.php?t=11796

Seite 1 von 1

Schnelle Cosinus-Berechnung

Verfasst: 27.01.2007 21:39
von Helle
Wer eine CPU mit SSE2-Unterstützung im PC hat und dem der Wertebereich von -90° bis +90° ausreicht kann Folgendes testen:

Code: Alles auswählen

;- Cosinus mittels Funktionsreihe und SSE2 (Double-Precision) berechnen
;- Empfohlener Wertebereich: -90° bis +90° (Rad = -1.57079632679 bis +1.57079632679)
;- "Helle" Klaus Helbing, 27.01.2007, PB v4.02
;- Cos(x) = 1 - (x^2)/2! + (x^4)/4! - (x^6)/6! + (x^8)/8! - (x^10)/10! + (x^12)/12! ... hier nicht mehr Glieder
;- Die Grundidee stammt aus dem Flatassembler-Forum

Global Rad.d = 0.987654321   ;sinnvollen Wertebereich s.o.
Global CosPB.d
Global CosSSE.d

;- Benötigt CPU mit SSE2-Unterstützung -------------------- 
Macro CosSSE2(a,b)
     LEA edi,a               ;hier Speicher-Adresse von Rad, Rad = x 
     !movsd xmm0,[edi]       ;1 = x    1 = unteres Quad-Word  2 = oberes Quad-Word
     !movhpd xmm0,[edi]      ;1 und 2 = x
     !lea esi,[v_RezFak]     ;Anfangs-Speicher-Adresse der reziproken Fakultäten  
     !movapd xmm1,xmm0       ;XMM1 = XMM0 
     !mulpd xmm0,xmm1        ;XMM0: 1 = x^2   2 = x^2 
     !movapd xmm1,xmm0       ;XMM1: 1 = x^2   2 = x^2
     
     !mulsd xmm0,xmm1        ;XMM0: 1 = x^4   2 = x^2 
     !mulpd xmm1,xmm1        ;XMM1: 1 = x^4   2 = x^4 
     !movupd xmm3,[esi]      ;Anfangsglied   XMM3: 1 = 1.0  2 = 0.0 
     !movapd xmm2,xmm0       ;XMM2 = XMM0 
     !mov ecx,-16*3          ;3*2 weitere Glieder, also 6 Iterationen (immer paarweise, daher 3 Schleifendurchläufe) 

!@@:                         ;Test-Werte für 1. und 2. Iteration
     !movupd xmm4,[esi+64+ecx]  ;1 = 0.04166...   2 = -0.5        ;64=(16*3)+16
     !mulpd xmm0,xmm4        ;1 = x^4*0.04166...  2 = x^2*-0.5 
     !addpd xmm3,xmm0        ;1 = 1.0+0.04166...  2 = 0.0-0.5 
     !mulpd xmm2,xmm1        ;XMM2: 1 = x^8   2 = x^6   ;nach Ende-Test brachte nichts
     !movapd xmm0,xmm2       ;XMM0 = XMM2 
     !add ecx,16
     !jnz @r    
    
     LEA edi,b               ;Rückgabewert-Variable
     !movhlps xmm4,xmm3      ;Ergebnis ist die Summe der beiden Quad-Words von XMM3
     !addsd xmm4,xmm3   
     !movsd qword[edi],xmm4 
EndMacro
;----------------------------------------------------------

;- Struktur wegen Speicherverwaltung ----------------------
Structure Cosinus 
  RF1.d
  RF2.d
  RF3.d
  RF4.d
  RF5.d 
  RF6.d   
  RF7.d 
  RF8.d 
EndStructure

Define.Cosinus RezFak
RezFak\RF1  =  1.0 
RezFak\RF2  =  0.0 
RezFak\RF3  =  0.041666666666666667    ; 1/4!
RezFak\RF4  = -0.5                     ;-1/2!
RezFak\RF5  =  2.480158730158730159e-5 ; 1/8!
RezFak\RF6  = -0.001388888888888889    ;-1/6!
RezFak\RF7  =  2.087675698786809898e-9 ; 1/12!
RezFak\RF8  = -2.755731922398589065e-7 ;-1/10!
;----------------------------------------------------------

;- Test SSE2
T1= ElapsedMilliseconds() 
For i = 0 To 99999999 
  CosSSE2(Rad, CosSSE)                 ;Variablen freier Wahl
Next 
T1 = ElapsedMilliseconds() - T1 

;- Test PB
T2= ElapsedMilliseconds() 
For i = 0 To 99999999 
  CosPB = Cos(Rad) 
Next 
T2 = ElapsedMilliseconds() - T2 

Cos$ = StrD(CosSSE) + " - SSE2 - " + Str(T1) + " ms" + #CRLF$ 
Cos$ + StrD(CosPB) + " - PB -     " + Str(T2) + " ms" + #CRLF$ 
MessageRequester("Cosinus Double-Precision", Cos$) 

End
Gruss
Helle

Edit: Sprungadresse anonymisiert

Verfasst: 28.01.2007 12:29
von Helle
Hier noch Versionen für Sinus und Tangens:

Code: Alles auswählen

;- Sinus mittels Funktionsreihe und SSE2 (Double-Precision) berechnen
;- Empfohlener Wertebereich: -90° bis +90° (Rad = -1.57079632679 bis +1.57079632679)
;- "Helle" Klaus Helbing, 27.01.2007, PB v4.02
;- Sin(x) = x - (x^3)/3! + (x^5)/5! - (x^7)/7! + (x^9)/9! - (x^11)/11! + (x^13)/13!... hier nicht mehr Glieder

Global Rad.d = 0.123456789   ;sinnvollen Wertebereich s.o.
Global SinPB.d
Global SinSSE.d

;- Benötigt CPU mit SSE2-Unterstützung -------------------- 
Macro SinSSE2(a,b)
     LEA edi,a               ;hier Speicher-Adresse von Rad, Rad = x 
     !movsd xmm0,[edi]       ;1 = x    1 = unteres Quad-Word  2 = oberes Quad-Word
     !movhpd xmm0,[edi]      ;1 und 2 = x
     !lea esi,[v_RezFak]     ;Anfangs-Speicher-Adresse der reziproken Fakultäten  
     !movapd xmm1,xmm0       ;XMM1 = XMM0 
     !movapd xmm5,xmm0       ;Zwischenspeicherung
     !mulpd xmm0,xmm1        ;XMM0: 1 = x^2   2 = x^2 
     !movapd xmm1,xmm0       ;XMM1: 1 = x^2   2 = x^2
     !mulsd xmm0,xmm1        ;XMM0: 1 = x^4   2 = x^2 
     !mulpd xmm1,xmm1        ;XMM1: 1 = x^4   2 = x^4 
     !mulpd xmm0,xmm5        ;XMM0: 1 = x^5   2 = x^3
     !movq xmm3,[edi]        ;Anfangsglied   XMM3: 1 = x   2 = 0.0 
     !movapd xmm2,xmm0       ;XMM2 = XMM0 
     !mov ecx,-16*3          ;3*2 weitere Glieder, also 6 Iterationen (immer paarweise, daher 3 Schleifendurchläufe) 

!@@:                         ;Test-Werte für 1. und 2. Iteration
     !movupd xmm4,[esi+48+ecx]  ;1 = 0.08333...   2 = -0.1666...    ;48=16*3
     !mulpd xmm0,xmm4        ;1 = x^5*0.08333...  2 = x^3*-0.1666... 
     !addpd xmm3,xmm0        ;1 = x+0.08333...    2 = 0.0-0.1666... 
     !mulpd xmm2,xmm1        ;XMM2: 1 = x^9   2 = x^7   ;nach Ende-Test brachte nichts
     !movapd xmm0,xmm2       ;XMM0 = XMM2 
     !add ecx,16
     !jnz @r    
    
     LEA edi,b               ;Rückgabewert-Variable
     !movhlps xmm4,xmm3      ;Ergebnis ist die Summe der beiden Quad-Words von XMM3
     !addsd xmm4,xmm3   
     !movsd qword[edi],xmm4 
EndMacro
;----------------------------------------------------------

;- Struktur wegen Speicherverwaltung ----------------------
Structure Sinus 
  RF1.d
  RF2.d
  RF3.d 
  RF4.d   
  RF5.d 
  RF6.d 
EndStructure

Define.Sinus RezFak
RezFak\RF1  =  0.008333333333333333    ; 1/5!
RezFak\RF2  = -0.166666666666666667    ;-1/3!
RezFak\RF3  =  2.755731922398589065e-6 ; 1/9!
RezFak\RF4  = -0.000198412698412698    ;-1/7!
RezFak\RF5  =  1.605904383682161460e-10; 1/13!
RezFak\RF6  = -2.505210838544171878e-8 ;-1/11!
;----------------------------------------------------------

;- Test SSE2
T1= ElapsedMilliseconds() 
For i = 0 To 99999999 
  SinSSE2(Rad, SinSSE)                 ;Variablen freier Wahl
Next 
T1 = ElapsedMilliseconds() - T1 

;- Test PB
T2= ElapsedMilliseconds() 
For i = 0 To 99999999 
  SinPB = Sin(Rad) 
Next 
T2 = ElapsedMilliseconds() - T2 

Sin$ = StrD(SinSSE) + " - SSE2 - " + Str(T1) + " ms" + #CRLF$ 
Sin$ + StrD(SinPB) + " - PB -     " + Str(T2) + " ms" + #CRLF$ 
MessageRequester("Sinus Double-Precision", Sin$) 

End

Code: Alles auswählen

;- Tangens mittels Funktionsreihe und SSE2 (Double-Precision) berechnen
;- Empfohlener Wertebereich: -45° bis +45° (Rad = -0.78539816339744831 bis +0.78539816339744831)
;- "Helle" Klaus Helbing, 28.01.2007, PB v4.02
;- Tan(x) = x + F3*(x^3) + F5*(x^5) + F7*(x^7) + F9*(x^9) + F11*(x^11) + F13*(x^13) + F15*(x^15) + F17*(x^17)... hier nicht mehr Glieder

Global Rad.d = 0.456789123   ;sinnvollen Wertebereich s.o.
Global TanPB.d
Global TanSSE.d

;- Benötigt CPU mit SSE2-Unterstützung -------------------- 
Macro TanSSE2(a,b)
     LEA edi,a               ;hier Speicher-Adresse von Rad, Rad = x 
     !movsd xmm0,[edi]       ;1 = x    1 = unteres Quad-Word  2 = oberes Quad-Word
     !movhpd xmm0,[edi]      ;1 und 2 = x
     !lea esi,[v_Faktor]     ;Anfangs-Speicher-Adresse der Faktoren  
     !movapd xmm1,xmm0       ;XMM1 = XMM0 
     !movapd xmm5,xmm0       ;Zwischenspeicherung
     !mulpd xmm0,xmm1        ;XMM0: 1 = x^2   2 = x^2 
     !movapd xmm1,xmm0       ;XMM1: 1 = x^2   2 = x^2
     !mulsd xmm0,xmm1        ;XMM0: 1 = x^4   2 = x^2 
     !mulpd xmm1,xmm1        ;XMM1: 1 = x^4   2 = x^4 
     !mulpd xmm0,xmm5        ;XMM0: 1 = x^5   2 = x^3
     !movq xmm3,[edi]        ;Anfangsglied   XMM3: 1 = x   2 = 0.0 
     !movapd xmm2,xmm0       ;XMM2 = XMM0 
     !mov ecx,-16*4          ;4*2 weitere Glieder, also 8 Iterationen (immer paarweise, daher 4 Schleifendurchläufe) 

!@@:                         ;Test-Werte für 1. und 2. Iteration
     !movupd xmm4,[esi+64+ecx]  ;1 = 0.13333...   2 = 0.33333...        ;64=16*4
     !mulpd xmm0,xmm4        ;1 = x^5*0.13333...  2 = x^3*0.33333... 
     !addpd xmm3,xmm0        ;1 = x+0.13333...    2 = 0.0+0.33333... 
     !mulpd xmm2,xmm1        ;XMM2: 1 = x^9   2 = x^7   ;nach Ende-Test brachte nichts
     !movapd xmm0,xmm2       ;XMM0 = XMM2 
     !add ecx,16
     !jnz @r    
    
     LEA edi,b               ;Rückgabewert-Variable
     !movhlps xmm4,xmm3      ;Ergebnis ist die Summe der beiden Quad-Words von XMM3
     !addsd xmm4,xmm3   
     !movsd qword[edi],xmm4 
EndMacro
;----------------------------------------------------------

;- Struktur wegen Speicherverwaltung ----------------------
Structure Tangens 
  F5.d
  F3.d
  F9.d 
  F7.d   
  F13.d 
  F11.d 
  F17.d
  F15.d
EndStructure

Define.Tangens Faktor                  ;Bernoulli lässt grüssen!
Faktor\F5  =  0.133333333333333333
Faktor\F3  =  0.333333333333333333
Faktor\F9  =  0.021869488536155206
Faktor\F7  =  0.053968253968253971
Faktor\F13 =  0.003592128036572481
Faktor\F11 =  0.008863235529902195
Faktor\F17 =  0.000590027440945586
Faktor\F15 =  0.001455834387051318
;----------------------------------------------------------

;- Test SSE2
T1= ElapsedMilliseconds() 
For i = 0 To 99999999 
  TanSSE2(Rad, TanSSE)                 ;Variablen freier Wahl
Next 
T1 = ElapsedMilliseconds() - T1 

;- Test PB
T2= ElapsedMilliseconds() 
For i = 0 To 99999999 
  TanPB = Tan(Rad) 
Next 
T2 = ElapsedMilliseconds() - T2 

Tan$ = StrD(TanSSE) + " - SSE2 - " + Str(T1) + " ms" + #CRLF$ 
Tan$ + StrD(TanPB) + " - PB -     " + Str(T2) + " ms" + #CRLF$ 
MessageRequester("Tangens Double-Precision", Tan$) 

End
Wenn geringere Genauigkeit reicht kann man auch die Anzahl der Iterationen reduzieren.

Gruss
Helle

Edit: Sprungadresse anonymisiert
Alle Zeiten sind UTC+02:00
Seite 1 von 1

Powered by phpBB® Forum Software © phpBB Limited

Deutsche Übersetzung durch phpBB.de