正余弦的近似计算

G-Spider · 发表于 2011-2-22 20:41:50

10# G-Spider
通过上面的思路，有效精度14位，而以上弧度特性也正是fft中的sin,cos表。而算一个sin或cos值只需要两次浮点乘（加，移位不记），比较满意。

#include <stdio.h>
#include <math.h>
#define PI 3.1415926535897932
#define _2PI (2.0*PI)
double lpsin[256+1],lpcos[256+1];
int main()
{
double N1,theta;
unsigned int i,j,tmpi,N,N0;
j=10;
N=(unsigned int)1<<j;
N0=N>>2; //256
N1=1.0/(double)N;
theta=N1*_2PI;
//lpsin=(double *)malloc( (N0+1)*sizeof(double) );
//lpcos=(double *)malloc( (N0+1)*sizeof(double) );
//if(lpsin==NULL || lpcos==NULL)
//printf("Memory Allocated Error\n");
lpsin[0]=0.0;
lpcos[0]=1.0;
lpsin[1]=sin(theta);
lpcos[1]=cos(theta);
for(i=2;i<=N0;i+=2)
{//i=2,4,6,....
tmpi=i>>1;
lpsin[i]= 2.0* lpsin[tmpi] * lpcos[tmpi];
lpcos[i]=1.0- 2.0* lpsin[tmpi] * lpsin[tmpi];
//======================================================
lpsin[i+1]= lpcos[i] * lpsin[1] + lpsin[i] * lpcos[1];
lpcos[i+1]=-lpsin[i] * lpsin[1] + lpcos[i] * lpcos[1];
}
for(i=0;i<=N0;i++)
printf("%18.14lf %18.14lf\n",*(lpsin+i),*(lpcos+i));
system("PAUSE");
return 0;
}

复制代码

wayne · 发表于 2011-2-23 10:30:37

11# G-Spider
math.h 里面有常量 M_PI ，精度比你的还高呢

G-Spider · 发表于 2011-2-23 10:41:14

11# G-Spider
接上，对于上面的弧度严递增情形，也只需要算0~pi/4
//sin(pi/2 -x)=cos(x); sin(N0-i)=cos(i); N0对应pi/2时的下标，i 对应x的下标，x->0~pi/4
//cos(pi/2-x)=sin(x); cos(N0-i)=sin(i);
这样既不会过大的传递误差（事实上，通过用sin()从0~pi/2，与上面的得到了相同的16位精度），而对于偶数项存在访存相邻，不会有严重的Cache的抖动。
且相比上面的，减少了一半的浮点乘，对于x,2x,3x,....Nx 同时计算sin,cos值只需要2N次浮点乘即可。

#include <stdio.h>
#include <math.h>
#define PI 3.1415926535897932
#define _2PI (2.0*PI)
double lpsin[256+1],lpcos[256+1];
int main()
{
double N1,theta;
unsigned int i,j,tmpi,N,N0;
j=10;
N=(unsigned int)1<<j; //N=1024
N0=N>>3; //N0=1024/8=128
N1=1.0/(double)N;
theta=N1*_2PI;
//lpsin=(double *)malloc( (N0+1)*sizeof(double) );
//lpcos=(double *)malloc( (N0+1)*sizeof(double) );
//if(lpsin==NULL || lpcos==NULL)
//printf("Memory Allocated Error\n");
lpsin[0]=0.0;
lpcos[0]=1.0;
lpsin[1]=sin(theta);
lpcos[1]=cos(theta);
for(i=2;i<=N0;i+=2)
{//i=2,4,6,....
tmpi=i>>1;
lpsin[i]= 2.0* lpsin[tmpi] * lpcos[tmpi];
lpcos[i]=1.0- 2.0* lpsin[tmpi] * lpsin[tmpi];
//======================================================
lpsin[i+1]= lpcos[i] * lpsin[1] + lpsin[i] * lpcos[1];
lpcos[i+1]=-lpsin[i] * lpsin[1] + lpcos[i] * lpcos[1];
}
N0<<=1;
for(;i<=N0;i+=2)
{//i=N0'+2,N0'+4,N0'+6,....,N0 (N0'=N0/2)
//sin(pi/2 -x)=cos(x); sin(N0-i)=cos(i);
//cos(pi/2-x)=sin(x); cos(N0-i)=sin(i);
tmpi=N0-i;
lpsin[i]=lpcos[tmpi];
lpcos[i]=lpsin[tmpi];
//======================================================
tmpi--;
lpsin[i+1]=lpcos[tmpi];
lpcos[i+1]=lpsin[tmpi];
}
printf(" rad sin(rad) cos(rad)\n");
for(i=0;i<=N0;i++)
printf("%18.16lf %18.16lf %18.16lf\n",i*N1*_2PI ,*(lpsin+i),*(lpcos+i));
system("PAUSE");
return 0;
}

复制代码

部分结果：

rad sin(rad) cos(rad)
0.0000000000000000 0.0000000000000000 1.0000000000000000
0.0061359231515426 0.0061358846491545 0.9999811752826011
0.0122718463030851 0.0122715382857199 0.9999247018391445
0.0184077694546277 0.0184067299058048 0.9998305817958233
0.0245436926061703 0.0245412285229123 0.9996988186962043
0.0306796157577128 0.0306748031766366 0.9995294175010933
0.0368155389092554 0.0368072229413588 0.9993223845883495
0.0429514620607980 0.0429382569349408 0.9990777277526455
0.0490873852123405 0.0490676743274180 0.9987954562051724
0.0552233083638831 0.0551952443496899 0.9984755805732947
0.0613592315154256 0.0613207363022086 0.9981181129001492
0.0674951546669682 0.0674439195636641 0.9977230666441915
...
1.5462526341887264 0.9996988186962043 0.0245412285229123
1.5523885573402689 0.9998305817958233 0.0184067299058048
1.5585244804918115 0.9999247018391445 0.0122715382857199
1.5646604036433540 0.9999811752826011 0.0061358846491545
1.5707963267948966 1.0000000000000000 0.0000000000000000
请按任意键继续. . .

复制代码

G-Spider · 发表于 2011-2-23 10:50:42

12# wayne
#13实际结果表明，结果已经可以了。对于32位不需要这高的精度，double型也就只能表示15~16位的有效数，实在觉得过意不去，载入fldpi 得到扩展双精度的pi值，但现在已经够用了。

younger110 · 发表于 2011-9-30 10:38:04

牛人，学习了。

账号		自动登录	找回密码
密码			欢迎注册

[讨论] 正余弦的近似计算