x86上128位二进制乘法最快速算法征解

无心人

使用C/C++代码和汇编
在指令集限于SSE2以下的情况下
(x86, MMX, SSE, SSE2代码)
对128X128 二进制整数运算写出最快速的代码

  函数形式

void  mul128(DWORD * result, DWORD * left, DWORD * right)

运算1000万次的累加时间在各机器上最少的为胜利

gxqcn

一年前，我曾发过千分求汇编优化：UInt96x96To192(...)， 吸引了不少高手参与，用到的汇编指令包括ALU、SSE2等（我本人的SSE2汇编就是从这个讨论入门的） 上周与楼主聊天时提到了这个。楼主想来个更大的：

void UInt128x128To256( UINT32 result[8], const UINT32 left[4], const UINT32 right[4] );

（我想像的函数形式，数组下标从小到大依次对应于变量的由低到高位） 楼主所给的函数原型不知是否与上述相同？ 对出入参数组的大小是否限定？即result是否限定于8个DWORD。。。

无心人

是 但最好写做指针格式 也不考虑内存分配啊等 只考虑4 * 4 =8算法

无心人

为了统一计时 函数计算 FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF * FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF = FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFE00000000000000000000000000000001 作为测试样例， 但函数应该能处理任何输入

无心人

看了你们的96位算法的讨论 发现两点问题 1、这么重要的问题, 代码里竟然没验证功能部分 2、没验证或者考虑MMX寄存器作为存储 PS: 可利用寄存器总共是下面几个 EAX, ECX, EDX, ECX ESI, EDI MM0-MM7 XMM0-XMM7

gxqcn

感觉这次擂台比上次的 UInt96x96To192() 更有意义，
为了大家方便交流和测试，我做了一个模板：仅一个.c文件再加一个dsp文件即可。

其中为了方便大家自测算法的正确性，我已完成了一个标准c的版本，
该版本的结果与我用HugeCalc的一致，并经过反复推敲论证过的，应该是无bug且高效的。
代码如下：

/************************************************************************/
/* UInt128x128To256.c                                                   */
/************************************************************************/

#include < wtypes.h >
#include < stdio.h >

#if 0
#   define _RAND_TEST   /* 随机测试 */
#endif

#ifndef _DEBUG          /* 测试次数 */
#  define _TEST_TIMES   10000000UL
#else
#  define _TEST_TIMES   1UL     /* 先DEBUG算法的正确性 */
#endif

typedef void ( *lpfn_UInt128x128To256 )( UINT32 * const result,
                                        const UINT32 * const left,
                                        const UINT32 * const right );

static BOOL s_bSupport_MMX = FALSE;
static BOOL s_bSupport_SSE = FALSE;
static BOOL s_bSupport_SSE2 = FALSE;

static UINT64 s_u64Frequency = 1;
static UINT64 s_u64Start, s_u64End;

static BYTE buffer[(4+4+8)*4+15];
static UINT32 * left = NULL;
static UINT32 * right = NULL;
static UINT32 * result = NULL;

void initParam( void )
{
    /* 16字节对齐，以便于 SSE/SSE2 优化 */
    left = (UINT32 *)( ((UINT32)( buffer+15 )) & -16 );
    right = left + 4;
    result = right + 4;
    /* 待测函数中不得假定三个指针的相对偏移量情况 */

#ifndef _RAND_TEST
    left[0] = 0xFFFFFFFF;
    left[1] = 0xFFFFFFFF;
    left[2] = 0xFFFFFFFF;
    left[3] = 0xFFFFFFFF;

    right[0] = 0xFFFFFFFF;
    right[1] = 0xFFFFFFFF;
    right[2] = 0xFFFFFFFF;
    right[3] = 0xFFFFFFFF;
#else
    /* 随机数确保可充满32bits */
#   define RAND_VAL()    (( (UINT32)(rand()) << 17 ) | rand())
    srand(GetTickCount());

    left[0] = RAND_VAL();
    left[1] = RAND_VAL();
    left[2] = RAND_VAL();
    left[3] = RAND_VAL();

    right[0] = RAND_VAL();
    right[1] = RAND_VAL();
    right[2] = RAND_VAL();
    right[3] = RAND_VAL();
#endif

    /* 测试 CPU 支持的指令集 */
    __asm
    {
        mov     eax, 1;
        cpuid;

        mov     ecx, 800000h;   /* 23 bit */
        and     ecx, edx;
        neg     ecx;
        sbb     ecx, ecx;
        neg     ecx;
        mov     dword ptr[ s_bSupport_MMX ], ecx;

        mov     ecx, 2000000h;  /* 25 bit */
        and     ecx, edx;
        neg     ecx;
        sbb     ecx, ecx;
        neg     ecx;
        mov     dword ptr[ s_bSupport_SSE ], ecx;

        mov     ecx, 4000000h;  /* 26 bit */
        and     ecx, edx;
        neg     ecx;
        sbb     ecx, ecx;
        neg     ecx;
        mov     dword ptr[ s_bSupport_SSE2 ], ecx;
    }

    QueryPerformanceFrequency((LARGE_INTEGER *)&s_u64Frequency );
}

/************************************************************************/
/* UInt128x128To256 ANSI C 版，经过了严格测试                           */
/************************************************************************/
void UInt128x128To256_ANSI_C32( UINT32 * const result,
                               const UINT32 * const left,
                               const UINT32 * const right )
{
    typedef union tag_UINT64
    {
        UINT64 u64Val;
        UINT32 u32LH[2];
    } UInt64;

    UInt64 u64_0x0, u64_0x1, u64_0x2, u64_0x3;
    UInt64 u64_1x0, u64_1x1, u64_1x2, u64_1x3;
    UInt64 u64_2x0, u64_2x1, u64_2x2, u64_2x3;
    UInt64 u64_3x0, u64_3x1, u64_3x2, u64_3x3;

    u64_0x0.u64Val = UInt32x32To64( left[0], right[0] );
    u64_0x1.u64Val = UInt32x32To64( left[0], right[1] );
    u64_0x2.u64Val = UInt32x32To64( left[0], right[2] );
    u64_0x3.u64Val = UInt32x32To64( left[0], right[3] );

    u64_1x0.u64Val = UInt32x32To64( left[1], right[0] );
    u64_1x1.u64Val = UInt32x32To64( left[1], right[1] );
    u64_1x2.u64Val = UInt32x32To64( left[1], right[2] );
    u64_1x3.u64Val = UInt32x32To64( left[1], right[3] );

    u64_2x0.u64Val = UInt32x32To64( left[2], right[0] );
    u64_2x1.u64Val = UInt32x32To64( left[2], right[1] );
    u64_2x2.u64Val = UInt32x32To64( left[2], right[2] );
    u64_2x3.u64Val = UInt32x32To64( left[2], right[3] );

    u64_3x0.u64Val = UInt32x32To64( left[3], right[0] );
    u64_3x1.u64Val = UInt32x32To64( left[3], right[1] );
    u64_3x2.u64Val = UInt32x32To64( left[3], right[2] );
    u64_3x3.u64Val = UInt32x32To64( left[3], right[3] );

    /*                                        FF FE  00 01 --[0][0]
                                       FF FE  00 01 --[0][1]
                                       FF FE  00 01 --[1][0]
                                FF FE  00 01 --[0][2]
                                FF FE  00 01 --[1][1]
                                FF FE  00 01 --[2][0]
                         FF FE  00 01 --[0][3]
                         FF FE  00 01 --[1][2]
                         FF FE  00 01 --[2][1]
                         FF FE  00 01 --[3][0]
                  FF FE  00 01 --[1][3]
                  FF FE  00 01 --[2][2]
                  FF FE  00 01 --[3][1]
           FF FE  00 01 --[2][3]
           FF FE  00 01 --[3][2]
    FF FE  00 01 --[3][3]  表示  FFFF FFFE  0000 0001*/
    u64_0x1.u64Val += u64_0x0.u32LH[1];

    u64_0x1.u64Val += u64_1x0.u64Val;
    u64_2x0.u32LH[1] += ( u64_0x1.u64Val < u64_1x0.u64Val );
    u64_0x2.u64Val += u64_0x1.u32LH[1];

    u64_0x2.u64Val += u64_1x1.u64Val;
    u64_3x0.u32LH[1] += ( u64_0x2.u64Val < u64_1x1.u64Val );
    u64_0x2.u64Val += u64_2x0.u64Val;
    u64_2x1.u32LH[1] += ( u64_0x2.u64Val < u64_2x0.u64Val );
    u64_0x3.u64Val += u64_0x2.u32LH[1];

    u64_0x3.u64Val += u64_1x2.u64Val;
    u64_3x1.u32LH[1] += ( u64_0x3.u64Val < u64_1x2.u64Val );
    u64_0x3.u64Val += u64_2x1.u64Val;
    u64_2x2.u32LH[1] += ( u64_0x3.u64Val < u64_2x1.u64Val );
    u64_0x3.u64Val += u64_3x0.u64Val;
    u64_2x3.u64Val += ( u64_0x3.u64Val < u64_3x0.u64Val );
    u64_1x3.u64Val += u64_0x3.u32LH[1];

    u64_1x3.u64Val += u64_2x2.u64Val;
    u64_3x2.u32LH[1] += ( u64_1x3.u64Val < u64_2x2.u64Val );
    u64_1x3.u64Val += u64_3x1.u64Val;
    u64_3x3.u64Val += ( u64_1x3.u64Val < u64_3x1.u64Val );
    u64_2x3.u64Val += u64_1x3.u32LH[1];

    u64_2x3.u64Val += u64_3x2.u64Val;
    u64_3x3.u32LH[1] += ( u64_2x3.u64Val < u64_3x2.u64Val );
    u64_3x3.u64Val += u64_2x3.u32LH[1];

    result[0] = u64_0x0.u32LH[0];
    result[1] = u64_0x1.u32LH[0];
    result[2] = u64_0x2.u32LH[0];
    result[3] = u64_0x3.u32LH[0];
    result[4] = u64_1x3.u32LH[0];
    result[5] = u64_2x3.u32LH[0];
    result[6] = u64_3x3.u32LH[0];
    result[7] = u64_3x3.u32LH[1];
}
/************************************************************************/
/* 测试函数代码粘贴区开始 begin{                                        */
/*----------------------------------------------------------------------*/
/* 待测函数命名规范(推荐)：UInt128x128To256_{1}_{2}(...)                */
/*     其中{1}代表用到的最高级指令集；                                  */
/*         {2}代表发帖楼层，如果写错了可以自行修订或请管理员帮助修改    */
/* 以后大家只需发自己待测函数的代码即可，不必再贴全测试代码             */
/************************************************************************/

_declspec(naked)
void UInt128x128To256_MMX_xxxF( UINT32 * const result,
                               const UINT32 * const left,
                               const UINT32 * const right )
{
    __asm
    {
        /* do something ... */


        ret;
    }
}

_declspec(naked)
void UInt128x128To256_SSE_xxxF( UINT32 * const result,
                               const UINT32 * const left,
                               const UINT32 * const right )
{
    __asm
    {
        /* do something ... */


        ret;
    }
}

_declspec(naked)
void UInt128x128To256_SSE2_11F( UINT32 * const result,
                               const UINT32 * const left,
                               const UINT32 * const right )
{
    __asm
    {
        /* 具体代码在 11#，需登陆才可见 */
        ret;
    }
}

/************************************************************************/
/* }end 测试函数代码粘贴区结束                                          */
/************************************************************************/

void testFun( const lpfn_UInt128x128To256 pFun,
             const LPCTSTR lpszFunName,
             const UINT32 u32TestTimes )
{
    UINT32 i;

    printf( "\nTest function: %s(..) %u times... \n",
        lpszFunName, u32TestTimes );

    QueryPerformanceCounter((LARGE_INTEGER *)&s_u64Start );

    for ( i = 0; i < u32TestTimes; ++i )
    {
        (*pFun)( result, left, right );
    }

    QueryPerformanceCounter((LARGE_INTEGER *)&s_u64End );

    i = (UINT32)(( s_u64End - s_u64Start ) * 1000000UL / s_u64Frequency );
    printf( "Elapsed time: %d.%03u ms\n", i / 1000, i % 1000 );

    printf( "  %08X %08X %08X %08X * %08X %08X %08X %08X\n"
        "= %08X %08X %08X %08X   %08X %08X %08X %08X\n",
        left[3], left[2], left[1], left[0],
        right[3], right[2], right[1], right[0],
        result[7], result[6], result[5], result[4],
        result[3], result[2], result[1], result[0] );
}

int main(int argc, char* argv[])
{
    initParam();

    /* 标准结果 */
    testFun( UInt128x128To256_ANSI_C32,
            "UInt128x128To256_ANSI_C32", /*1*/ _TEST_TIMES );

    /* MMX 版本测试 */
    if ( s_bSupport_MMX )
    {
        testFun( UInt128x128To256_MMX_xxxF,
                "UInt128x128To256_MMX_xxxF", _TEST_TIMES );

        /* test other functions:
        testFun( ... );
        */
    }

    /* SSE 版本测试 */
    if ( s_bSupport_SSE )
    {
        testFun( UInt128x128To256_SSE_xxxF,
                "UInt128x128To256_SSE_xxxF", _TEST_TIMES );

        /* test other functions:
        testFun( ... );
        */
    }

    /* SSE2 版本测试 */
    if ( s_bSupport_SSE2 )
    {
        testFun( UInt128x128To256_SSE2_11F,
                "UInt128x128To256_SSE2_11F", _TEST_TIMES );

        /* test other functions:
        testFun( ... );
        */
    }

    printf( "\n" );
    system( "pause" );
    return 0;
}

完整的测试模板包如下： UInt128x128To256.zip (3.75 KB, 下载次数: 12)

2008-3-18 19:19 上传

点击文件名下载附件
UInt128x128To256 测试模板

（包含一个.c文件，及dsp文件）

请有兴趣的朋友在本地机上编译运行；如发现任何问题，请及时反馈，以便修正，谢谢！

无心人

一共四个字 l0, l0, l2, l3 r0, r1, r2, r3 esi = l edi = r ebx = result

无心人

下面是乘法三角形 l3 r3 l2 r2 l1 r1 l0 r0 l2 r3 l1 r2 l0 r1 l3 r2 l2 r1 l1 r0 l1 r3 l0 r2 l3 r1 l2 r0 l0 r3 l3 r0 左面是最高位 每个数字宽度为32位 从左面开始 相邻符号表示两个数字乘积 比如l3 r3表示 l3 * r3

无心人

可以判断 从最右上角乘并向左推进是最节约的方法

无心人

现在的问题是 存在最多4个四字加 这么做? 假设mm0 + mm1 + mm2 + mm3 xor mm4, mm4 xor mm5, mm5 xor mm6. mm6 xor mm7, mm7 movd mm4, mm0 pshrq mm0, 32 movd mm5, mm1 pshrq mm1, 32 addq mm4, mm5 movd mm6, mm2 psrq mm2, 32 addq mm4, mm6 movd mm7, mm3 psrq mm3, 32 addq mm4, mm7 movd [result + x], mm4 pshrq mm4, 32 addq mm0, mm1 addq mm2, mm3 addq mm0, mm2 addq mm0, mm4 现在结果在mm0:[result+x]