Schnelle Cosinus-Berechnung

Hier könnt Ihr gute, von Euch geschriebene Codes posten. Sie müssen auf jeden Fall funktionieren und sollten möglichst effizient, elegant und beispielhaft oder einfach nur cool sein.
Antworten
Benutzeravatar
Helle
Beiträge: 566
Registriert: 11.11.2004 16:13
Wohnort: Magdeburg

Schnelle Cosinus-Berechnung

Beitrag 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
Zuletzt geändert von Helle am 28.01.2007 14:30, insgesamt 1-mal geändert.
Nach oben
Benutzeravatar
Helle
Beiträge: 566
Registriert: 11.11.2004 16:13
Wohnort: Magdeburg

Beitrag 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
Nach oben
Antworten

Zurück zu „Code, Tipps und Tricks“