sinus funktion vereinfachen bzw. schneller berechnen lassen



  • dein code ist doof. viel zu renundant. eine schleife, die bis 6 zählt, wäre doch echt erlaubt gewesen.

    ich hacke den mal rein

    #include <cmath>
    #include <iostream>
    using namespace std;
    
    double z0[40000];
    double const PI=3.1415926535897932384626433832795;
    double const PI1=PI*2;
    double const PI2=PI*2;
    double const PI3=PI*2;
    double const PI4=PI*2;
    double const PI5=PI*2;
    double const PI6=PI*2;
    double C=300000;
    int a1=100,a2=150;
    int b1=200,b2=250;
    int c1=300,c2=350;
    int d1=400,d2=450;
    int e1=500,e2=550;
    int f1=600,f2=650;
    double const PI_1=PI;
    double const PI_2=PI;
    double const PI_3=PI;
    double const PI_4=PI;
    double const PI_5=PI;
    double const PI_6=PI;
    int const Anzahl=6;
    
    double LA=10;
    double LB=10;
    double LC=10;
    double LD=10;
    double LE=10;
    double LF=10;
    
    inline double sqrt2(double x){
    	return sqrt(x);
    }
    
    inline double sin2(double x){
    	return sin(x);
    }
    
    void test(int s,double t){
    	for(int y=0;y<200;y++)        // 200 Zeilen 
    	{ 
    		for(int x=0;x<200;x++)    // 200 Spalten 
    		{ 
    			z0[s]=floor(((    
    				sin ( PI1 * ( t * C - sqrt2 ( (x - a1) * (x - a1) + (y - a2) * (y - a2) ) / LA ) + PI_1 ) 
    				+    sin ( PI2 * ( t * C - sqrt2 ( (x - b1) * (x - b1) + (y - b2) * (y - b2) ) / LB ) + PI_2 ) 
    				+    sin ( PI3 * ( t * C - sqrt2 ( (x - c1) * (x - c1) + (y - c2) * (y - c2) ) / LC ) + PI_3 ) 
    				+    sin ( PI4 * ( t * C - sqrt2 ( (x - d1) * (x - d1) + (y - d2) * (y - d2) ) / LD ) + PI_4 ) 
    				+    sin ( PI5 * ( t * C - sqrt2 ( (x - e1) * (x - e1) + (y - e2) * (y - e2) ) / LE ) + PI_5 ) 
    				+    sin ( PI6 * ( t * C - sqrt2 ( (x - f1) * (x - f1) + (y - f2) * (y - f2) ) / LF ) + PI_6 )     
    					) * 255/Anzahl) + 0.5 );  
    		}
    	}
    }
    
    typedef unsigned long long u64;
    
    u64 rdtsc() {
      long long tmp;
      asm volatile("rdtsc" : "=A" (tmp));
      return tmp;
    }
    
    int main(){
    	u64 time=-rdtsc();
    	for(int s=0;s<10;++s){
    		test(s,1.0*s/100);
    	}
    	time+=rdtsc();
    	cout<<time/400000000.0<<endl;
    	return 0;
    }
    

    und wie du siehst, werden die messergebnisse unzuverlässig. weil es mir viel zu aufwendig ist bei dieser schneittstelle ein anständiges umfeld zu simulieren.

    naja, messwert

    3.97

    mal gucken, was gregors tabelle bringt.

    inline double sqrt2(double x){
    	return sqrt(x);
    }
    
    double* initSinTab(){
    	static double sinTab[1024];
    	for(int i=0;i<1024;++i)
    		sinTab[i]=sin(2*PI/i);
    	return sinTab;
    }
    double* sinTab=initSinTab();
    
    inline double sin2(double x){
    //	return sin(x);
    	if(x<0) return -sin2(-x);
    	x=fmod(x,2*PI);
    	return sinTab[int(x/(2*PI)*1024)];
    }
    

    3.37

    naja, fast nix.

    mal sehe, was _hypoz zu bieten hat

    z0[s]=floor(((    
    				sin ( PI1 * ( t * C - _hypot(x-a1,y-a2) / LA ) + PI_1 )
    

    4.05

    schlecht. zu viele doubles involviert, was?

    mal LA bis LF als kehrwerte speichern, denn multiplizieren ist feiner als dividieren.

    double LA=1.0/10;
    double LB=1.0/10;
    ...
    				sin ( PI1 * ( t * C - sqrt2 ( (x - a1) * (x - a1) + (y - a2) * (y - a2) ) * LA ) + PI_1 ) 
    				+    sin ( PI2 * ( t * C - sqrt2 ( (x - b1) * (x - b1) + (y - b2) * (y - b2) ) * LB ) + PI_2 ) 
    ...
    

    3.11

    für die übersicht ein eigenes hypot machen.

    inline double hypot(int x,int y){
    	return sqrt2(x*x+y*y);
    }
    ...
    void test(int s,double t){
    	for(int y=0;y<200;y++)        // 200 Zeilen 
    	{ 
    		for(int x=0;x<200;x++)    // 200 Spalten 
    		{ 
    			z0[s]=floor(((    
    					sin ( PI1 * ( t * C - hypot(x-a1,y-a2) * LA ) + PI_1 ) 
    				+    sin ( PI2 * ( t * C - hypot(x-b1,y-b2) * LB ) + PI_2 ) 
    				+    sin ( PI2 * ( t * C - hypot(x-c1,y-c2) * LB ) + PI_2 ) 
    				+    sin ( PI2 * ( t * C - hypot(x-d1,y-d2) * LB ) + PI_2 ) 
    				+    sin ( PI2 * ( t * C - hypot(x-e1,y-e2) * LB ) + PI_2 ) 
    				+    sin ( PI2 * ( t * C - hypot(x-f1,y-f2) * LB ) + PI_2 ) 
    					) * 255/Anzahl) + 0.5 );  
    		}
    	}
    }
    

    3.10

    code kleinmachen tut doch gar nicht so weh.

    t*C hat der compiler sicherlich schon aus den schleifen herausgezogen. zur sicherheit mache ich es per hand.

    double tC=t*C;
    	for(int y=0;y<200;y++)        // 200 Zeilen 
    	{ 
    		for(int x=0;x<200;x++)    // 200 Spalten 
    		{ 
    			z0[s]=floor(((    
    					sin ( PI1 * ( tC - hypot(x-a1,y-a2) * LA ) + PI_1 )
    

    3.05

    PI1 könnte vorab mit tC und LA plutimiziert werden.

    3.03

    double LA=1.0/10*PI1;
    double LB=1.0/10*PI2;
    double LC=1.0/10*PI3;
    double LD=1.0/10*PI4;
    double LE=1.0/10*PI5;
    double LF=1.0/10*PI6;
    ...
    void test(int s,double t){
    	double tC1=t*C*PI1;
    	double tC2=t*C*PI2;
    	double tC3=t*C*PI3;
    	double tC4=t*C*PI4;
    	double tC5=t*C*PI5;
    	double tC6=t*C*PI6;
    	for(int y=0;y<200;y++)        // 200 Zeilen 
    	{ 
    		for(int x=0;x<200;x++)    // 200 Spalten 
    		{ 
    			z0[s]=floor(((    
    					sin ( ( tC1 - hypot(x-a1,y-a2) * LA ) + PI_1 ) 
    				+    sin ( ( tC2 - hypot(x-b1,y-b2) * LB ) + PI_2 ) 
    				+    sin ( ( tC3 - hypot(x-c1,y-c2) * LC ) + PI_2 ) 
    				+    sin ( ( tC4 - hypot(x-d1,y-d2) * LD ) + PI_2 ) 
    				+    sin ( ( tC5 - hypot(x-e1,y-e2) * LE ) + PI_2 ) 
    				+    sin ( ( tC6 - hypot(x-f1,y-f2) * LF ) + PI_2 ) 
    					) * 255/Anzahl) + 0.5 );  
    		}
    	}
    }
    

    jetzt noch schauen, om man ne schnellere sqrt() im netz findet. eine, die für integers ist?

    oder auch nachgucken!

    andere idee:

    schau mal die ausgabe von

    double x=1,y=0;
    double d=0.001;
    for(;;){
       cout<<x<<' '<<y<<endl;
       x+=y*d;
       y-=x*d;
    }
    

    an.

    sinuskurve ganz ohne sin. an den startwerten und d kann man frequenz, phase und amplitude verstellen. vielleicht sollte man 200*200*6 soche dinger benutzen?

    oder gibt es ein feines dynamisches system, das eine welle in 2d erzeugt? und davon 6 parallel machen und dann erst addieren?



  • ( (x - a1) * (x - a1) + (y - a2) * (y - a2) )
    

    das kannst du doch inkrementell ausrechnen,d.h.

    (y - a2) * (y - a2) brauchst du nur zu berechnen wenn sich y ändert.

    Außerdem:

    Es sei
    wert_alt = (x - a1) * (x - a1)
    dann macht "for" irgendwann
    x++
    somit ist dann
    wert_new = (x+1-a1) * (x+1-a1)
    = 2(x-a)+1 + wert_alt
    = ((x-a)<<1)+ 1 + wert_alt
    

    [EDIT]
    Ach ja, nach Möglichkeit casten zwischen int und float vermeiden
    eventuell sogar den Schleifenzähler als float.

    Viele Grüße
    Fischi



  • speichere zu jedem punkt die abstände zu den zentren in einem double abstand[200][256][8].



  • MgcReal MgcFastFunction::Sin0 (MgcReal fT)
    {
        // assert:  0 <= fT <= PI/2
        // maximum absolute error = 1.6415e-04
        // speedup = 1.91
    
        MgcReal fTSqr = fT*fT;
        MgcReal fResult = 7.61e-03;
        fResult *= fTSqr;
        fResult -= 1.6605e-01;
        fResult *= fTSqr;
        fResult += 1.0;
        fResult *= fT;
        return fResult;
    }
    //----------------------------------------------------------------------------
    MgcReal MgcFastFunction::Sin1 (MgcReal fT)
    {
        // assert:  0 <= fT <= PI/2
        // maximum absolute error = 2.3279e-09
        // speedup = 1.40
    
        MgcReal fTSqr = fT*fT;
        MgcReal fResult = -2.39e-08;
        fResult *= fTSqr;
        fResult += 2.7526e-06;
        fResult *= fTSqr;
        fResult -= 1.98409e-04;
        fResult *= fTSqr;
        fResult += 8.3333315e-03;
        fResult *= fTSqr;
        fResult -= 1.666666664e-01;
        fResult *= fTSqr;
        fResult += 1.0;
        fResult *= fT;
        return fResult;
    }
    

    aus der magic engine
    mfg



  • also erstmal danke dass ich euch so viele gedanken zu meiner frage macht, ich werd mal schaun wsa ich davon verwenden kann, vieles versteh ich etz auf anhieb gar net (nach 1mal durchlesen), bin aber auch grad erst aufgestanden...

    zu volkard wollt ich noch sagen, dass ich die werte von z0 auf jeden fall in einem integer feld speichere um platz zu sparen, vlt lässt sich ja auch schneller mit rechnen...
    und zwar gleich aufm heap:

    int *z0 = new int[40000];
    

    ansonsten werd ich mir mal ein paar gedanken über eure vorschläge machen und ein bisschen rumtesten... was ich jetzt schon mal gut find, is die wellenlängen gleich als kehrwert zu speichern, das scheint ja bei dir sehr viel gebracht zu haben und des hört sich auch sehr logisch an... ansonstens dauerts einfach noch ein bisschen bis ich hinter des ganze komm... ich sag bescheid inwiefern sich mein code ändert.

    achja und der compiler ist visual c++ 6

    Sirstefen



  • warum nicht ne LUT auf den Sinus vom Quadrat des Abstandes zum Wellenzentrum?
    das wird ein integer lookup weil ja dx und dy auch nur int sind oder erlaubst du float für Wellenzentrum?
    die grösse würd nur dx_max2+dy_max2 sein.
    Mann kann so jede Welle vorberechnen in Radiusrichtung und dan direckt aus den deltas den zugehörigen höhenwert ablesen



  • die wellenzentren sind nur integer werte, halt von 0 bis 200 möglich

    aber ich versteh net ganz was du meinst....



  • Sirstefen schrieb:

    die wellenzentren sind nur integer werte, halt von 0 bis 200 möglich

    aber ich versteh net ganz was du meinst....

    Die Intensität der WElle ist nur abhängig vom Abstand zum zentrum (bei konzentrischen Wellen) da kann man sich pro Welle einfach nen Vektor mit allen Intensitäten ausrechnen, in welche man für jeden Abstand die dazugehörige intensität speichert. nun muss man beim lookup nicht unbedingt den Abstand selber nehmen sonder kann gleich das Quadrat benutzen, man spart die Wurzel.



  • Edit: gegenstandslos
    

    hier mal eine arbeitsversion, etwa 20mal schneller, als das original (auf p4 optimiert aber ohne sse). sqrt wird durch eine LUT ersetzt, ausserdem befinden sich in der innersten schleife nur noch simple operationen bis auf sin. floor ist auch eine sehr teuere operation, mir ist keine round funktion bekannt, die analog zu floor arbeitet, aber eben normal rundet. wegen der beschränkungen des inline assmblers von vc kann man diese auch nicht so ohne weiteres performant hinzufügen (obwohl sie sehr simpel ist), mit gcc ist das kein problem. daher habe ich die gesamte for schleife mittels inline assembler formuliert. eine performante integer sinus version (mit festkommazahlen) wäre sehr wünschenswert. dann könnte man gänzlich auf gleitkommaberechnungen verzichten. dadurch, dass nun jede welle einzeln berechnent wird, ergeben sich weitere möglichkeiten zum optimieren. (zumal es jetzt auch schneller geht, wenn eine quelle nicht existiert). man könnte hier von den symmetrieeigenschaften der welle gebrauch machen. im idealfall (zentrum in der mitte), reduziert sich dann der aufwand auf 1/8.

    double brauchen wir nicht, da zum schluss sowieso gerundet wird, genügt float völlig. dann reicht es aber, die fpu gleich anzuweisen, nur mit veringerter genauigkeit zu arbeiten, daher das controlfp.



  • Edit: gegenstandslos
    

    hier das ganze mal mit festkommaarithmetik. ist noch etwas schneller als die vorhergehende version. ich hab aber keine ahnung ob sie korrekt funktioniert 🙂
    sie könnte wesentlich schneller sein, wenn die grösse der LUTs reduziert wird. idealerweise sollten beide wenigstens in den L2 cache passen, was im moment nicht der fall ist ( 1MB+310KB ). Ausserdem kann man den einzelnen Wellen noch relative wichte, also eine intensität beigeben.



  • also erstmal vielen dank für deine mühe!!

    hast du beim ersten beispiel nur die zeit getestet oder auch ob des ganze wirklich auch korrekt berechnet ist? (also natürlich net 100% korrekt, is ja in der physik eh nie, aber es sollte halt sagen wir mal net mehr als 30% oder so vom ergebnis abweichen)...
    wenn des so wär dann würd ich mir nämlich mal die mühe machen mich damit auseinander zu setzen und des ganze in mein programm zu implementieren (was halt sehr lange dauern würde, weil ich des ganze mit der winapi (also halt für windows) mache und da ewig viele buttons und variablen hab um des zu starten, außerdem müsste ich erstmal alles verstehen was du so bei deiner rechnung machst... und müsste mich da mal einarbeiten)
    und ich hatte jetzt schon des problem dass bei volkards voschlag aus dem dividieren ein multiplizieren zu machen das ganze nimmer funktioniert hat und nur noch ein schwarzer bildschirm, statt eben dem wellenmuster ausgegeben wurde. obwohl ich von der formel her alles richtig eingegeben habe und ich auch sonst nix geändert hab... hab keine ahnung was da des problem is...

    aber wenns von den werten her eben bei dir funktionieren würde und des ganze wirklich 20 mal (oder vlt auch nur 5-10 mal) schneller wär, dann würd ich des echt mal testen...

    Sirstefen



  • Sirstefen schrieb:

    und ich hatte jetzt schon des problem dass bei volkards voschlag aus dem dividieren ein multiplizieren zu machen das ganze nimmer funktioniert hat und nur noch ein schwarzer bildschirm

    du hast nicht double LA=1.0/10; geschrieben, sondern double LA=1/10; ,oder?



  • hier eine fehlerbereinigte version, sie braucht (auf P4) etwa 25 takte je punkt und welle. damit ist sie etwa 12mal schneller als das original. allerdings vermute ich ohnehin, das die langsame darstellung nicht durch die berechnung (die war im original schon schnell genug für 10 bilder/s bei 800MHz) sondern das zeichnen verursacht wird.

    #include <cmath>
    #include <iostream>
    #include <cassert>
    
    typedef unsigned short u16;
    typedef unsigned int u32;
    typedef unsigned long long u64;
    
    const int num_sources = 6;
    const int num_rows = 200;
    const int num_columns = 200;
    
    const int num_values = num_rows * num_columns;
    const int max_sqrt_arg = ( num_rows - 1 ) * ( num_rows - 1 ) + ( num_columns - 1 ) * ( num_columns - 1 );
    
    int z[ num_values ];
    
    double const PI = 3.1415926535897932384626433832795;
    
    double C = 1;
    
    int WZ[ num_sources ][ 2 ] = { { 10, 10 }, { 40, 10 }, { 70, 10 }, { 100, 10 }, { 130, 10 }, { 160, 10 } };
    bool active[ num_sources ] = { true, true, true, true, true, true };
    double phase[ num_sources ] = { PI, PI, PI, PI, PI, 0 };
    double intensity[ num_sources ] = { 2, 1, 1, 1, 1, 1 };
    
    double L[ 6 ] = { 50, 50, 50, 50, 50, 50 };
    
    u16 sqrtTab[ max_sqrt_arg + 1 ];
    short sinTab[ 65536 ];
    
    u16* initSqrtTab()
    {
        for( int i = 0; i <= max_sqrt_arg; ++i ) sqrtTab[ i ] = static_cast< u16 >( 128 * sqrt( double( i ) ) );
        return sqrtTab;
    }
    
    short* initSinTab()	// jede Periode wird in 65536 Teile geteilt, damit kann das notwendige %
    {					// beim abruf durch schnelleres & ersetzt werden
    	for( int i = 0; i < 65536; ++i )
    	{				// gespeichert wird der das 32767 fache des sinus
    		sinTab[ i ] = static_cast< short >( floor( 32767.4999 * sin( 2 * PI * double( i ) / 65536.0 ) + 0.5 ) );
    	}
    	return sinTab;
    }
    
    u16* sqrtTabp = initSqrtTab();
    short* sinTabp = initSinTab();
    
    u32 sqr2(const int x)
    {
    	return u32( x * x );
    }
    
    u32 sqrt2(const u32 x)
    {
    	assert( x <= max_sqrt_arg );
    	return sqrtTab[ x ];
    }
    
    int sin2(const u32 x)
    {
    	return sinTab[ x % 65536 ];
    }
    
    void test(double t)
    {
    	memset( z, 0, sizeof( z ) );
    	double total_intensity = 0;
    	for( int k = 0; k < num_sources; ++k ) if( active[ k ] ) total_intensity += intensity[ k ];
    	for( int k = 0; k < 6; ++k )
    	{
    		if( active[ k ] )
    		{
    			int rel_weight = int( 255.9999 * intensity[ k ] / total_intensity );
    			u32 current_phase;
    			{
    				double cp = fmod( phase[ k ] / ( 2 * PI ) + C * t / L[ k ], 1 );
    				if( cp < 0 ) cp += 1;
    				current_phase = static_cast< u32 >( 65536.0 * cp );
    			}
    			u32 Lk;
    			{
    				double lk = - 65536.0 / L[ k ];
    				if( lk < 0 )
    				{
    					lk = - lk;
    					current_phase = 65536u - current_phase;
    				}
    				Lk = static_cast< u32 >( lk );
    			}
    			u32 s = 0;
    			for( int y = 0; y < num_rows; ++y )
    			{
    				u32 xySqr = sqr2( WZ[ k ][ 0 ] ) + sqr2( y - WZ[ k ][ 1 ] );
    				int xWZ = - WZ[ k ][ 0 ];
    				for( int x = 0; x < num_rows; ++ x )
    				{
    					z[ s++ ] += rel_weight * sin2( current_phase + Lk * sqrt2( xySqr ) / 128 );
    					xySqr += xWZ + xWZ + 1;
    					++xWZ;
    				}
    			}
    		}
    	}
    	for( u32 s = 0; s < num_values; ++s ) z[ s ] /= 32768;
    }
    
    u64 gtime=0;
    
    inline void starttimer()
    {
    	__asm
    	{
    		rdtsc
    		mov		dword ptr gtime, eax
    		mov		dword ptr gtime + 4, edx
    	}
    }
    
    inline void endtimer()
    {
    	__asm
    	{
    		rdtsc
    		sub		eax, dword ptr gtime
    		sbb		edx, dword ptr gtime + 4
    		mov		dword ptr gtime, eax
    		mov		dword ptr gtime + 4, edx
    	}
    }
    
    int main()
    {
    	for( int s = 0; s < 100; ++s ) test( 1.0f * s / 100 );
    	int k = 10000;
    	u64 besttime = ~0ULL;
    	do
    	{
    		starttimer();
    	    test( 0.01f * k );
    		endtimer();
    		if( gtime < besttime )
    		{
    			besttime = gtime;
    			k = 10000;
    		}
    	}
    	while( --k );
    	std::cout << double( besttime ) / double( num_values ) / double( num_sources ) << std::endl;
        return EXIT_SUCCESS;
    }
    


  • volkard schrieb:

    Sirstefen schrieb:

    und ich hatte jetzt schon des problem dass bei volkards voschlag aus dem dividieren ein multiplizieren zu machen das ganze nimmer funktioniert hat und nur noch ein schwarzer bildschirm

    du hast nicht double LA=1.0/10; geschrieben, sondern double LA=1/10; ,oder?

    ja genau... danke


Anmelden zum Antworten