SIMD函数整理：01 《PC平台新技术MMX（上册）：开发编程指南》第8章 MMX编码技术

最新推荐文章于 2019-12-19 15:14:19 发布
zyl910
最新推荐文章于 2019-12-19 15:14:19 发布
阅读量2.4k
点赞数
分类专栏：数学&优化技巧 #All - 所有 SIMD 文章标签：编程平台 dst matrix 扩展互联网
本文链接：https://blog.csdn.net/zyl910/article/details/7516127
版权
#All - 所有同时被 3 个专栏收录
61 篇文章 1 订阅
订阅专栏
数学&优化技巧
24 篇文章 1 订阅
订阅专栏
SIMD
13 篇文章 5 订阅
订阅专栏
一、来源
　　来源：《PC平台新技术MMX（上册）：开发编程指南》第8章 MMX编码技术
　　书籍信息——
http://e.360buy.com/30027396.html
PC平台新技术MMX（上册）：开发编程指南
作者：吴乐南编
出版社：东南大学出版社
ＩＳＢＮ：9787810502528
出版时间：1997-10-01
页数：149
字数：237000
所属分类：
电子书 > 计算机与互联网 > 编程语言与程序设计
电子书 > 计算机与互联网 > 计算机工具书
二、整理后的代码
　　代码——
#include <Windows.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <conio.h>
#include <assert.h>

// MMX, SSE, SSE2
#include <emmintrin.h>

// 紧缩无符号字 解包为 两组紧缩无符号双字
// 章节：8.1 数据拆封/8.1.1 无符号数拆封
//
// result: 两个零扩展的32位双字，来自源的两个低端字。
// mm1_dst_hi: 两个零扩展的32位双字，来自源的两个高端字。
// mm0_src: 源值（紧缩16位无符号数）。
inline __m64 md_unpack_mud4muw(__m64& mm1_dst_hi, const __m64 mm0_src)
{
	__m64 muwZero = _mm_setzero_si64();	// [MMX]赋值为0
	mm1_dst_hi = _mm_unpackhi_pi16(mm0_src, muwZero);	// 把两个高端字拆封到两个32位双字中。[MMX]高位解包.字到双字
	return       _mm_unpacklo_pi16(mm0_src, muwZero);	// 把两个低端字拆封到两个32位双字中。[MMX]低位解包.字到双字
}

// 紧缩带符号字 解包为 两组紧缩带符号双字
// 章节：8.1 数据拆封/8.1.2 带符号数拆封
//
// result: 两个符号扩展的32位双字，来自源的两个低端字。
// mm1_dst_hi: 两个符号扩展的32位双字，来自源的两个高端字。
// mm0_src: 源值（紧缩16位带符号数）。
inline __m64 md_unpack_mid4miw(__m64& mm1_dst_hi, const __m64 mm0_src)
{
	// 注：其实并不需要读取mm1_dst_hi，但为了符合语法，只能这样写。
	mm1_dst_hi = _mm_srai_pi32(_mm_unpackhi_pi16(mm1_dst_hi, mm0_src), 16); // 把源数据的两个高端字拆分到 第1字与第3字（即两个紧缩双字的高16位），再紧缩双字算术右移16位。使源数据的两个高端字扩展为2个32位带符号双字。
	return       _mm_srai_pi32(_mm_unpacklo_pi16(mm0_src,    mm0_src), 16); // 把源数据的两个低端字拆分到 第1字与第3字（即两个紧缩双字的高16位），再紧缩双字算术右移16位。使源数据的两个低端字扩展为2个32位带符号双字。
}

// 两组紧缩带符号双字 交叉饱和紧缩为 紧缩带符号字
// 章节：8.2 数据紧缩/8.2.1 带饱和的交叉紧缩
// 例如：将 {[B1,B0], [A1,A0]} 交叉紧缩为 {[B1',A1',B0',A0']}
// 注：紧缩（_mm_packs_pi32）是将 {[B1,B0], [A1,A0]} 转为 {[B1',B0',A1',A0']}
//
// result: 紧缩16位带符号数。第0字和第2字来自mm0_lo的带符号饱和双字，第1字和第3字来自mm1_hi的带符号饱和双字。
// mm0_lo: 低位源值（A）。
// mm1_hi: 高位源值（B）。
inline __m64 md_pack_s_cross_miw4mid(__m64 mm0_lo, __m64 mm1_hi)
{
	mm1_hi = _mm_packs_pi32(mm1_hi, mm1_hi);	// 紧缩并且符号饱和。即变为[B1',B0',B1',B0']。[MMX]饱和打包.双字到字
	mm0_lo = _mm_packs_pi32(mm0_lo, mm0_lo);	// 紧缩并且符号饱和。即变为[A1',A0',A1',A0']。
	return _mm_unpacklo_pi16(mm0_lo, mm1_hi);	// 交叉操作数的低16位。[MMX]低位解包.字到双字
}

// 两组紧缩无符号双字 交叉环绕紧缩为 紧缩无符号字
// 章节：8.2 数据紧缩/8.2.2 不带饱和的交叉紧缩
// 例如：将 {[B1,B0], [A1,A0]} 交叉紧缩为 {[B1',A1',B0',A0']}
//
// result: 紧缩16位无符号数。第0字和第2字来自mm0_lo的无符号双字，第1字和第3字来自mm1_hi的无符号双字。
// mm0_lo: 低位源值（A）。
// mm1_hi: 高位源值（B）。
inline __m64 md_pack_w_cross_muw4mud(__m64 mm0_lo, __m64 mm1_hi)
{
	mm1_hi = _mm_slli_pi32(mm1_hi, 16);	// 将每个双字的低16位左移至高16位
	mm0_lo = _mm_and_si64(mm0_lo, _mm_set_pi16(0, (short)0xFFFF, 0, (short)0xFFFF));	// 用0屏蔽每个双字的最高16位
	return _mm_or_si64(mm0_lo, mm1_hi);	// 合并两个操作数
}

// 2x2矩阵转置.紧缩双字
// 章节：8.3 非交叉拆分
// 例如：将2x2矩阵 [[A1,A0] [B1,B0]] 转置为 [[B0,A0] [B1,A1]]。
// 
// [A1 A0]    [B0 A0]
// [B1 B0] -> [B1 A1]
// msb<-lsb
//
// mm0_row0: 2x2矩阵的第0行（A）。
// mm1_row1: 2x2矩阵的第1行（B）。
inline void md_matrix_transpose_2x2_mmd(__m64& mm0_row0, __m64& mm1_row1)
{
	__m64 tmp = mm0_row0;	// 备份第0行
	mm0_row0 = _mm_unpacklo_pi32(mm0_row0, mm1_row1);	// 高32位为mm1_row1的低32位（B0），低32位为源mm0_row0的低32位（A0）。[MMX]低位解包.双字到四字
	mm1_row1 = _mm_unpackhi_pi32(tmp     , mm1_row1);	// 高32位为mm1_row1的高32位（B1），低32位为源mm0_row0的高32位（A1）。[MMX]高位解包.双字到四字
}

// 复数与常量相乘（紧缩字->紧缩双字）
// 章节：8.4 复数与常量相乘
//
// result: 复数乘法的结果，高32位是实部，低32位是虚部。
// mm0_src: 被乘数（[?,?,Dr,Di]）。
// mm1_c: 已调整好顺序的常量乘数（[Cr,-Ci,Ci,Cr]）。
inline __m64 md_complex_mul_c_mid4miw(__m64 mm0_src, const __m64 mm1_c)
{
	mm0_src = _mm_unpacklo_pi32(mm0_src, mm0_src);	// 产生 [Dr,Di,Dr,Di]。[MMX]低位解包.双字到四字
	return _mm_madd_pi16(mm0_src, mm1_c);	// 操作结果是 [(Dr*Cr-Di*Ci), (Dr*Ci+Di*Cr)]。[MMX]乘后二加.带符号16位至带符号32位
}

// 无符号紧缩字节的绝对差
// 章节：8.5 数的绝对差\8.5.1 无符号数的绝对差
//
// result: 无符号紧缩字节的绝对差。伪代码——result[i]=abs(mm0[i] - mm1[i])。
// mm0: 源操作数A。
// mm1: 源操作数B。
inline __m64 md_absolute_deviation_mub(const __m64 mm0, const __m64 mm1)
{
	return _mm_or_si64(_mm_subs_pu8(mm0, mm1), _mm_subs_pu8(mm1, mm0));
		// 1. "_mm_subs_pu8(mm0, mm1)": 计算差值
		// 2. "_mm_subs_pu8(mm1, mm0)": 以另一种途径计算差值
		// 3. "_mm_or_si64(...,  ...)": 合并结果
}

// 带符号紧缩字的绝对差
// 章节：8.5 数的绝对差\8.5.2 带符号数的绝对差
//
// result: 带符号紧缩字的绝对差。伪代码——result[i]=abs(mm0[i] - mm1[i])。
// mm0: 源操作数A。
// mm1: 源操作数B。
inline __m64 md_absolute_deviation_miw(const __m64 mm0, const __m64 mm1)
{
	__m64 miwMaskGt = _mm_cmpgt_pi16(mm0, mm1);	// 产生 A>B 的屏蔽值
	__m64 miwXor = _mm_and_si64(_mm_xor_si64(mm0, mm1), miwMaskGt);	// 产生交换屏蔽值（仅在A>B时的XOR(A,B)值）。即当A>B时，该字是XOR(A,B)；而A<=B时，该字是是0。
	__m64 miwMin = _mm_xor_si64(mm0, miwXor);	// 当A>B时就用xor交换，产生最小值
	__m64 miwMax = _mm_xor_si64(mm1, miwXor);	// 当B<=A时就用xor交换，产生最大值
	return _mm_sub_pi16(miwMax, miwMin);	// 绝对差 = 最大值 - 最小值
}

// 带符号紧缩字的绝对值
// 章节：8.6 绝对值
//
// result: 带符号紧缩字的绝对值。伪代码——result[i]=abs(mm0[i])。
// mm0: 源操作数。
inline __m64 md_abs_miw(const __m64 mm0)
{
	__m64 miwSign = _mm_srai_pi16(mm0, 15);	// 将符号位转为掩码。使每个字为全0（对于非负数）或全1（对于负数）。注：补码下的“全1”代表数值“-1”，减法碰到“-1”就形成了“加一”。
	return _mm_subs_pi16(_mm_xor_si64(mm0, miwSign), miwSign);	// 为了获得绝对值，仅对负数求相反数。补码求相反数规则——原码取反再加一。
}

// 将带符号紧缩字限制在[iLow,iHigh]区间
// 章节：8.7 数值的截取/8.7.1 对任意有符号数范围截取符号数/[0]
//
// result: 限制后的带符号紧缩字。伪代码——result[i]=(mm0[i]<iLow)?iLow:( (mm0[i]>iHigh)?iHigh:mm0[i] )。
// mm0: 源操作数。
inline __m64 md_clamp_miw(const __m64 mm0, short iLow, short iHigh)
{
	const __m64 miwMinInt16 = _mm_set1_pi16((short)0x8000);	// 带符号16位的最小值
	__m64 tmp = _mm_add_pi16(mm0, miwMinInt16);	// 利用环绕加法，将带符号数 偏移至 无符号数的空间。
	tmp = _mm_adds_pu16(tmp, _mm_set1_pi16( (short)(0xFFFF-(iHigh+0x8000)) ));	// 限制最高值
	tmp = _mm_subs_pu16(tmp, _mm_set1_pi16( (short)(0xFFFF-(iHigh+0x8000)+(iLow+0x8000)) ));	// 限制最低值
	return _mm_add_pi16(tmp, _mm_set1_pi16( iLow ));	// 恢复偏移
}

// 将无符号紧缩字限制在[uLow,uHigh]区间
// 章节：8.7 数值的截取/8.7.2 对任意有符号数范围截取符号数
//
// result: 限制后的带符号紧缩字。伪代码——result[i]=(mm0[i]<uLow)?uLow:( (mm0[i]>uHigh)?uHigh:mm0[i] )。
// mm0: 源操作数。
inline __m64 md_clamp_muw(const __m64 mm0, unsigned short uLow, unsigned short uHigh)
{
	__m64 tmp = _mm_adds_pu16(mm0, _mm_set1_pi16( (short)(0xFFFFU-uHigh) ));	// 限制最高值
	tmp       = _mm_subs_pu16(tmp, _mm_set1_pi16( (short)(0xFFFFU-uHigh+uLow) ));	// 限制最低值
	return _mm_add_pi16(tmp, _mm_set1_pi16( uLow ));	// 恢复偏移
}

// 返回常数：0
// 章节：8.8 生成常量/[0]在MM0产生0寄存器
inline __m64 md_setzero_mmq()
{
	__m64 tmp=_mm_setzero_si64();	// 其实并不需要赋值，但为了符合语法，只能这样写。
	return _mm_xor_si64(tmp, tmp);
	// 其实Intrinsics函数中有这样的函数——
	// return _mm_setzero_si64();
}

// 返回常数：全1
// 章节：8.8 生成常量/[1]在寄存器MM1中置全1，它在每一个紧缩数据类型的值域中都是-1
inline __m64 md_setfull_mmq()
{
	__m64 tmp=_mm_setzero_si64();	// 其实并不需要赋值，但为了符合语法，只能这样写。
	return _mm_cmpeq_pi8(tmp, tmp);
}

// 返回常数：每个紧缩字节为1
// 章节：8.8 生成常量/[2]在每一个紧缩字节[或紧缩字](或紧缩双字)的值域中产生常数1
inline __m64 md_set_1_mib()
{
	__m64 mibZero = _mm_setzero_si64();
	__m64 mibNegativeOne = _mm_cmpeq_pi8(mibZero, mibZero);
	return _mm_sub_pi8(mibZero, mibNegativeOne);
}

// 返回常数：每个紧缩字为pow(2,n)-1
// 章节：8.8 生成常量/[3]在每一个紧缩字(或紧缩双字)的值域中产生带符号常数pow(2,n)-1
inline __m64 md_set_pow2n_sub1_miw(int n)
{
	assert((n>=1) && (n<=16));
	__m64 mibZero = _mm_setzero_si64();
	__m64 mibFull = _mm_cmpeq_pi8(mibZero, mibZero);
	return _mm_srli_pi16(mibFull, 16-n);
}

// 返回常数：每个紧缩字为-pow(2,n)
// 章节：8.8 生成常量/[4]在每一个紧缩字(或紧缩双字)的值域中产生带符号常数-pow(2,n)
inline __m64 md_set_neg_pow2n_miw(int n)
{
	assert((n>=0) && (n<=15));
	__m64 mibZero = _mm_setzero_si64();
	__m64 mibFull = _mm_cmpeq_pi8(mibZero, mibZero);
	return _mm_slli_pi16(mibFull, n);
}

// 验证
void doTest(int cnt)
{
	__m64 t0,t1,t2;
	int i;

	// 紧缩无符号字 解包为 两组紧缩无符号双字
	printf("md_unpack_mud4muw:\n");
	t0 = _mm_set_pi32(0x01234567, 0x89ABCDEF);
	printf("[%.8X%.8X] -> ", t0.m64_u32[1], t0.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t2 = md_unpack_mud4muw(t1, t0);
	}
	printf("[%.8X%.8X],[%.8X%.8X]\n", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	printf("\n");

	// 紧缩带符号字 解包为 两组紧缩带符号双字
	printf("md_unpack_mid4miw:\n");
	t0 = _mm_set_pi32(0x01234567, 0x89ABCDEF);
	printf("[%.8X%.8X] -> ", t0.m64_u32[1], t0.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t2 = md_unpack_mid4miw(t1, t0);
	}
	printf("[%.8X%.8X],[%.8X%.8X]\n", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	printf("\n");

	// 两组紧缩带符号双字 交叉饱和紧缩为 紧缩带符号字
	printf("md_pack_s_cross_miw4mid:\n");
	t1 = _mm_set_pi32(0x00001111, 0x000F2222);
	t2 = _mm_set_pi32(0xFFFFCCCC, 0xFFFFDDDD);
	printf("[%.8X%.8X],[%.8X%.8X] -> ", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t0 = md_pack_s_cross_miw4mid(t2, t1);
	}
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);
	printf("\n");

	// 两组紧缩无符号双字 交叉环绕紧缩为 紧缩无符号字
	printf("md_pack_w_cross_muw4mud:\n");
	t1 = _mm_set_pi32(0x00001111, 0x000F2222);
	t2 = _mm_set_pi32(0xFFFFCCCC, 0xFFFFDDDD);
	printf("[%.8X%.8X],[%.8X%.8X] -> ", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t0 = md_pack_w_cross_muw4mud(t2, t1);
	}
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);
	printf("\n");

	// 2x2矩阵转置.紧缩双字
	printf("md_matrix_transpose_2x2_mmd:\n");
	t1 = _mm_set_pi32(0x00001111, 0x000F2222);
	t2 = _mm_set_pi32(0xFFFFCCCC, 0xFFFFDDDD);
	printf("[%.8X%.8X],[%.8X%.8X] -> ", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		md_matrix_transpose_2x2_mmd(t1, t2);
	}
	printf("[%.8X%.8X],[%.8X%.8X]\n", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	printf("\n");

	// 复数与常量相乘（紧缩字->紧缩双字）
	printf("md_complex_mul_c_mid4miw:\n");
	t1 = _mm_set_pi16(0,0, 1, 1);	// 1+i
	t2 = _mm_set_pi16(3,-2, 2,3);	// 3+2i.	(1+i)*(3+2i) = 1+5i
	printf("[%.8X%.8X],[%.8X%.8X] -> ", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t0 = md_complex_mul_c_mid4miw(t1, t2);
	}
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);
	printf("\n");

	// 无符号紧缩字节的绝对差
	printf("md_absolute_deviation_mub:\n");
	t1 = _mm_set_pi8(1,2,3,4,5,6,7,8);
	t2 = _mm_set_pi8(8,7,6,5,4,3,2,1);
	printf("[%.8X%.8X],[%.8X%.8X] -> ", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t0 = md_absolute_deviation_mub(t1, t2);
	}
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);
	printf("\n");

	// 带符号紧缩字的绝对差
	printf("md_absolute_deviation_miw:\n");
	t1 = _mm_set_pi16(-1, 1, 3, 5);
	t2 = _mm_set_pi16( 2, 2, 2, 2);
	printf("[%.8X%.8X],[%.8X%.8X] -> ", t1.m64_u32[1], t1.m64_u32[0], t2.m64_u32[1], t2.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t0 = md_absolute_deviation_miw(t1, t2);
	}
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);
	printf("\n");

	// 带符号紧缩字的绝对值
	printf("md_abs_miw4miw:\n");
	t0 = _mm_set_pi16(-1, 1, 3, -5);
	printf("[%.8X%.8X] -> ", t0.m64_u32[1], t0.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t1 = md_abs_miw(t0);
	}
	printf("[%.8X%.8X]\n", t1.m64_u32[1], t1.m64_u32[0]);
	printf("\n");

	// 将带符号紧缩字限制在[iLow,iHigh]区间
	printf("md_clamp_miw:\n");
	t0 = _mm_set_pi16(-15, 1, 254, 257);
	printf("[%.8X%.8X] -> ", t0.m64_u32[1], t0.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t1 = md_clamp_miw(t0, -1, 255);
	}
	printf("[%.8X%.8X]\n", t1.m64_u32[1], t1.m64_u32[0]);
	printf("\n");

	// 将无符号紧缩字限制在[uLow,uHigh]区间
	printf("md_clamp_muw:\n");
	t0 = _mm_set_pi16(1, 254, 257, 32769U);
	printf("[%.8X%.8X] -> ", t0.m64_u32[1], t0.m64_u32[0]);
	for(i=0; i<cnt; ++i)
	{
		t1 = md_clamp_muw(t0, 16, 255);
	}
	printf("[%.8X%.8X]\n", t1.m64_u32[1], t1.m64_u32[0]);
	printf("\n");

	// 返回常数：0
	printf("md_setzero_mmq:\t");
	t0 = md_setzero_mmq();
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);

	// 返回常数：全1
	printf("md_setfull_mmq:\t");
	t0 = md_setfull_mmq();
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);

	// 返回常数：每个紧缩字节为1
	printf("md_set_1_mib:\t");
	t0 = md_set_1_mib();
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);

	// 返回常数：每个紧缩字为pow(2,n)-1
	printf("md_set_pow2n_sub1_miw:\t");
	t0 = md_set_pow2n_sub1_miw(8);
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);

	// 返回常数：每个紧缩字为pow(2,n)-1
	printf("md_set_neg_pow2n_miw:\t");
	t0 = md_set_neg_pow2n_miw(15);
	printf("[%.8X%.8X]\n", t0.m64_u32[1], t0.m64_u32[0]);

}

int main(int argc, char* argv[])
{
	doTest((rand()&1) + 1);	// 用一个随机数作为循环次数，避免编译器优化循环
	return 0;
}