正余弦的近似计算

G-Spider

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

1. #include <stdio.h>
2. #include <math.h>
4. #define PI 3.1415926535897932
5. #define _2PI (2.0*PI)
7. double lpsin[256+1],lpcos[256+1];
8. int main()
9. {
10. double N1,theta;
11. unsigned int i,j,tmpi,N,N0;
13. j=10;
14. N=(unsigned int)1<<j;
15. N0=N>>2; //256
16. N1=1.0/(double)N;
18. theta=N1*_2PI;
19. //lpsin=(double *)malloc( (N0+1)*sizeof(double) );
20. //lpcos=(double *)malloc( (N0+1)*sizeof(double) );
21. //if(lpsin==NULL || lpcos==NULL)
22. //printf("Memory Allocated Error\n");
23. lpsin[0]=0.0;
24. lpcos[0]=1.0;
25. lpsin[1]=sin(theta);
26. lpcos[1]=cos(theta);
28. for(i=2;i<=N0;i+=2)
29. {//i=2,4,6,....
30. tmpi=i>>1;
32. lpsin[i]= 2.0* lpsin[tmpi] * lpcos[tmpi];
33. lpcos[i]=1.0- 2.0* lpsin[tmpi] * lpsin[tmpi];
34. //======================================================
35. lpsin[i+1]= lpcos[i] * lpsin[1] + lpsin[i] * lpcos[1];
36. lpcos[i+1]=-lpsin[i] * lpsin[1] + lpcos[i] * lpcos[1];
37. }
39. for(i=0;i<=N0;i++)
40. printf("%18.14lf %18.14lf\n",*(lpsin+i),*(lpcos+i));
42. system("PAUSE");
43. return 0;
44. }

复制代码

wayne

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

G-Spider

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次浮点乘即可。

1. #include <stdio.h>
2. #include <math.h>
4. #define PI 3.1415926535897932
5. #define _2PI (2.0*PI)
7. double lpsin[256+1],lpcos[256+1];
8. int main()
9. {
10. double N1,theta;
11. unsigned int i,j,tmpi,N,N0;
13. j=10;
14. N=(unsigned int)1<<j; //N=1024
15. N0=N>>3; //N0=1024/8=128
16. N1=1.0/(double)N;
18. theta=N1*_2PI;
19. //lpsin=(double *)malloc( (N0+1)*sizeof(double) );
20. //lpcos=(double *)malloc( (N0+1)*sizeof(double) );
21. //if(lpsin==NULL || lpcos==NULL)
22. //printf("Memory Allocated Error\n");
23. lpsin[0]=0.0;
24. lpcos[0]=1.0;
25. lpsin[1]=sin(theta);
26. lpcos[1]=cos(theta);
28. for(i=2;i<=N0;i+=2)
29. {//i=2,4,6,....
30. tmpi=i>>1;
32. lpsin[i]= 2.0* lpsin[tmpi] * lpcos[tmpi];
33. lpcos[i]=1.0- 2.0* lpsin[tmpi] * lpsin[tmpi];
34. //======================================================
35. lpsin[i+1]= lpcos[i] * lpsin[1] + lpsin[i] * lpcos[1];
36. lpcos[i+1]=-lpsin[i] * lpsin[1] + lpcos[i] * lpcos[1];
37. }
39. N0<<=1;
40. for(;i<=N0;i+=2)
41. {//i=N0'+2,N0'+4,N0'+6,....,N0 (N0'=N0/2)
42. //sin(pi/2 -x)=cos(x); sin(N0-i)=cos(i);
43. //cos(pi/2-x)=sin(x); cos(N0-i)=sin(i);
44. tmpi=N0-i;
45. lpsin[i]=lpcos[tmpi];
46. lpcos[i]=lpsin[tmpi];
47. //======================================================
48. tmpi--;
49. lpsin[i+1]=lpcos[tmpi];
50. lpcos[i+1]=lpsin[tmpi];
51. }
53. printf(" rad sin(rad) cos(rad)\n");
54. for(i=0;i<=N0;i++)
55. printf("%18.16lf %18.16lf %18.16lf\n",i*N1*_2PI ,*(lpsin+i),*(lpcos+i));
57. system("PAUSE");
58. return 0;
59. }

复制代码

部分结果：

2. rad sin(rad) cos(rad)
3. 0.0000000000000000 0.0000000000000000 1.0000000000000000
4. 0.0061359231515426 0.0061358846491545 0.9999811752826011
5. 0.0122718463030851 0.0122715382857199 0.9999247018391445
6. 0.0184077694546277 0.0184067299058048 0.9998305817958233
7. 0.0245436926061703 0.0245412285229123 0.9996988186962043
8. 0.0306796157577128 0.0306748031766366 0.9995294175010933
9. 0.0368155389092554 0.0368072229413588 0.9993223845883495
10. 0.0429514620607980 0.0429382569349408 0.9990777277526455
11. 0.0490873852123405 0.0490676743274180 0.9987954562051724
12. 0.0552233083638831 0.0551952443496899 0.9984755805732947
13. 0.0613592315154256 0.0613207363022086 0.9981181129001492
14. 0.0674951546669682 0.0674439195636641 0.9977230666441915
16. ...
18. 1.5462526341887264 0.9996988186962043 0.0245412285229123
19. 1.5523885573402689 0.9998305817958233 0.0184067299058048
20. 1.5585244804918115 0.9999247018391445 0.0122715382857199
21. 1.5646604036433540 0.9999811752826011 0.0061358846491545
22. 1.5707963267948966 1.0000000000000000 0.0000000000000000
23. 请按任意键继续. . .

复制代码

G-Spider

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

younger110

牛人，学习了。