技术宅的结界

 找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 177|回复: 4
收起左侧

0xAA55的cos轮子

[复制链接]

11

主题

100

帖子

738

积分

用户组: 大·技术宅

UID
3808
精华
1
威望
16 点
宅币
557 个
贡献
44 次
宅之契约
0 份
在线时间
90 小时
注册时间
2018-5-6
发表于 2018-10-10 23:12:49 | 显示全部楼层 |阅读模式

欢迎访问技术宅的结界,请注册或者登录吧。

您需要 登录 才可以下载或查看,没有帐号?立即注册→加入我们

x
本帖最后由 watermelon 于 2018-10-12 22:04 编辑

前两天站长在群里发了一个是他自己写的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画了一个简单的二维的对比图,很是简陋,请见谅。

[C] 纯文本查看 复制代码
//此程序核心函数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_6  x_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 ... =SSE&expand=825
工程和测试结果和图片以及py脚本在附件中:


对比图.png

0xAA55的re_cos_sse.zip

131.26 KB, 下载次数: 2

评分

参与人数 1威望 +2 宅币 +16 贡献 +4 收起 理由
0xAA55 + 2 + 16 + 4 很棒

查看全部评分

菜鸟一枚,直接指正,不必留情

1007

主题

2232

帖子

5万

积分

用户组: 管理员

一只技术宅

UID
1
精华
199
威望
263 点
宅币
16684 个
贡献
33359 次
宅之契约
0 份
在线时间
1595 小时
注册时间
2014-1-26
发表于 2018-10-12 21:08:41 | 显示全部楼层
//下面一句就是-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相加的结果再加到一起。

11

主题

100

帖子

738

积分

用户组: 大·技术宅

UID
3808
精华
1
威望
16 点
宅币
557 个
贡献
44 次
宅之契约
0 份
在线时间
90 小时
注册时间
2018-5-6
 楼主| 发表于 2018-10-12 21:36:11 | 显示全部楼层
0xAA55 发表于 2018-10-12 21:08
这一句其实是让-x²/2!和-x¹⁰/10!相加,让x⁴/4!和x¹²/12!相加,让-x⁶/6!和-x¹⁴/14!相加,让x⁸/8! ...

哦哦,谢谢站长,我在重新算算
菜鸟一枚,直接指正,不必留情

11

主题

100

帖子

738

积分

用户组: 大·技术宅

UID
3808
精华
1
威望
16 点
宅币
557 个
贡献
44 次
宅之契约
0 份
在线时间
90 小时
注册时间
2018-5-6
 楼主| 发表于 2018-10-12 21:56:52 | 显示全部楼层
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),果然是的。

第一句那个我现在明白没有说清,话说站长真的是严谨(严格)啊
我马上改一下再。
菜鸟一枚,直接指正,不必留情

1007

主题

2232

帖子

5万

积分

用户组: 管理员

一只技术宅

UID
1
精华
199
威望
263 点
宅币
16684 个
贡献
33359 次
宅之契约
0 份
在线时间
1595 小时
注册时间
2014-1-26
发表于 2018-10-16 01:47:50 | 显示全部楼层
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中的一个的方式,来达到优化的效果。

本版积分规则

QQ|申请友链|Archiver|手机版|小黑屋|技术宅的结界 ( 滇ICP备16008837号|网站地图

GMT+8, 2018-11-20 13:41 , Processed in 0.089908 second(s), 16 queries , Gzip On, Memcache On.

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表