哪位给个sin/cos的优化版?(C/C++)

fallening 2009-08-28 03:03:53
程序中调用了20亿次的样子,要将近80多秒(1.3G cpu)的运行时间
哪位收藏有更好的?
偶运行的程序大致如下(实际使用的修修补补过,不过并没有很大的变化,也就是少做一两次乘法的样子)
如果哪位有基于查表法之类的,请不要吝啬



/*
* 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;
}
...全文
611 点赞 收藏 76
写回复
76 条回复
切换为时间正序
当前发帖距今超过3年,不再开放新的回复
发表回复
中華大帝 2011-05-11
sin_tay这个函数,在360度绘制圆的时候,超过了180度,就开始偏离轨道了。最终绘制的图像不是个圆,而是一个类似问号(?)的图形,这个怎么解决?

难道还要加上
if(角度<180)
else
...
这样么?这样的话,效率也不会太高了啊。
回复
fallening 2009-08-30
呵呵,多谢

我引入了巨量的全局变量,虽然把代码搞得非常丑陋,但是却成功地将运行时间降到了10s以下

结了:)

回复
linren 2009-08-30
上面发错链接了……
https://forum.csdn.net/BList/CUDA
回复
linren 2009-08-30
[Quote=引用 71 楼 fallening 的回复:]
cuda???
[/Quote]如果是CPU的话
在linux/unix下可以用pthread来并行……
不过或许还是GPU的并行好一些?

CUDA就是一种利用GPU来做一些CPU做的事情……
并且还能做的更好的这样一种概念……

CUDA有专版:
http://community.csdn.net/
不过没有怎么接触过……
回复
fallening 2009-08-30
cuda???
回复
churenxh 2009-08-30
学习了
回复
linren 2009-08-29
我的电脑的GPU是ATI的……
CUDA……
这个真没有(泪)……
回复
linren 2009-08-29
[Quote=引用 63 楼 fallening 的回复:]
cpu缓存怎么清空?
我把一部分调用sin/cos函数的代码贴出来你帮我看看

C/C++ code#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

C/C++ code//Date: Tue Jun 23 13:53:43 CST 2009////
#include<linear_least_square_fit.hpp>
#include<lu_solve.hpp>

#include<vector>usingnamespace std;


vector<double> Linear_Least_Square_Fit::linear_least_square_fit()const
{const unsignedint M= X.size();const unsignedint N= x.size();double alfa[M][M];double beta[M];for ( unsignedint k=0; k< M;++k )for ( unsignedint j=0; j< M;++j )
{double sum=0.0;for ( unsignedint i=0; i< N;++i )
{
sum+= X[j](x[i])* X[k](x[i])/ ( sigma[i]*sigma[i] );
}
alfa[k][j]= sum;
}for ( unsignedint k=0; k< M;++k )
{double sum=0;for ( unsignedint 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;
}
其中的那个Lu_Solve调用了blas解线性方程,已经无法再优化了
[/Quote]
没有看明白这段程序……

只是看懂了X[k](x)可能是:
sin(x), cos(x), sin(Ux), cos(Ux), sin(x)cos(Ux), cos(x)sin(Ux)
这几种类型的函数……

不过感觉上……
xingzhe2001说的“并行”可能就是楼主需要的答案了^_^
回复
xuboamei 2009-08-29
如果从提高程序的效率的角度,还是查表最快了,把sin cos算出来,下次要用直接取来
回复
linren 2009-08-29
[Quote=引用 61 楼 fallening 的回复:]
呃,但是在偶这边,自己的sine可比系统中的sin快了不止一点半点啊
[/Quote]
性能优化可能只能具体情况具体分析……
C库的差异
编译器的差异
操作系统的差异
进程的运行环境
……
以上都会对结果产生很大的影响……

[Quote=引用 62 楼 xingzhe2001 的回复:]
可是用了GetTickCount以后占用时间的大头就变成了GetTickCount, sin占用的时间只是个小头,所以两个时间看其来差别很小。这也是我们平常优化的时候为什么要优化占用时间最长的函数。

用每次循环清空cpu指令和数据缓存的方法看看?
[/Quote]
GetTickCount和gettimeofday确实会占用一定的时间……

#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


不过可以通过上面的办法
用减法来算出一个大致的时间……
回复
fallening 2009-08-29
这边倒是有个六核的sun的机器,不过我的程序是为了发给别人用的,他们那边貌似机器都很老的样子
回复
xingzhe2001 2009-08-29
看起来适合并行阿,用多核多线程看看?
你们有钱的话用CUDA做并行也可以。
回复
xingzhe2001 2009-08-29
[Quote=引用 63 楼 fallening 的回复:]
cpu缓存怎么清空?
我把一部分调用sin/cos函数的代码贴出来你帮我看看
[/Quote]

清空缓存是为了更慢,代替GetTickCount消耗cpu。
回复
fallening 2009-08-29
上边调用的时候,N>>M,开销大头是
sum += X[j](x[i]) * X[k](x[i]) / ( sigma[i]*sigma[i] );


sum += y[i]*X[k](x[i]) / ( sigma[i]*sigma[i] );

回复
fallening 2009-08-29
cpu缓存怎么清空?
我把一部分调用sin/cos函数的代码贴出来你帮我看看

#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;
}


其中的那个Lu_Solve调用了blas解线性方程,已经无法再优化了
回复
xingzhe2001 2009-08-29
可是用了GetTickCount以后占用时间的大头就变成了GetTickCount, sin占用的时间只是个小头,所以两个时间看其来差别很小。这也是我们平常优化的时候为什么要优化占用时间最长的函数。

用每次循环清空cpu指令和数据缓存的方法看看?
回复
fallening 2009-08-29
呃,但是在偶这边,自己的sine可比系统中的sin快了不止一点半点啊
回复
linren 2009-08-29
(上接59楼)
【说明】
59楼中的是用虚拟机在Ubuntu下进行的测试……
因为是虚拟机的缘故
所以把1000000000次改成了10000000次……
回复
linren2 2009-08-29
#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$

回复
linren 2009-08-29
(上接57楼)

如果是在linux/unix下的话……
针对57楼中的小技巧
可以使用gettimeofday来替换GetTickCount……
回复
发动态
发帖子
数据结构与算法
创建于2007-08-27

3.2w+

社区成员

数据结构与算法相关内容讨论专区
申请成为版主
社区公告
暂无公告