0xAA55的cos轮子
本帖最后由 watermelon 于 2018-10-12 22:04 编辑前两天站长在群里发了一个是他自己写的cos函数，我一想这个真的是太好的机会了，大佬的源码肯定要拿来研读一下啊，肯定对增长功力是大补啊
于是我就偷偷的把一段函数从QQ中copy一下（CV大法好:D）,粘贴到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画了一个简单的二维的对比图，很是简陋，请见谅。
//此程序核心函数r_cos_sse由0xAA55编写
#include <xmmintrin.h> //SSE指令集头文件，此头文件里包含MMX头文件
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
typedef float real_t; //将float类型定义为real_t
#define r_pi 3.1415926 //定义圆周率π为3.1415926
#define mathsimd_func(r,n) r n //宏定义
//这里进行了宏替换，直接将mathsimd_func(r,n)替换为r n
//即real_t r_cos_sse(real_t x){...} 返回值为real_t,函数名为r_cos_sse,形参为real_t x
mathsimd_func(real_t, r_cos_sse)(real_t x)
{
//__m128 类似于4个float，占128bit
//后面就会看到mfac2_4_6_8是对应x乘方的系数的分母
__m128 mfac2_4_6_8 = _mm_set_ps( //将下面四个单精度的浮点值赋值给mfac2_4_6_8
-2.0f,
24.0f,
-720.0f,
40320.0f);
//后面就会看到mfac2_4_6_8是对应x乘方的系数
__m128 mfac10_12_14_16 = _mm_set_ps( //将下面四个单精度的浮点值赋值给mfac2_4_6_8
-3628800.0f,
479001600.0f,
-87178291200.0f,
20922789888000.0f);
__m128 m1111 = _mm_set1_ps(1); //给m1111赋值四个单精度浮点数1
__m128 m111x; //一些变量的声明
__m128 m11xx;
__m128 m1xxx;
__m128 mxxxx;
__m128 mxp2_4_6_8;
__m128 mxp10_12_14_16;
__m128 mxp8_8_8_8;
__m128 m1234;
__m128 mr;
real_t r;
//0xAA55语：这一句就是单纯地为了提高精度而写的，因为泰勒级数展开的精度取决于算的次数和X的绝对值
//所以我让X的值“折返”来利用cos的对称性
x = x - (float)floor(x / (r_pi * 2)) * r_pi * 2;
if (x > r_pi) x = r_pi * 2 - x;
mxxxx = _mm_load1_ps(&x); //将变化过后的4个x值写入
m111x = _mm_move_ss(m1111, mxxxx); //将mxxxx中的最低位的x移入m111x，将m1111中的高三位float移入m111x，所以变量取名为m111x（真形象，:D)
m11xx = _mm_movelh_ps(mxxxx, m1111); //顾名思义
m1xxx = _mm_movelh_ps(mxxxx, m111x); //顾名思义
//将mxxxx和m111x中对应位置的数进行相乘，每一位的结果放入mxp2_4_6_8中
//以下_mm_mul_ps函数同理,通过查询说明也可以得知 mxp_2_4_6_8中的每一位
//下面四句程序用纸和笔来写一下就可以清晰看出 // 2 4 6 8
mxp2_4_6_8 = _mm_mul_ps(mxxxx, m111x); // x x x x²
mxp2_4_6_8 = _mm_mul_ps(mxp2_4_6_8, m11xx); // x x x² x³
mxp2_4_6_8 = _mm_mul_ps(mxp2_4_6_8, m1xxx); // x x² x³ x_4
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_6x_8
//利用vs中的查看定义可以得出_MM_SHUFFLE的含义,再查询说明，可以得出mxp8_8_8_8中每一项就是mxp2_4_6_8中的8所对应的那一项
mxp8_8_8_8 = _mm_shuffle_ps(mxp2_4_6_8, mxp2_4_6_8, _MM_SHUFFLE(0, 0, 0, 0));
mxp10_12_14_16 = _mm_mul_ps(mxp2_4_6_8, mxp8_8_8_8);//计算精度到了x的16次方除以16的阶乘那一项
//下面一句就是-x²/2! + x_4/4! - x_6/6! + x_8/8! .....
//更加准确地格式上是：这一句其实是让-x²/2!和-x¹⁰/10!相加，让x⁴/4!和x¹²/12!相加，让-x⁶/6!和-x¹⁴/14!相加，让x⁸/8!和x¹⁶/16!相加，然后存入m1234里面。
//之后再想办法吧m1234的各项都加起来就得到结果了。
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));
//同理_MM_SHUFFLE(2,3,0,1)算出来然后再进行_mm_shuffle_ps(m1234,m1234,imm8)就是取出来m1234的第二项，第三项，第0项，第一项（从右向左数）的每一项的组合
//_mm_add_ps(m1234, _mm_shuffle_ps(m1234, m1234, _MM_SHUFFLE(2, 3, 0, 1)));
//这句让1和2相加，3和4相加。之后再考虑把3和4相加的结果再加到一起。
mr = _mm_add_ps(m1234, _mm_shuffle_ps(m1234, m1234, _MM_SHUFFLE(2, 3, 0, 1)));
//下面一句就是cos(x) = 1-x²/2! + x_4/4! - x_6/6! + x_8/8! .....的泰勒展开的式子，mr就是结果
mr = _mm_add_ss(m1111, _mm_add_ss(mr, _mm_shuffle_ps(mr, mr, _MM_SHUFFLE(1, 0, 3, 2))));
_mm_store_ss(&r, mr); //0xAA55语：_mm_store_ss这条指令相当于告诉CPU：抽空把数值往内存里存一下。
_mm_sfence(); //_mm_sfence则相当于告诉CPU：存完了再继续运行,类似于文件操作中的fflush(stdio.h中的函数)
//_compiler_barrier; //源程序中原有指令
//asm volatile ("" ::: "memory"); //gcc
_ReadWriteBarrier(); //0xAA55语：vc中末尾的_ReadWriteBarrier的作用是保证编译器一定先生成_mm_sfence();的对应指令，再生成return r;的对应指令
//否则算出来的数值还没来得及存入r呢，它就从r这个变量取出数值用来返回了。
return r; //返回结果
}
int main(void)
{
float test = 0; //存放每次结果，防止math.h中的cos被release优化掉
SetThreadAffinityMask(GetCurrentThread(), 1); //用来设置特定的线程与哪几个CPU核心相关
//进行0xAA55的r_cos_sse
LARGE_INTEGER begin_r_cos_sse;
LARGE_INTEGER end_r_cos_sse;
LARGE_INTEGER freq_r_cos_sse;
QueryPerformanceFrequency(&freq_r_cos_sse);
QueryPerformanceCounter(&begin_r_cos_sse);
for (int i = 0; i < 100; i++) //外循环100次
{
for (int j = 0; j < 1048576; j++) //内循环进行1048576次运算
{
test = r_cos_sse((float)j);
}
}
QueryPerformanceCounter(&end_r_cos_sse);
printf("r_cos_sse中的test:%f\n", test);
printf("运行r_cos_sse时间:%f秒\n", (float)(end_r_cos_sse.QuadPart - begin_r_cos_sse.QuadPart) / (float)freq_r_cos_sse.QuadPart);
printf("*****************************************\n"); //分割线
//进行math.h中cos的测试
LARGE_INTEGER begin_cos;
LARGE_INTEGER end_cos;
LARGE_INTEGER freq_cos;
QueryPerformanceFrequency(&freq_cos);
QueryPerformanceCounter(&begin_cos);
for (int i = 0; i < 100; i++)
{
for (int j = 0; j < 1048576; j++) //内循环进行1048576次运算
{
test = cos((float)j);
//cos((float)j);
}
}
QueryPerformanceCounter(&end_cos);
printf("cos中的test:%f\n", test);
printf("运行cos时间:%f秒\n", (float)(end_cos.QuadPart - begin_cos.QuadPart) / (float)(freq_cos.QuadPart));
system("pause");
return 0;
}
运行时间对比图：横轴为测试的次数，一共测试了10次，纵轴为时间，单位秒
Intel Intrinsics Guide:https://software.intel.com/sites/landingpage/IntrinsicsGuide/#techs=SSE&expand=825
工程和测试结果和图片以及py脚本在附件中：
//下面一句就是-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相加的结果再加到一起。 0xAA55 发表于 2018-10-12 21:08
这一句其实是让-x²/2!和-x¹⁰/10!相加，让x⁴/4!和x¹²/12!相加，让-x⁶/6!和-x¹⁴/14!相加，让x⁸/8! ...
哦哦，谢谢站长，我在重新算算 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),果然是的。
第一句那个我现在明白没有说清，话说站长真的是严谨（严格）啊:handshake
我马上改一下再。 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中的一个的方式，来达到优化的效果。
页:
[1]