sinus funktion vereinfachen bzw. schneller berechnen lassen



  • hi,
    ich hab folgende berechnung in meinem programm und die funktioniert auch bisher perfekt, allerdings muss ich sie 40000 durchführen weil ich ein fenster mit 200*200 Pixeln habe und ich die Berechnung für jeden einzelnen Punkt durchführen muss. Da ich aber die Simulation nicht nur als Standbild haben möchte, sondern sogar als eine art video, müsste das ganze auch noch 25 mal pro sekunde funktionieren und am besten auch auf einer 800 MHz CPU. Ich habs zwar noch nicht getestet, aber ich seh dass das Standbild auf meinem Athlon XP 2400+ (mit 2000MHz) schon ca. 1 Sekunde benötigt. Auf meinem Zweit PC (Duron mit 700 MHz) braucht sie etwa 2-3 Sekunden.

    Das ist die Formel:

    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 - sqrt ( (x - a1) * (x - a1) + (y - a2) * (y - a2) ) / LA ) + PI_1 ) 
    						+	sin ( PI2 * ( t * C - sqrt ( (x - b1) * (x - b1) + (y - b2) * (y - b2) ) / LB ) + PI_2 )
    						+	sin ( PI3 * ( t * C - sqrt ( (x - c1) * (x - c1) + (y - c2) * (y - c2) ) / LC ) + PI_3 )
    						+	sin ( PI4 * ( t * C - sqrt ( (x - d1) * (x - d1) + (y - d2) * (y - d2) ) / LD ) + PI_4 )
    						+	sin ( PI5 * ( t * C - sqrt ( (x - e1) * (x - e1) + (y - e2) * (y - e2) ) / LE ) + PI_5 )
    						+	sin ( PI6 * ( t * C - sqrt ( (x - f1) * (x - f1) + (y - f2) * (y - f2) ) / LF ) + PI_6 )    
    																												    ) * 255/Anzahl) + 0.5 );
    

    zur Erklärung:
    es handelt sich um die Berechnung der Intensität von Wellen für eine Interferenz Simulation mit 6 Wellenzentren, die beliebig aktivierbar/deaktivierbar (bis auf 0 Wellen, das ist verboten) sind... (das ganze ist die umgeformte Wellenfunktion die von der Zeit, der Wellenlänge und dem Abstand zum Wellenzentrum abhängig ist.)
    z0 ist ein int Feld mit eben 40000 Werten (was also max. ca. 80KB verbraucht, was also vom Speicherplatz her kein problem sein dürfte, vor allem, weil ich das Feld nach jedem Durchgang lösche)
    floor ist (soweit ich weiß) abrunden, deswegen am ende + 0,5
    PI1 bis PI6 ist jeweils einfach der doppelte PI Wert, der je nachdem ob die Welle aktiviert ist, eben 6,283... oder 0 ist
    t ist die Zeit, die bisher eben einfach nur 1 Sekunde ist
    C ist die Lichtgeschwindigkeit
    sqrt die Wurzel (für den Pythagoras zur Bestimmung des Abstandes des Pixels (x,y) vom Wellenzentrzm A(a1,a2) bis F(f1,f2))
    LA bis LF sind die Wellenlängen, die aufgrund des fensters von 200 pixeln zwischen 1 und 999 beliebig ganzzahlig vom benutzer wählbar ist.
    PI_1 bis PI_2 sind Änderungsfaktoren, die die Werte 0, 0.5*PI, PI oder 1,5*PI annehmen können, falls der Benutzer es wünscht, dass die Welle nicht mit der Intensität 0 sondern z.b. mit einem Maximum oder einem Minimum beginnt.
    " * 255 " deswegen, weil ich das ganze in farben darstelle, negative wert in rot (und eben umso heller, umso niedriger der Wert ist), positive blau (dasselbe wie bei rot)
    Anzahl ist ein Integer Wert, in dem einfach gespeichert ist, wie viele Wellen aktviert sind.
    Die einzelnen Sinus Werte werden addiert, damit die maximale Intensität der Anzahl der aktivierten Wellen entspricht (also max. 6) und damit max. 255 als farbwert gibt bzw. -6 und -255 für Rot (was dann durch ne if funktion nach der berechnung zugeordnet wird und wodurch die roten negativen werte auch wieder positiv werden...)

    So damit müsste eigentlich jedem so halbwegs klar sein, was das ganze bedeuten soll.

    Ich möchte nicht auf irgendwelche Features wie z.b. die Aktivierbarkeit/Deaktivierbarkeit verschiedener Wellen (PI1 bis PI6 und Anzahl), die Anfangsintensität (PI_1 bis PI_6)...

    Weiß jemand trotzdem wie ich die Berechnung vereinfachen kann oder so umstellen kann, damit sie schneller berechnet werden kann (evlt durch Algorithmen oder sonstige Sachen)?

    Wär dankbar für Hile Sirstefen



  • Leg dir einfach eine Tabelle für 1000 Werte zwischen 0 und 2Pi an. Du verzichtest dann auf die Berechnung und nimmst (nachdem du einen Wert auf dieses Intervall zurückgerechnet hast) einfach das Ergebnis des nächsten vorberechneten Wertes. Nachgucken geht nämlich schneller als rechnen. 🙂



  • 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


Log in to reply