slawek <s...@h...pl> napisał(a):

> Jak chcesz się bawić - to ok.

Pobawiłem się na czterech funkcjach. Aproksymuję wielomianem 2-go stopnia i
trapezami. Nie mam pewności czy w aproksymacji wielomianem nie mam błędu.
Na pierwszy rzut oka stosunek dokładności do czasu wygląda niemal tak
samo dla wielomianu i trapezów. Ale zapewne aproksymację wielomianem
da się wydajniej zakodować. Kod (pod winde) umieszczam poniżej. Da się
go łatwo rozbudować na inne funkcje i inne sposoby całkowania.

> Na przykład za pomocą FFT? To chyba samo w sobie ciekawe pytanie.
Czemu nie?

Pozdrawiam

typedef double tf; // type float

#define d_prec_min (1<<2)
#define d_prec_max (1<<19)
#define d_cnt_test (500)

// Pomiar czasu
class Times
{
private:
__int64 v_start;
public:
void start(void)
{
FILETIME tmp1,tmp2,tmp3,tmp4;
GetThreadTimes( GetCurrentThread() , &tmp1 , &tmp2 , &tmp3 , &tmp4 );
v_start = *((__int64*)&tmp4);
}
int end(int cnt_test)
{
FILETIME tmp1,tmp2,tmp3,tmp4;
__int64 v_end;
GetThreadTimes( GetCurrentThread() , &tmp1 , &tmp2 , &tmp3 , &tmp4 );
v_end = *((__int64*)&tmp4);
return (int)((v_end-v_start)/10/cnt_test);
}
};

// Baza dla testowanych funkcji
class BaseFunc
{
friend class BaseTest;
virtual tf exactValue(void) const = 0; // dokładna wartość
virtual tf start(void) const = 0; // początek zakresu
virtual tf end(void) const = 0; // koniec zakresu
virtual tf func(const tf x) const = 0; // funkcja
};

// sin(x)/x
class Si : public BaseFunc
{
public:
virtual tf exactValue(void) const { return
((tf)1.608194499263757628296323926138231688753059204
3316413876215625399831498); }
virtual tf start(void) const { return ((tf)0.01); }
virtual tf end(void) const { return ((tf)15); }
virtual tf func(const tf x) const { return ((tf)(sin(x)/x)); }
};

// cos(x)/x
class Ci : public BaseFunc
{
public:
virtual tf exactValue(void) const { return
((tf)4.074258198656752511848286445453690152144932604
98723926233961338683); }
virtual tf start(void) const { return ((tf)0.01); }
virtual tf end(void) const { return ((tf)15); }
virtual tf func(const tf x) const { return ((tf)(cos(x)/x)); }
};


// sin(x)/x^2
class Sii : public BaseFunc
{
public:
virtual tf exactValue(void) const { return
((tf)5.026408327223266870735360867169039335334048495
8641615625132381406263886425955984186496153422221141
36);
}
virtual tf start(void) const { return ((tf)0.01); }
virtual tf end(void) const { return ((tf)25); }
virtual tf func(const tf x) const { return ((tf)(sin(x)/x/x)); }
};

// 3xxx-2xx+5x+4
class Xxx : public BaseFunc
{
public:
virtual tf exactValue(void) const { return
((tf)7123.333333333333333333333333333333333333333333
3333333333333); }
virtual tf start(void) const { return ((tf)0.0); }
virtual tf end(void) const { return ((tf)10); }
virtual tf func(const tf x) const { return ((tf)(3*x*x*x-2*x*x+5*x+4)); }
};

// Baza dla testów
class BaseTest
{
private:
BaseFunc *bf; // testowana funkcja
tf *vals; // wartości
int cnt_test; // ilość testów (pomiarów czasu)
int prec_min; // minimalna dokładność
int prec_max; // maksymalna dokładność
tf step; // stały krok pomiaru
private:
void makeValue(const int prec)
{
vals = (tf*)realloc(vals,sizeof(tf)*prec);
step = (bf->end() - bf->start()) / (prec-1);
tf x = bf->start() - step;
for(int i=0 ; i<prec ; i++ )
vals[i] = bf->func( x += step );
}
protected:
virtual tf integral( const tf *vals , const tf step , int cnt ) = 0; // całka
public:
BaseTest(BaseFunc *bf,int prec_min=d_prec_min,int prec_max=d_prec_max,int
cnt_test=d_cnt_test):bf(bf),prec_min(max(prec_min,2)
),prec_max(prec_max),cnt_test(cnt_test),vals(NULL){}
virtual ~BaseTest() { free(vals); }
void go(void)
{
int prec = prec_min/2;
do
{
Times t;
tf comps;
prec *= 2;
makeValue( prec );
t.start();
for(int i=0 ; i<cnt_test ; i++ )
comps = integral( vals , step , prec );
printf("%10d %6d
%13.10lf\n",(int)prec,(int)t.end(cnt_test),(double)f
abs(comps-bf->exactValue()));
}
while( prec < prec_max );
}
};

class Trapezy : public BaseTest
{
protected:
virtual tf integral(const tf *vals , const tf step , int cnt)
{
tf sum = 0;
cnt--;
for(int i=1;i<cnt;i++)
sum += vals[i];
sum += (vals[0] + vals[cnt]) / 2;
return sum * step;
}
public:
Trapezy(BaseFunc *bf,int prec_min=d_prec_min,int prec_max=d_prec_max,int
cnt_test=d_cnt_test):BaseTest(bf,prec_min,prec_max,c
nt_test){}

};

class Polyn_2 : public BaseTest
{
private:
void abc( const tf y1, const tf y2, const tf y3, tf *a, tf *b, tf *c )
{
*c = y1;
*b = (-4*y2+4*y1+y3-y1) / (-2);
*a = (y2 - y1 - *b ) ;
}
protected:
virtual tf integral(const tf *vals , const tf step , int cnt)
{
tf a,b,c;
tf sum = 0;
for(int i=0;i<cnt-2;i++)
{
abc(vals[i],vals[i+1],vals[i+2],&a,&b,&c);
sum += a/3 + b/2 + c;
/*
if( fabs(a*0*0 + b*0 + c - vals[i] ) > 0.000000001 )
abort();
if( fabs(a + b + c - vals[i+1]) > 0.000000001 )
abort();
if( fabs( a*2*2 + b*2 + c - vals[i+2]) > 0.000000001 )
abort();
*/
}
abc(vals[cnt-3],vals[cnt-2],vals[cnt-1],&a,&b,&c);
sum += a*2*2*2 / 3 + b*2*2 / 2 + c*2 - a/3 - b/2 - c ;


return sum * step;
}
public:
Polyn_2(BaseFunc *bf,int prec_min=d_prec_min,int prec_max=d_prec_max,int
cnt_test=d_cnt_test):BaseTest(bf,prec_min,prec_max,c
nt_test){}

};

int main(int argc, char* argv[])
{
Si si;
Trapezy t1(&si);
t1.go();
Polyn_2 p1(&si);
p1.go();

Ci ci;
Trapezy t2(&ci);
t2.go();
Polyn_2 p2(&ci);
p2.go();

Sii sii;
Trapezy t3(&sii);
t3.go();
Polyn_2 p3(&sii);
p3.go();

Xxx xxx;
Trapezy t4(&xxx);
t4.go();
Polyn_2 p4(&xxx);
p4.go();

return 0;
}





--
Wysłano z serwisu Usenet w portalu Gazeta.pl -> http://www.gazeta.pl/usenet/