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

http://blog.csdn.net/zyl910/article/details/7516127

一、来源

　　来源：《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）。
inlinevoid 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;
}