数学研发论坛

 找回密码
 欢迎注册
查看: 36694|回复: 138

[擂台] x86上128位二进制乘法最快速算法征解

[复制链接]
发表于 2008-3-3 20:59:48 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

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

  函数形式

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

运算1000万次的累加时间在各机器上最少的为胜利
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-3-3 21:37:33 | 显示全部楼层
一年前,我曾发过千分求汇编优化:UInt96x96To192(...)
吸引了不少高手参与,用到的汇编指令包括ALU、SSE2等(我本人的SSE2汇编就是从这个讨论入门的)

上周与楼主聊天时提到了这个。楼主想来个更大的:
void UInt128x128To256( UINT32 result[8], const UINT32 left[4], const UINT32 right[4] );
(我想像的函数形式,数组下标从小到大依次对应于变量的由低到高位)

楼主所给的函数原型不知是否与上述相同?
对出入参数组的大小是否限定?即result是否限定于8个DWORD。。。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-3 21:45:20 | 显示全部楼层

但最好写做指针格式
也不考虑内存分配啊等

只考虑4 * 4 =8算法
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-3 21:47:58 | 显示全部楼层
为了统一计时

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

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-4 08:09:30 | 显示全部楼层
看了你们的96位算法的讨论
发现两点问题
1、这么重要的问题, 代码里竟然没验证功能部分
2、没验证或者考虑MMX寄存器作为存储


PS:
可利用寄存器总共是下面几个
EAX, ECX, EDX, ECX
ESI, EDI
MM0-MM7
XMM0-XMM7
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-3-5 19:41:41 | 显示全部楼层

测试程序模板

感觉这次擂台比上次的 UInt96x96To192() 更有意义,
<br>为了大家方便交流和测试,我做了一个模板:仅一个.c文件再加一个dsp文件即可。
<br>
<br>其中为了方便大家自测算法的正确性,我已完成了一个标准c的版本,
<br>该版本的结果与我用HugeCalc的一致,并经过反复推敲论证过的,应该是无bug且高效的。
<br>
代码如下:
<pre>
/************************************************************************/
/* 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;
}</pre>
推荐


<br>完整的测试模板包如下: UInt128x128To256.zip (3.75 KB, 下载次数: 12)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-5 19:55:11 | 显示全部楼层
一共四个字
l0, l0, l2, l3
r0, r1, r2, r3

esi = l
edi = r
ebx = result
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-5 20:06:11 | 显示全部楼层
下面是乘法三角形

  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
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-5 20:11:59 | 显示全部楼层
可以判断
从最右上角乘并向左推进是最节约的方法
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-3-5 20:23:44 | 显示全部楼层
现在的问题是
存在最多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]
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-9-18 03:05 , Processed in 0.062073 second(s), 19 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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