0xAA55的cos轮子

watermelon

前两天站长在群里发了一个是他自己写的cos函数，我一想这个真的是太好的机会了，大佬的源码肯定要拿来研读一下啊，肯定对增长功力是大补啊
于是我就偷偷的把一段函数从QQ中copy一下（CV大法好）,粘贴到vs中，然后根据群里的一些语言进行相应的修改，当天晚上就运行了一下，发现
编译不过。给0xAA55前辈发了一个消息，第二天一早，站长回了我的消息同时也是非常耐心地给我解释了很多，最终程序能够正常编译了，然后他给
了我一个网址，是Intel Intrinsics Guide，他给我说了查询的方法，我就开始按照上面的指令开始翻译。废话说的有点多了，开始正题。

1.简介
        0xAA55写的cos的原理是基于泰勒展开式，既然有了这个大的方向，那么翻译就会很有目的了哈哈。一开始感觉SSE的指令集非常牛x，第一次见完全
看不懂。不够也没有关系，好在Guide对于每一条指令都有description和operation，非常容易理解，然后我就翻译核心函数r_cos_sse翻译了三个小时
终于把泰勒展开的原理弄清楚了，SSE中的__m128这个类型是一次存储4个float单精度数（32位），从右往左分别是0,1,2,3位float数，很多指令的参数和返回值就是__m128，
这样很多指令就可以依次进行4个数的加减乘除运算，效率就会增加很多。至于r_cos_sse是怎么实现泰勒展开的，请看下面的源程序和小弟写的备注。

2.运行时候的问题
        在编译完程序的时候，站长说要我运行一下看看测试结果，运行时候发现math.h中的cos比r_cos_sse要快很多，站长给我说我的测试方法和编译器的优化参数有问题
并且不能用debug版，因为用debug版的话，每一条_mm指令都会写内存，所以要用release版，还有编译器的参数优化，除了启用内置函数关闭以外其他的都要开启。好的我照做。
怀着忐忑的心情在15：40去上课，
        晚上吃完了饭我就开始在宿舍反复运行那个用于测试的主函数，很是奇怪，我发现运行cos函数的时间是0，没错，一直是0。这个就很郁闷了。我就开始在网上查VS的
参数优化，发现没有什么特别的。无奈把那个release版本的程序拖进IDA pro看看，发现它居然在release版本下没有进行cos的循环，反汇编的结果只有一个测试r_cos_sse的双重循环
这个就是问题了，我上网一查，发现vs的release版本在编译运行程序时候会自动把没有用处的东西给忽略掉。我天，我就给math.h的cos函数的测试写成float test = cos((float)j);
最后还得再printf("%f\n",test);否则你会发现release还是会略去cos的测试的双重循环。然后程序果然照预期进行了，测试的结果是，二者分别在100*1048576次运算的条件下，r_cos_sse
的计算速度几乎是math.h中cos函数计算速度的2.2倍。并且在开始的r_cos_sse的宏定义可以在IDA中看出，他是按照宏的定义，在调用时直接把代码放在了调用区域，在这方面可以说对运行速度也是有贡献的。

3.感想
        非常高兴可以翻译0xAA55大佬的核心函数和在他的指导下完成的理解，同时在编写测试主函数的时候也是给了我大量的提示，我也是从中接触到了SSE指令集和他的优越性，理解vs的编译时优化参数的设置，
接触到了volatile关键字（这个关键字在0xAA55的前两天的关于自旋锁的设计中有提及），更重要的是解决问题的能力。好东西不能一个人分享，同时我也不能过于自私，所以我在征得站长的同意后
决定发出这篇帖子，小弟水平有限，有些语句可能理解不是很到位，希望各位大佬可以直接给小弟指正，不必留情。最后小弟用python画了一个简单的二维的对比图，很是简陋，请见谅。

1. //此程序核心函数r_cos_sse由0xAA55编写
   
2. #include <xmmintrin.h>                        //SSE指令集头文件，此头文件里包含MMX头文件
   
3. #include <stdio.h>
   
4. #include <stdlib.h>
   
5. #include <math.h>
   
6. #include <windows.h>
   
7. 
   
8. typedef float real_t;                        //将float类型定义为real_t
   
9. #define r_pi 3.1415926                        //定义圆周率π为3.1415926
   
10. #define mathsimd_func(r,n) r n        //宏定义
    
11. 
    
12. //这里进行了宏替换，直接将mathsimd_func(r,n)替换为r n
    
13. //即real_t r_cos_sse(real_t x){...} 返回值为real_t,函数名为r_cos_sse,形参为real_t x
    
14. mathsimd_func(real_t, r_cos_sse)(real_t x)
    
15. {
    
16.         //__m128 类似于4个float，占128bit
    
17.         //后面就会看到mfac2_4_6_8是对应x乘方的系数的分母
    
18.         __m128 mfac2_4_6_8 = _mm_set_ps(                //将下面四个单精度的浮点值赋值给mfac2_4_6_8
    
19.                 -2.0f,
    
20.                 24.0f,
    
21.                 -720.0f,
    
22.                 40320.0f);
    
23.         //后面就会看到mfac2_4_6_8是对应x乘方的系数
    
24.         __m128 mfac10_12_14_16 = _mm_set_ps(        //将下面四个单精度的浮点值赋值给mfac2_4_6_8
    
25.                 -3628800.0f,
    
26.                 479001600.0f,
    
27.                 -87178291200.0f,
    
28.                 20922789888000.0f);
    
29.         __m128 m1111 = _mm_set1_ps(1);                        //给m1111赋值四个单精度浮点数1
    
30.         __m128 m111x;                                                        //一些变量的声明
    
31.         __m128 m11xx;
    
32.         __m128 m1xxx;
    
33.         __m128 mxxxx;
    
34.         __m128 mxp2_4_6_8;
    
35.         __m128 mxp10_12_14_16;
    
36.         __m128 mxp8_8_8_8;
    
37.         __m128 m1234;
    
38.         __m128 mr;
    
39.         real_t r;
    
40. 
    
41.         //0xAA55语：这一句就是单纯地为了提高精度而写的，因为泰勒级数展开的精度取决于算的次数和X的绝对值
    
42.         //所以我让X的值“折返”来利用cos的对称性
    
43.         x = x - (float)floor(x / (r_pi * 2)) * r_pi * 2;
    
44.         if (x > r_pi) x = r_pi * 2 - x;
    
45. 
    
46.         mxxxx = _mm_load1_ps(&x);                                        //将变化过后的4个x值写入
    
47.         m111x = _mm_move_ss(m1111, mxxxx);                        //将mxxxx中的最低位的x移入m111x，将m1111中的高三位float移入m111x，所以变量取名为m111x（真形象，:D)
    
48.         m11xx = _mm_movelh_ps(mxxxx, m1111);                //顾名思义
    
49.         m1xxx = _mm_movelh_ps(mxxxx, m111x);                //顾名思义
    
50. 
    
51. 
    
52.         //将mxxxx和m111x中对应位置的数进行相乘，每一位的结果放入mxp2_4_6_8中
    
53.         //以下_mm_mul_ps函数同理,通过查询说明也可以得知                                                mxp_2_4_6_8中的每一位
    
54.         //下面四句程序用纸和笔来写一下就可以清晰看出        //                                              2    4    6    8
    
55.         mxp2_4_6_8 = _mm_mul_ps(mxxxx, m111x);                        //                                      x    x    x    x²       
    
56.         mxp2_4_6_8 = _mm_mul_ps(mxp2_4_6_8, m11xx);                //                                      x    x    x²   x³
    
57.         mxp2_4_6_8 = _mm_mul_ps(mxp2_4_6_8, m1xxx);                //                                      x    x²   x³   x_4
    
58.         mxp2_4_6_8 = _mm_mul_ps(mxp2_4_6_8, mxp2_4_6_8);  //   x_4,x_6,x_8代表次方                    x²  x_4   x_6  x_8
    
59. 
    
60.         //利用vs中的查看定义可以得出_MM_SHUFFLE的含义,再查询说明，可以得出mxp8_8_8_8中每一项就是mxp2_4_6_8中的8所对应的那一项
    
61.         mxp8_8_8_8 = _mm_shuffle_ps(mxp2_4_6_8, mxp2_4_6_8, _MM_SHUFFLE(0, 0, 0, 0));
    
62. 
    
63.         mxp10_12_14_16 = _mm_mul_ps(mxp2_4_6_8, mxp8_8_8_8);//计算精度到了x的16次方除以16的阶乘那一项
    
64. 
    
65.         //下面一句就是-x²/2! + x_4/4! - x_6/6! + x_8/8! .....
    
66.         //更加准确地格式上是：这一句其实是让-x²/2!和-x¹⁰/10!相加，让x⁴/4!和x¹²/12!相加，让-x⁶/6!和-x¹⁴/14!相加，让x⁸/8!和x¹⁶/16!相加，然后存入m1234里面。
    
67.         //之后再想办法吧m1234的各项都加起来就得到结果了。
    
68.         m1234 = _mm_add_ps(_mm_div_ps(mxp2_4_6_8, mfac2_4_6_8),
    
69.                 _mm_div_ps(mxp10_12_14_16, mfac10_12_14_16));
    
70. 
    
71.         //同理_MM_SHUFFLE(2,3,0,1)算出来然后再进行_mm_shuffle_ps(m1234,m1234,imm8)就是取出来m1234的第二项，第三项，第0项，第一项（从右向左数）的每一项的组合
    
72.         //_mm_add_ps(m1234, _mm_shuffle_ps(m1234, m1234, _MM_SHUFFLE(2, 3, 0, 1)));
    
73.         //这句让1和2相加，3和4相加。之后再考虑把3和4相加的结果再加到一起。
    
74.         mr = _mm_add_ps(m1234, _mm_shuffle_ps(m1234, m1234, _MM_SHUFFLE(2, 3, 0, 1)));
    
75. 
    
76.         //下面一句就是cos(x) = 1-x²/2! + x_4/4! - x_6/6! + x_8/8! .....的泰勒展开的式子，mr就是结果
    
77.         mr = _mm_add_ss(m1111, _mm_add_ss(mr, _mm_shuffle_ps(mr, mr, _MM_SHUFFLE(1, 0, 3, 2))));
    
78. 
    
79.         _mm_store_ss(&r, mr);                                //0xAA55语：_mm_store_ss这条指令相当于告诉CPU：抽空把数值往内存里存一下。
    
80.         _mm_sfence();                                                //_mm_sfence则相当于告诉CPU：存完了再继续运行,类似于文件操作中的fflush(stdio.h中的函数)
    
81.         //_compiler_barrier;                                //源程序中原有指令
    
82.         //asm volatile ("" ::: "memory");        //gcc
    
83.         _ReadWriteBarrier();                //0xAA55语：vc中末尾的_ReadWriteBarrier的作用是保证编译器一定先生成_mm_sfence();的对应指令，再生成return r;的对应指令
    
84.         //否则算出来的数值还没来得及存入r呢，它就从r这个变量取出数值用来返回了。
    
85.         return r;                                                        //返回结果
    
86. }
    
87. 
    
88. 
    
89. 
    
90. int main(void)
    
91. {
    
92.         float test = 0;                                        //存放每次结果，防止math.h中的cos被release优化掉
    
93.         SetThreadAffinityMask(GetCurrentThread(), 1);        //用来设置特定的线程与哪几个CPU核心相关
    
94.         //进行0xAA55的r_cos_sse
    
95.         LARGE_INTEGER begin_r_cos_sse;
    
96.         LARGE_INTEGER end_r_cos_sse;
    
97.         LARGE_INTEGER freq_r_cos_sse;
    
98.         QueryPerformanceFrequency(&freq_r_cos_sse);
    
99.         QueryPerformanceCounter(&begin_r_cos_sse);
    
100. 
     
101.         for (int i = 0; i < 100; i++)                        //外循环100次
     
102.         {
     
103.                 for (int j = 0; j < 1048576; j++)        //内循环进行1048576次运算
     
104.                 {
     
105.                         test = r_cos_sse((float)j);       
     
106.                 }
     
107.         }
     
108. 
     
109.         QueryPerformanceCounter(&end_r_cos_sse);
     
110.         printf("r_cos_sse中的test:%f\n", test);
     
111.         printf("运行r_cos_sse时间:%f秒\n", (float)(end_r_cos_sse.QuadPart - begin_r_cos_sse.QuadPart) / (float)freq_r_cos_sse.QuadPart);
     
112. 
     
113.         printf("*****************************************\n");                        //分割线
     
114. 
     
115.         //进行math.h中cos的测试
     
116.         LARGE_INTEGER begin_cos;
     
117.         LARGE_INTEGER end_cos;
     
118.         LARGE_INTEGER freq_cos;
     
119.         QueryPerformanceFrequency(&freq_cos);
     
120.         QueryPerformanceCounter(&begin_cos);
     
121.         for (int i = 0; i < 100; i++)
     
122.         {
     
123.                 for (int j = 0; j < 1048576; j++)        //内循环进行1048576次运算
     
124.                 {
     
125.                         test = cos((float)j);
     
126.                         //cos((float)j);
     
127.                 }
     
128.         }
     
129.         QueryPerformanceCounter(&end_cos);
     
130.         printf("cos中的test:%f\n", test);
     
131.         printf("运行cos时间:%f秒\n", (float)(end_cos.QuadPart - begin_cos.QuadPart) / (float)(freq_cos.QuadPart));
     
132. 
     
133.         system("pause");
     
134.         return 0;
     
135. }
     
136. 
     
137. 
     

复制代码

运行时间对比图：横轴为测试的次数，一共测试了10次，纵轴为时间，单位秒




Intel Intrinsics Guide:  https://software.intel.com/sites ... =SSE&expand=825
工程和测试结果和图片以及py脚本在附件中：

0xAA55

//下面一句就是-x²/2! + x_4/4! - x_6/6! + x_8/8! .....
m1234 = _mm_add_ps(_mm_div_ps(mxp2_4_6_8, mfac2_4_6_8),
    _mm_div_ps(mxp10_12_14_16, mfac10_12_14_16));

这一句其实是让-x²/2!和-x¹⁰/10!相加，让x⁴/4!和x¹²/12!相加，让-x⁶/6!和-x¹⁴/14!相加，让x⁸/8!和x¹⁶/16!相加，然后存入m1234里面。

之后再想办法吧m1234的各项都加起来就得到结果了。

_mm_add_ps(m1234, _mm_shuffle_ps(m1234, m1234, _MM_SHUFFLE(2, 3, 0, 1)));
这句让1和2相加，3和4相加。之后再考虑把3和4相加的结果再加到一起。

watermelon

0xAA55 发表于 2018-10-12 21:08
这一句其实是让-x²/2!和-x¹⁰/10!相加，让x⁴/4!和x¹²/12!相加，让-x⁶/6!和-x¹⁴/14!相加，让x⁸/8! ...

哦哦，谢谢站长，我在重新算算

watermelon

0xAA55 发表于 2018-10-12 21:08
这一句其实是让-x²/2!和-x¹⁰/10!相加，让x⁴/4!和x¹²/12!相加，让-x⁶/6!和-x¹⁴/14!相加，让x⁸/8! ...

真的，其中_mm_add_ps(m1234,_mm_shuffle_ps(m1234,m1234,_MM_SHUFFLE(2,3,0,1)));
就是_mm_add_ps(m1234,m2143),果然是的。

第一句那个我现在明白没有说清，话说站长真的是严谨（严格）啊
我马上改一下再。

0xAA55

SIMD的优化在于它可以让CPU用一条单一的指令对多个数据进行同时的计算，典型的就如代码里面的例子，一条指令对4个浮点数做同时的加、减、乘、除等运算。
而不直接写汇编的原因则在于，使用C语法风格的intrinsic既可以让编译器自动选择合适的寄存器，又可以让编译器能有一定的条件进行自动的合理的优化，比如函数自动内联等。
而且，代码里面的__m128这种变量类型它其实就是代表了一个SSE寄存器（xmm0-xmm7中的一个），它只在debug的时候出于调试的目的才会被写入到内存——你可以在VS等调试器里面通过监视等方式查看xmm寄存器理论上的数值。它在release的情况下是直接操作xmm0-xmm7寄存器的。

此外，VS编译器的“使用内置函数”优化，可以让cos()的行为变为自动使用SSE2的指令的方式——通过比较一个名为__use_sse2_mathfcns的全局变量来判断CPU是否支持SSE2，从而跳转到使用该指令集的对应位置继续执行，而非只使用FPU指令集的fcos指令进行余弦的计算。sse2与sse相比主要的差异在于sse2提供了大量的针对双double变量类型二维向量的加速指令集，而非sse的针对4float变量类型四维向量的加速指令集，存在精度上的差异。在不追求精度只追求效率、并且不想每次都判断__use_sse2_mathfcns的值后再进行计算的话，可以直接做个cos函数的函数指针，然后根据CPUID的查询情况自动选择使用cos（封装为输入float返回float）或r_cos_sse中的一个的方式，来达到优化的效果。