4. SSE exp() denkfehler??

SSE exp() denkfehler??

  • Z

    in diesem teil des forums is ja nicht viel los 😃

    also ich habe das problem gelöst.

    --------------
    edit:
    Aber leider nimmt die laufzeit wegen der schleifen für große exponenten linear zu. Damit ist für große und kleine Zahlen die laufzeit leider fürn popo.....
    Während die exp() der math.h quasi eine konstante laufzeit aufweißt,....

    -------------

    Ich werde gleich mal den benchmark plus korregierte version posten, aber ich habe noch ne frage zur abschätzung des fehlers,..

    Wie genau ist exp() ?

    grüße

    0
  • A

    zeusosc schrieb:

    in diesem teil des forums is ja nicht viel los 😃

    Na, dann mach deine Funktionen "portabel", damit man sie unter Linux bauen und damit rumspielen kann 🙂

    0
  • Z

    ... sagen wir mal so,...
    da ich ja jetzt einen guten lösungsansatz für eine konstante laufzeit habe, versuche ich die mathematische und konzeptionelle beschreibung so zu gestalten das jeder der will auch gerne auf ein anderes system (auf x64, avr,m430) portieren kann, wenn es äquivalente mnemonics gibt xD

    0
  • Z

    Happy :xmas1: ,

    also ich hatte mir das ja nochmal angeschaut usw...
    und habe einen evtl. schon bekannten ansatz genommen, ich hoffe ihr versteht
    diesen auf anhieb und könnt mich dann auch sofort korregieren 😃

    Zuerst schreibe ich die Aufgabe ein bissl um:

    ex=y⇔y=2x⋅Log2(e)e^{x}=y \Leftrightarrow y= 2^{x\cdot \text{Log}_2(e)}

    Da der VS2008 Compiler ein bissl doof ist,(der schneidet bei der 7ten
    nachkommastelle gerne etwas ab) habe ich die double wert für

    \text{Log}_2(e)$$ per hand ausgerechnet: ```cpp static const __int64 l2de =0x3FF71547652B82FE; //1,4426950408889634 ``` Nun kommt trick17: ich zerlege den exponenten $$ x\cdot\text{Log}_2(e)

    in zwei Komponenten, einmal den nachkommaanteil und den Ganzzahlanteil:

    2x⋅Log_2(e)=2F_L+F_H=2F_L⋅2FH2^{x\cdot \text{Log}\_2(e)}=2^{F\_L + F\_H}=2^{F\_L}\cdot 2^{F_H}

    Dabei bezeichnet $$F_L,F_H$$ den "FracLow" und "FracHigh" des
    exponenten.

    Dann haben die Faktoren folgende Eigenschaften:
    A) $$2^{F_L} \in \left[1,2\right] \Leftrightarrow F_L \in \left[0,1\right]$$
    und für die double schreibweise:

    ⇒(1+MF_L252)⋅2E_FL\Rightarrow \left(1+\frac{M_{F\_L}}{2^{52}}\right)\cdot 2^{E\_{F_L}}

    😎 $$2^{F_H}\in \mathbb{N}_\pm \Leftrightarrow F_H \in \mathbb{N}_\pm$$
    und für die double schreibweise:

    2F_H∈N_±⇒MF_H=0⇒(1+0252)⋅2E_F_H⋅(±1_FH)2^{F\_H} \in \mathbb{N}\_\pm \Rightarrow M_{F\_H}=0 \Rightarrow \left(1+\frac{0}{2^{52}}\right)\cdot 2^{E\_{F\_H}}\cdot (\pm 1\_{F_H})

    Die Folgerung währe:

    y=2F_L⋅2F_H⇒(1+MF_L252)⋅(±1_F_H)⋅2E_F_L+E_FHy=2^{F\_L}\cdot 2^{F\_H} \Rightarrow \left(1+\frac{M_{F\_L}}{2^{52}}\right)\cdot (\pm 1\_{F\_H}) \cdot 2^{E\_{F\_L} + E\_{F_H}}

    Würde man sich wohl mit 2-4 stellen genauigkeit nach dem komma genügen,
    so würde es reichen den Fraclow als mantisse zu nehmen
    (lineare regression von $$ \sqrt[n]{2}$$), ansonsten halt ne
    Reihenentwicklung, so wie ich es gemacht habe....

    Ich habe vergessen: $$ (\pm 1_{F_H}) \rightarrow 1 \wedge {E_{F_H}}\in \mathbb{N}_\pm$$
    das vorzeichen von dem exponenten fließt ja in den Ganzkommaanteil

    F_H$$ und somit in $$E_{F_H}$$ ein.... Comments are welcome 😉 :xmas2:
    0
  • Z

    Soo, ich denke ich bin fertig, mit ner akzeptablen abstand zur exp funktion der math.h...

    Hier einmal die Ausgabe der testdurchläufe:
    (//Compiled with /arch:SSE2)

    Cycles:                 1e+009
Runtime at all: 80073.3
Runtime per Cycle: 8.00733e-005

-------          exp()   -------
Cycles:                 1e+009
Runtime:                        39642.5
Runtime per Cycle:      3.96425e-005
Result:                 6.05732e+021

-------          _expf()         -------
Cycles:                 1e+009
Runtime:                        37629.8
Runtime per Cycle:      3.76298e-005
Result:                 6.05732e+021

-----------------------------------------------------
Absolute Timedifference: 2012.65
Arithmated middle mantissa distance between exp() and _expf():  0.000355568
Please press enter....

    Hier die Testfunktionen (in c/c++)

    //Compiled with /arch:SSE2

std::string dts(double x)
{
	std::stringstream strRet;
	strRet<<x;
	return strRet.str();
};

std::string _mark(std::string strName,double e,double t, double div, double cycles)
{
	std::string strRet="------- \t "+strName+"\t -------\r\n";
	strRet+="Cycles:			"+dts(cycles)+"\r\n";
	strRet+="Runtime:			"+dts(t)+"\r\n";
	strRet+="Runtime per Cycle:	"+dts(t*div)+"\r\n";
	strRet+="Result:			"+dts(e)+"\r\n";
	return strRet;
};

std::string _print_bench_result(double er1, double er2, double t1, double t2,double loops_div, double cycles,double ts)
{

	std::string strRet;
	strRet+="Cycles:			"+dts(cycles)+"\r\n";
	strRet+="Runtime at all:	"+dts(ts)+"\r\n";
	strRet+="Runtime per Cycle: "+dts(ts*loops_div)+"\r\n\r\n";
	strRet+=_mark("exp()",er1,t1,loops_div,cycles)+"\r\n";
	strRet+=_mark("_expf()",er2,t2,loops_div,cycles)+"\r\n\r\n";
	strRet+="-----------------------------------------------------\r\n";
	strRet+="Absolute Timedifference: "+dts(fabs(t1-t2))+"\r\n";
	return strRet;
};

int _tmain(int argc, _TCHAR* argv[])
	{
_stde::_cpu_ticks_struct cpu,cpu1; //winapi QueryPerformanceCounter
cpu.GetMilliseconds();

static double ergebniss1=0;
static double ergebniss2=0;
static double x= 0.155555555555f;

static const unsigned int end=1000000;
static const unsigned int end2=1000;
static const double cycles=(((double)end)*((double) end2));
static const double divi=1.0f/cycles;
unsigned int i=0;
unsigned int j=0;
double time_exp=0.0f;
double time_expf_test=0.0f;
double diff=0;
static double adder=50.0f/end;
cpu_1.GetMilliseconds();
while(j<end)
{
	adder=adder*(-1.0f);
	x=x*(-1.0f);
	x=x+adder;

	i=0;
	cpu.GetMilliseconds();
	while(i<end2)
	{
			ergebniss1=exp(x);
		i++;
	}
	time_exp+=cpu.GetMilliseconds();

	i=0;
	cpu.GetMilliseconds();
	while(i<end2)
	{
			ergebniss2=_expf(x);
		i++;
	};
	time_expf_test+=cpu.GetMilliseconds();
int n=0;
double blub=frexp(ergebniss2,&n);
diff+= (ergebniss2-ergebniss1)/pow(2.0f,n);

j++;
};
diff=diff/cycles;
std::cout<<_print_bench_result(ergebniss1,ergebniss2,time_exp,time_expf_test,divi,cycles,cpu_1.GetMilliseconds());
std::cout<<"Arithmated middle mantissa distance between exp() and _expf():\t"<<diff<<"\r\nPlease press enter....";
std::cin>>j;

return 0;
};

    sooo und nun mein hilight

    inline double __fastcall _expf(double y)
{

	double debug_temp=0.0f;
__int64 blub=1;
int a1=0;

	__asm{

		movsd	xmm1, y				// y to xmm1
		movsd   xmm0, l2de		// log2(e) at xmm0
		//multiply it
		mulsd	xmm0,xmm1			//build the exponent +-y*log2(e)
		//convert
		cvttsd2si eax,xmm0
		cvtsi2sd  xmm2,eax

//movq blub,xmm1
//movsd debug_temp,xmm1
//movq blub,xmm0	
//movsd debug_temp,xmm0
//movq blub,xmm2
//movsd debug_temp,xmm2
		//now we splitt into frachigh(xmm1) and frac low part(xmm0)
		subsd	xmm0,xmm2
		//now we take our new polynom

		movsd xmm3,_pb6
		mulsd xmm3,xmm0

		addsd xmm3,_pb5
		mulsd xmm3,xmm0

		addsd xmm3,_pb4
		mulsd xmm3,xmm0

		addsd xmm3,_pb3
		mulsd xmm3,xmm0

		addsd xmm3,_pb2
		mulsd xmm3,xmm0

		addsd xmm3,_pb1
		mulsd xmm3,xmm0

		addsd xmm3,_pb0

//##############################################################
//xmm3,xmm2,xmm0,eax are reserved
//here we convert the frac high part to the exponentialform
		//movsd		xmm1,SignPosd

		//store the sign
		movsd		xmm1,xmm2

		psrlq		xmm1,63 //shift sign to right
		movq		blub,xmm1
		mov			ebx,dword ptr[blub]
		//test		ebx,ebx
		//jnz					

		//now we extract the exponent of xmm3
		movsd		xmm2,ExponentPosd
		andpd		xmm2,xmm3
		psrlq		xmm2,52
		movsd		blub,xmm2
		mov			ecx,dword ptr[blub]
		sub			ecx,0x3fe
	//	int 3
		//proof if sgn was set
		test		ebx,ebx
		jz			SHORT _no_sgn
		jmp			SHORT _sgn

_no_sgn:
		add		ecx,eax	//just add the exponents
		jmp SHORT _fin_sum_exponents
_sgn:
		sub		ecx,eax //just sub the exponents
_fin_sum_exponents:

		add ecx,0x3fe //we have to add the zeropoint again
		and	ecx,0x7ff //we have to make sure that only 7ff can be written

		//move back
		mov dword ptr[blub],ecx

		movsd	xmm2,invExponentPosd
		movsd	xmm1,blub
		andpd	xmm2,xmm3
		movsd   blub,xmm2
//int 3
		psllq	xmm1,52 //shift back left
		movsd   blub,xmm1
//int 3
		orpd	xmm2,xmm1 //and set the exponent

//##############################################################

		movsd y,xmm2	

	}

return y;
};

    Also wenn ihr noch Optimierungspotenzial seht, bitte uuuuunbedingt bescheid sagen....
    ----------------------------------------------------------------------------
    Edit:
    hier noch ein paar präzise konstanten:

    static const __int64 l2de	=0x3FF71547652B82FE; //1.4426950408889634
static const __int64 _pb0	=0x3FF0000000000000;//1.000000f
static const __int64 _pb1	=0x3FE62E429E0A41A2;//0.693147;
static const __int64 _pb2	=0x3FCEBFF47735C183;//0.240233;
static const __int64 _pb3	=0x3FAC6653F1EF9C05;//0.0554682f;
static const __int64 _pb4	=0x3F83E27EE5B92CDD;//0.00970935
static const __int64 _pb5	=0x3F53EBF0289C8156;//0.00121592
static const __int64 _pb6	=0x3F2DB52DDAEAD8E3;//0.000226652
static const __int64 _fhDiv	=0x10000000000000;// 2^52
/*1. + 0.693147 x + 0.240233 x^2 + 0.0554682 x^3 + 0.00970935 x^4 + 
 0.00121592 x^5 + 0.000226652 x^6*/

    Und die BitMasken:

    static const unsigned __int64 SignPosd			=0x8000000000000000;
static const unsigned __int64 ExponentPosd		=0x7ff0000000000000;
static const unsigned __int64 invExponentPosd	=0x800fffffffffffff;
static const unsigned __int64 MantissaPosd		=0x000FFFFFFFFFFFFF;
    0
  • ?

    hi,
    das ist gut! - habe es noch etwas verkürzt (masm-syntax). Könntest es ja mal mit deiner Version vergleichen (und hier posten ;-)). Wäre auch toll, wenn du mal die Formel für die Polynom-Koeffizienten zeigst.

    ; return through xmm0
expf proc x:REAL8

.const
	align 16
	l2de	dq 3FF71547652B82FEh	; 1.4426950408889634 
	_pb6	dq 3F2DB52DDAEAD8E3h	; 0.000226652 
	_pb5	dq 3F53EBF0289C8156h	; 0.00121592 
	_pb4	dq 3F83E27EE5B92CDDh	; 0.00970935
	_pb3	dq 3FAC6653F1EF9C05h	; 0.0554682f
	_pb2	dq 3FCEBFF47735C183h	; 0.240233
	_pb1	dq 3FE62E429E0A41A2h	; 0.693147
	_pb0	dq 3FF0000000000000h	; 1.000000f 
.code

	movsd xmm1,l2de        ; log2(e) at xmm0 
	mulsd xmm1,x           ; build the exponent +-y*log2(e) 

	cvttsd2si eax,xmm1 
	cvtsi2sd  xmm0,eax 

	lea eax,[eax*8]		;
	lea eax,[eax+eax]	; shl eax,4

	subsd xmm1,xmm0 

	movsd xmm0,_pb6 
	mulsd xmm0,xmm1 

	addsd xmm0,_pb5 
	mulsd xmm0,xmm1 

	addsd xmm0,_pb4 
	mulsd xmm0,xmm1 

	addsd xmm0,_pb3 
	mulsd xmm0,xmm1 

	addsd xmm0,_pb2 
	mulsd xmm0,xmm1 

	addsd xmm0,_pb1 
	mulsd xmm0,xmm1 

	addsd xmm0,_pb0 

	pextrw ecx,xmm0,3
	mov edx,ecx
	and edx,0800fh
	jns @1
	neg eax	
@1: 
	and ecx,07ff0h
	lea ecx,[ecx+eax]
	and ecx,07ff0h
	lea ecx,[ecx+edx]
	pinsrw xmm0,ecx,3
	; xmm0 holds result

	ret

expf endp

    grüße, masm

    0
  • Z

    Jo ich habe die verkürzte fassung noch nicht ausprobiert (mache ich aber heute noch)...

    Zu den Koeffizienten,..
    Das Polynom ist wie oben beschrieben die entwicklung von

    \sqrt[\hat{x}]{2}=2^{\frac{1}{\hat{x}}}=2^{x}\forall x\in \left[0,..,1\right]\subset\mathbb{R}

    Diese erhält man entweder durch eine Taylorapproximation um den Entwicklungspunkt

    x_0=0$$ oder durch einen Polynomialen fit in Origin oder Mathematica für die genannte Funktion.... Um den Punkt 0.5 hat man zu große abstände an den Rändern 1,2 von der funktion,... greetz \-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\- Edit: ich hatte schon über eine lookup table nachgedacht, aber 2^52 werte reinklatschen,... dann noch die 64 bit addressierung...... und die ganze rumrechnerei erstmal :D, Da müsste ich n makro in Mathematica schreiben das mir die werte exakt umrechnet...(ok, letztes bit ungenauigkeit ist ja nach ieee 7x erlaubt....) 🤡
    0
  • Z

    geilomat masm,...

    deine modifikation ist ja mal richtig performant 😉

    Cycles:                 1e+009
Runtime at all: 70561.2
Runtime per Cycle: 7.05612e-005

-------          _expf_masm()   -------
Cycles:                 1e+009
Runtime:                        30499.4
Runtime per Cycle:      3.04994e-005
Result:                 6.05732e+021

-------          _expf()         -------
Cycles:                 1e+009
Runtime:                        37343.9
Runtime per Cycle:      3.73439e-005
Result:                 6.05732e+021

-----------------------------------------------------
Absolute Timedifference: 6844.58
Arithmated middle mantissa distance between _expf_masm() and _expf():  0.000355567
Please press enter....

    leider weiß ich nicht woher die differenz der mantisse kommen soll, denn vorher habe ich drei zufallszahlen überprüft und die waren identisch bis zur letzten stelle....

    Ich habe noch eine idee wie man evtl. eine höhere Präzision und geringere Laufzeit rausholt,... da würde die berechnung des fraclow von 13 auf 7 anweisungen schrumpfen...

    greetz n thx
    ---------------------------------------------------------------------------
    Edit:

    Also ich habe den maximalen Absoluten fehler der 2^fraclow komponente durch
    neue koeffizienten von rund $$|1.2E-7|$$
    auf gerundet $$|5.5463E-9|$$ reduzieren können...

    p(x)=1. + 0.693147 x + 0.240231 x^2 + 0.0554788 x^3 + 0.00968722 x^4 +
    0.00123681 x^5 + 0.000219337 x^6,..

    Mit $$2^x-p(x)$$ kann man die differenzfunktion plotten die die abweichung zeigt...

    0
  • N

    entschuldigung, ich habe da nur ein paar Verständnis Fragen:

    wie war nochmal der eigentliche theoretisch/praktische Hintergrund?
    welche Funktion, bzw. welcher Algo soll genau beschleunigt werden? bzw. welcher Algo ist wo genau zu langsam? (nicht der mathematische, der ist sowieso immer zu langsam 😉 , der in welchem asm-code).

    Also z.B. ist der (zu langsame) Ausgangspunkt ein asm-Programm wie:

    ;standardrundungsverhalten der Fpu
fld exponent
fldl2e
fmulp
fld st(0)
frndint                 ;falls nötig, Control-Word anpassen
fsub st(1), st(0)
fxch
f2xm1
fld1
faddp st(1), st(0)
fscale
fstp st(1)
;Ergebnis in st(0)

    ?

    0
  • Z

    Hi nachtfeuer,

    also wenn nicht mit der option /arch:SSE2 /arch:SSE unter VS2008 kompiliert wird,
    wird die FPU für die berechnung der exponentialfunktion (math.h) genutzt. Wird hingegen mit dieser Compilerflag kompiliert, wird die Kalkulation der exponentialfunktion mittels SSE vorgenommen, was eine beschleunigung bei rund 1e9 durchgängen von rund 7sec bedeutet.

    Da ich die berechnung nun etwas weiter beschleunigen wollte , habe ich mir überlegt wie die berechnung sein müsste und versucht diese per hand zu optimieren.
    (Danke an masm, der nochmal aus meinen source riiiichtig performance rausgeholt hat 😃 )

    Nun, die methode wie die exponentialfunktion berechnet wird habe ich gepostet, wie groß die abweichung der fraclow komponente ist habe ich gepostet.

    Nun habe ich eine funktion die bei 1e9 durchgängen gegenüber der fpu variante rund 17 sec spart,... und das bei einer genauigkeit von rund 8 stellen der Mantisse...

    greetz
    --------------------------------
    Edit:
    Benchmark @ Intel Core2Duo P8400 @ 2.66Ghz -VS2008 Pro-DebugMode

    0
Anmelden zum Antworten