33,008
社区成员
发帖
与我相关
我的任务
分享
/*
* cosine_sine.cpp
*
* Created on: 2009-8-28
* Author: Administrator
*/
//1.27323954 = 4/pi
//0.405284735 =-4/(pi^2)
/*********************************************************
* high precision sine/cosine
*********************************************************/
double shift(const double x) {
//always wrap input angle to -PI..PI
while (x < -3.14159265)
x += 6.28318531;
while (x > 3.14159265)
x -= 6.28318531;
return x;
}
double sine(const double x) {
double sin = 0;
x = shift( x );
//compute sine
if (x < 0) {
sin = 1.27323954 * x + .405284735 * x * x;
if (sin < 0)
sin = .225 * (sin * -sin - sin) + sin;
else
sin = .225 * (sin * sin - sin) + sin;
} else {
sin = 1.27323954 * x - 0.405284735 * x * x;
if (sin < 0)
sin = .225 * (sin * -sin - sin) + sin;
else
sin = .225 * (sin * sin - sin) + sin;
}
return sin;
}
double cosine(const double x) {
double cos = 0;
//compute cosine: sin(x + PI/2) = cos(x)
x += 1.57079632;
if (x > 3.14159265)
x -= 6.28318531;
if (x < 0) {
cos = 1.27323954 * x + 0.405284735 * x * x;
if (cos < 0)
cos = .225 * (cos *-cos - cos) + cos;
else
cos = .225 * (cos * cos - cos) + cos;
}
else
{
cos = 1.27323954 * x - 0.405284735 * x * x;
if (cos < 0)
cos = .225 * (cos *-cos - cos) + cos;
else
cos = .225 * (cos * cos - cos) + cos;
}
return cos;
}
#include <stdio.h>
#include <math.h>
#include <windows.h>
#define PI 3.1415926535
inline double sin_tay(double x)
{
double x2 = x*x;
const double a1 = 0.99999748719198;
const double a3 = -0.166651680984443;
const double a4 = 0.00830951691334793;
const double a5 = -0.000184472198573026;
return x*( a1+ ( a3 + ( a4 + a5*x2 ) * x2 ) * x2);
}
class CTimer
{
public:
__forceinline CTimer( void ) {
QueryPerformanceFrequency( (PLARGE_INTEGER)&m_nFreq );
QueryPerformanceCounter( (PLARGE_INTEGER)&m_nStart ); }
__forceinline void Reset( void ) {
QueryPerformanceCounter( (PLARGE_INTEGER)&m_nStart ); }
__forceinline double Stop( void ) {
QueryPerformanceCounter( (PLARGE_INTEGER)&m_nCur );
return double( m_nCur - m_nStart ) / double( m_nFreq ); }
__forceinline bool SelfTest( void ){
return ( 0 != QueryPerformanceFrequency( (PLARGE_INTEGER)&m_nFreq ) ); }
private:
__int64 m_nCur;
__int64 m_nFreq;
__int64 m_nStart;
};
int main(){
double a=3*PI/8,b;
int i;
CTimer tmr;
tmr.Reset();
for(i=0;i<1000000000;i++){
GetTickCount();
}
printf("use time: %f s\n\n",tmr.Stop());
tmr.Reset();
for(i=0;i<1000000000;i++){
b=sin_tay(a);
}
printf("use time: %f s\n",tmr.Stop());
printf("sin_tay(%f) = %f\n\n",a,b);
tmr.Reset();
for(i=0;i<1000000000;i++){
b=sin_tay(a);
GetTickCount();
}
printf("use time: %f s\n",tmr.Stop());
printf("sin_tay(%f) = %f\n\n",a,b);
return 0;
}
use time: 6.378010 s
use time: 50.226837 s
sin_tay(1.178097) = 0.923879
use time: 59.875487 s
sin_tay(1.178097) = 0.923879
Press any key to continue
sum += X[j](x[i]) * X[k](x[i]) / ( sigma[i]*sigma[i] );
sum += y[i]*X[k](x[i]) / ( sigma[i]*sigma[i] );
#ifndef _LINEAR_LEAST_SQUARE_FIT_HPP_INCLUDED
#define _LINEAR_LEAST_SQUARE_FIT_HPP_INCLUDED
//
//to link:
// g++ -o xx xx.o -L/home/feng/lib -lgsl -lgslcblas -llinear_least_square_fit -llu_solve
//
//General Linear Least Square Fit Algorithm
//
//Tips:
// to fit a set of data points(x_{i}, y_{i}) to a model that
// is just a linear combination of any M specified
// functions of x -- X(x), the model is
// y(x) = \sum_{k=0}^{M} a_{k} X_{k}(x)
//
//Inputs:
// x -- x_{i}
// y -- y_{i}
// X -- (X{k})(x)
//
//Output: vector<double> linear_least_square_fit()const
// the fitted parameters
// a -- a_{k}
//
//Last Edit:
// Wed Jan 28 11:29:27 CST 2009
//
#include <cassert>
#include <vector>
#include <algorithm>
using std::vector;
class Linear_Least_Square_Fit
{
public:
Linear_Least_Square_Fit(
const vector<double>& _x,
const vector<double>& _y,
const vector<double(*)(double)>& _X
)
: x(_x), y(_y), X(_X)
{
assert( x.size() == y.size() );
assert( x.size() >= X.size() );
sigma.resize( x.size() );
fill_n( sigma.begin(), x.size(), 1.0 );
}
void set_sigma( const vector<double>& _sigma )
{
assert( sigma.size() == _sigma.size() );
sigma = _sigma;
}
vector<double> linear_least_square_fit() const;
private:
vector<double> x;
vector<double> y;
vector<double(*)(double)> X;
vector<double> sigma; //optional
};
#endif //_LINEAR_LEAST_SQUARE_FIT_HPP_INCLUDED
//Date: Tue Jun 23 13:53:43 CST 2009
//
//
#include <linear_least_square_fit.hpp>
#include <lu_solve.hpp>
#include <vector>
using namespace std;
vector<double> Linear_Least_Square_Fit::linear_least_square_fit()const
{
const unsigned int M = X.size();
const unsigned int N = x.size();
double alfa[M][M];
double beta[M];
for ( unsigned int k = 0; k < M; ++k )
for ( unsigned int j = 0; j < M; ++j )
{
double sum = 0.0;
for ( unsigned int i = 0; i < N; ++i )
{
sum += X[j](x[i]) * X[k](x[i]) / ( sigma[i]*sigma[i] );
}
alfa[k][j] = sum;
}
for ( unsigned int k = 0; k < M; ++k )
{
double sum = 0;
for ( unsigned int i = 0; i < N; ++i )
{
sum += y[i]*X[k](x[i]) / ( sigma[i]*sigma[i] );
}
beta[k] = sum;
}
double* fake_alfa = (double*) alfa;
const vector<double> a( fake_alfa, fake_alfa+M*M );
const vector<double> b( beta, beta+M );
const vector<double> ans = Lu_Solve( a, b ).lu_solve();
return ans;
}
#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#define PI 3.1415926535
inline double sin_tay(double x)
{
double x2 = x*x;
const double a1 = 0.99999748719198;
const double a3 = -0.166651680984443;
const double a4 = 0.00830951691334793;
const double a5 = -0.000184472198573026;
return x*( a1+ ( a3 + ( a4 + a5*x2 ) * x2 ) * x2);
}
int main(){
double a=3*PI/8,b;
int i;
struct timeval t1,t2;
float f;
gettimeofday(&t1,NULL);
for(i=0;i<10000000;i++){
b=sin(a);
gettimeofday(NULL,NULL);
}
gettimeofday(&t2,NULL);
f=(t2.tv_sec-t1.tv_sec)*1000000+t2.tv_usec-t1.tv_usec;
f/=1000000;
printf("use time: %f s\n",f);
printf("sin(%f) = %f\n\n",a,b);
gettimeofday(&t1,NULL);
for(i=0;i<10000000;i++){
b=sin_tay(a);
gettimeofday(NULL,NULL);
}
gettimeofday(&t2,NULL);
f=(t2.tv_sec-t1.tv_sec)*1000000+t2.tv_usec-t1.tv_usec;
f/=1000000;
printf("use time: %f s\n",f);
printf("sin_tay(%f) = %f\n\n",a,b);
return 0;
}
linren@linren:~/work/select$ cc 5.c -lm -o 5
linren@linren:~/work/select$ ./5
use time: 4.766741 s
sin(1.178097) = 0.923880
use time: 4.476374 s
sin_tay(1.178097) = 0.923879
linren@linren:~/work/select$