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

无心人

共143条指令
全SSE2程序指令数必须小于该数字才有效率

无心人

运行时间高于MMX寄存器版本
指令数多16条
void UInt128x128To256_SSE2_42F( UINT32 * const result,
                                                           const UINT32 * const left,
                                                           const UINT32 * const right )
{
//使用SSE寄存器， 用SSE2指令， 但只利用低64位
        __asm
        {
                mov esi, left
                mov edi, right
                mov ebx, result
//0:0
        movd xmm0, dword ptr [esi]
                movd xmm2, dword ptr [edi]
                pmuludq xmm0, xmm2
                movd [ebx], xmm0
                psrlq xmm0, 32
// 0:1 1:0
                movd xmm1, dword ptr [esi]
                movd xmm3, dword ptr [edi+4]
                pmuludq xmm1, xmm3
                movd xmm2, dword ptr [edi]
                movd xmm4, dword ptr [esi+4]
                pmuludq xmm2, xmm4
                mov eax, -1
//xmm0+xmm1+xmm2
                paddq xmm0, xmm1//xmm0=xmm0+xmm1 xmm0肯定小于 2^64
                movd xmm1, eax
                pand xmm1, xmm0
                psrlq xmm0, 32 //分解xmm0 = xmm0:xmm1
                movd xmm3, eax
                pand xmm3, xmm2
                psrlq xmm2, 32 //分解xmm2=xmm2:xmm3
                paddq xmm1, xmm3 //xmm1=xmm1+xmm3
                paddq xmm0, xmm2 //xmm0=xmm0+xmm2
                movd dword ptr [ebx+4], xmm1 //xmm1低位存储
                psrlq xmm1, 32 //xmm1高位
                paddq xmm0, xmm1 //得到进位
//1:1 0:2 2:0
                movd xmm1, dword ptr [esi+4]
                movd xmm4, dword ptr [edi+4]
                pmuludq xmm1, xmm4
                movd xmm2, dword ptr [esi]
                movd xmm5, dword ptr [edi+8]
                pmuludq xmm2, xmm5
                movd xmm3, dword ptr [edi]
                movd xmm7, dword ptr [esi+8]
                pmuludq xmm3, xmm7
//xmm0+xmm1+xmm2+xmm3
                movd xmm6, eax
                pand xmm6, xmm0
                psrlq xmm0, 32 //分解xmm0=xmm0:xmm6
                movd xmm5, eax
                pand xmm5, xmm1
                psrlq xmm1, 32 //分解xmm1=xmm1:xmm5
                paddq xmm6, xmm5 //xmm6=xmm6+xmm5
                paddq xmm0, xmm1 //xmm0=xmm0+xmm1
                movd xmm4, eax
                pand xmm4, xmm2
                psrlq xmm2, 32 //分解xmm2=xmm2:xmm4
                paddq xmm6, xmm4 //xmm6=xmm6+xmm4
                paddq xmm0, xmm2 //xmm0=xmm0+xmm2
                movd xmm7, eax
                pand xmm7, xmm3
                psrlq xmm3, 32 //分解xmm3=xmm3:xmm7
                paddq xmm6, xmm7 //xmm6=xmm6+xmm7
                paddq xmm0, xmm3 //xmm0=xmm0+xmm3
                movd dword ptr [ebx+8], xmm6 //xmm6低位存储
                psrlq xmm6, 32  
                paddq xmm0, xmm6 //得到进位
//0:3 1:2 2:1 3:0
                movd xmm1, dword ptr [esi]
                movd xmm5, dword ptr [edi+12]
                pmuludq xmm1, xmm5
                movd xmm2, dword ptr [esi+4]
                movd xmm6, dword ptr [edi+8]
                pmuludq xmm2, xmm6
                movd xmm3, dword ptr [esi+8]
        movd xmm7, dword ptr [edi+4]
                pmuludq xmm3, xmm7
                movd xmm4, dword ptr [esi+12]
                movd xmm6, dword ptr [edi]
                pmuludq xmm4, xmm6
//xmm0+xmm1+xmm2+xmm3+xmm4
                movd xmm5, eax
                pand xmm5, xmm1
                psrlq xmm1, 32 //分解xmm1=xmm1:xmm5
                paddq xmm0, xmm5 //xmm0=xmm0+xmm5
                movd xmm6, eax
                pand xmm6, xmm2
                psrlq xmm2, 32 //分解xmm2=xmm2:xmm6
                paddq xmm0, xmm6 //xmm0=xmm0+xmm6
                paddq xmm1, xmm2 //xmm1=xmm1+xmm2
                movd xmm7, eax
                pand xmm7, xmm3
                psrlq xmm3, 32 //分解xmm3=xmm3:xmm7
                paddq xmm0, xmm7 //xmm0=xmm0+xmm7
                paddq xmm1, xmm3 //xmm1=xmm1+xmm3
                movd xmm5, eax
                pand xmm5, xmm4
                psrlq xmm4, 32 //分解xmm4=xmm4:xmm5
                paddq xmm0, xmm5 //xmm0=xmm0+xmm5
                paddq xmm1, xmm4 //xmm1=xmm1+xmm4
                movd dword ptr [ebx+12], xmm0 //xmm0低位存储
                psrlq xmm0, 32
                paddq xmm0, xmm1 //新进位
//1:3 2:2 3:1
                movd xmm2, dword ptr [esi+4]
                movd xmm4, dword ptr [edi+12]
                pmuludq xmm2, xmm4
                movd xmm1, dword ptr [esi+8]
                movd xmm5, dword ptr [edi+8]
                pmuludq xmm1, xmm5
                movd xmm3, dword ptr [esi+12]
                movd xmm6, dword ptr [edi+4]
                pmuludq xmm3, xmm6
//xmm0+xmm1+xmm2+xmm3
                movd xmm7, eax
                pand xmm7, xmm0
                psrlq xmm0, 32 //分解xmm0=xmm0:xmm7
                movd xmm6, eax
                pand xmm6, xmm1
                psrlq xmm1, 32 //分解xmm1=xmm1:xmm6
                paddq xmm7, xmm6 //xmm7=xmm7+xmm6
                paddq xmm0, xmm1 //xmm0=xmm0+xmm1
                movd xmm5, eax
                pand xmm5, xmm2
                psrlq xmm2, 32 //分解xmm2=xmm2:xmm5
                paddq xmm7, xmm5 //xmm7=xmm7+xmm5
                paddq xmm0, xmm2 //xmm0=xmm0+xmm2
                movd xmm4, eax
                pand xmm4, xmm3
                psrlq xmm3, 32 //分解xmm3=xmm3:xmm4
                paddq xmm7, xmm4 //xmm7=xmm7+xmm4
                paddq xmm0, xmm3 //xmm0=xmm0+xmm3
                movd dword ptr [ebx+16], xmm7 //xmm7低位存储
                psrlq xmm7, 32
                paddq xmm0, xmm7 //得到进位
//2:3 3:2
                movd xmm1, dword ptr [esi+8]
                movd xmm3, dword ptr [edi+12]
                pmuludq xmm1, xmm3
                movd xmm2, dword ptr [esi+12]
                movd xmm4, dword ptr [edi+8]
                pmuludq xmm2, xmm4
//xmm0+xmm1+xmm2
                movd xmm7, eax
                pand xmm7, xmm0
                psrlq xmm0, 32 //分解xmm0=xmm0:xmm7
                movd xmm6, eax
                pand xmm6, xmm1
                psrlq xmm1, 32 //分解xmm1=xmm1:xmm6
                paddq xmm7, xmm6 //xmm7=xmm7+xmm6
                paddq xmm0, xmm1 //xmm0=xmm0+xmm1
                movd xmm5, eax
                pand xmm5, xmm2
                psrlq xmm2, 32 //分解xmm2=xmm2:xmm5
        paddq xmm7, xmm5 //xmm7=xmm7+xmm5
                paddq xmm0, xmm2 //xmm0=xmm0+xmm2
                movd dword ptr [ebx+20], xmm7 //xmm7低位存储
                psrlq xmm7, 32
                paddq xmm0, xmm7 //得到进位
//3:3
                movd xmm1, dword ptr [esi+12]
                movd xmm2, dword ptr [edi+12]
                pmuludq xmm1, xmm2
//xmm0+xmm1
                movd xmm6, eax
        pand xmm6, xmm0
                psrlq xmm0, 32 //分解xmm0=xmm0:xmm6
                movd xmm7, eax
                pand xmm7, xmm1
                psrlq xmm1, 32 //分解xmm1=xmm1:xmm7
                paddq xmm6, xmm7 //xmm6=xmm6+xmm7
                paddq xmm0, xmm1 //xmm0=xmm0+xmm1
                movd dword ptr [ebx+24], xmm6 //xmm6低位存储
                psrlq xmm6, 32
                paddq xmm0, xmm6 //进位
                movd dword ptr [ebx+28], xmm0 //存储
//                ret
        }
}

无心人

1000万次
C版本          5132338
MMX版SSE2    1585527
XMM版SSE2        1753418

全SSE2双乘版暂时没写
上面两个版本还可优化
==========================
另外 GxQ给的模板存在错误
时间单位是us而不是ms

gxqcn

谢谢指正错误！我晚上回家后去修正那个附件。

由于想让测试程序编译好后可在各种机台上直接运行，
所以测试模版中特意增加了指令集支持与否的自动测试。

而算法是否可被正常运行，关键是所用指令是否被支持，所以我将楼主的两个测试函数名做了修正。

在我的机器上测试结果如下：

1. Test function: UInt128x128To256_ANSI_C32(..) 10000000 times...
   
2. Elapsed time: 696.680ms
   
3.   FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF * FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF
   
4. = FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE   00000000 00000000 00000000 00000001
   
5. 
   
6. Test function: UInt128x128To256_SSE2_40F(..) 10000000 times...
   
7. Elapsed time: 483.896ms
   
8.   FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF * FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF
   
9. = FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE   00000000 00000000 00000000 00000001
   
10. 
    
11. Test function: UInt128x128To256_SSE2_42F(..) 10000000 times...
    
12. Elapsed time: 826.723ms
    
13.   FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF * FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF
    
14. = FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE   00000000 00000000 00000000 00000001
    
15. 
    
16. 请按任意键继续. . .

复制代码

测试环境：Windows XP SP2，AMD Athlon 64 Processor 3200+，1GB DDR - 200MHz

无心人

嘿嘿
想把你CPU给顺来
比我这里快4倍啊

另外C版本你的快很多
编译器优化？？？

无心人

我没利用
_declspec(naked)

利用了是否要加
push ebp
mov ebp, esp
push ebx
push esi
push edi
......
,,,,,,
pop edi
pop esi
pop ebx
pop ebp
ret

gxqcn

原帖由 无心人 于 2008-3-18 14:32 发表
嘿嘿
想把你CPU给顺来
比我这里快4倍啊

另外C版本你的快很多
编译器优化？？？

我的编译环境是 VC6.0 + ICC10.0.027，编译选项全部为默认，未作任何优化。
在我机器上C版本快很多，估计又与 AMD 的 adc 指令较快相关。

gxqcn

应该是：
        _asm
        {
                push        ebx;
                push        ebp;
                push        esi;
                push        edi;

                //...

                pop                edi;
                pop                esi;
                pop                ebp;
                pop                ebx;
                ret;
        }
寄存器 esp 也要保持平衡。

无心人

哦， 差不了三瓜俩枣的， 就不加了

我刚又看了你们的96bit的讨论
现在问你
MMX寄存器的版本和那里的比较差多少？

刚仔细搜索了SSE指令
发现感兴趣的指令全是SSE4的
我倒啊

就PSHUFD PAND PADDQ PMULUDQ几个可利用
乘加全是16位的就没32位的乘加
否则可大大简化运算

无心人

SSE2代码终于在回家前赶完，明天调试
比预想的简单

mov esi, dword ptr [left]
  mov edi, dword ptr [right]
  mov ebx, dword ptr [result]
  movq xmm0, qword ptr [esi]
  pmuludq xmm0, qword ptr [edi]
  pcmpeqq xmm7, xmm7
   psrlq xmm7, 32 //00000000FFFFFFFF00000000FFFFFFFF
  movd [ebx], xmm0 //0:0结果低位保存
  pshufd xmm1, qword ptr [esi], 00010000b
  pshufd xmm2, qword ptr [edi], 00000010b
  pmuludq xmm1, xmm2 //xmm1= 1:0/0:1
  movq xmm4, xmm0
        pshufd xmm4, xmm4, 11111101b
  pshufd xmm3, xmm1, 00010000b
  pand xmm3, xmm7
  pshufd xmm2, xmm1, 00110010b
  pand xmm2, xmm7
  paddq xmm2, xmm3 //0:1+1:0=xmm2
  paddq xmm2, xmm4
  psrldq xmm0, 8
//xmm0=1:1
  movd [ebx+4], xmm2
  psrldq xmm2, 4
  pshufd xmm3, xmm2, 11111001b
  paddq xmm2, xmm3
//2:0 0:2
  pshufd xmm4, qword ptr [esi], 00100000b
  pshufd xmm5, qword ptr [edi], 00000010b
  pmuludq xmm4, xmm5 //xmm4=2:0/0:2
  movq xmm2, xmm2 //xmm2=进位
  pshufd xmm3, xmm4, 00010000b
  pand xmm3, xmm7
  pshufd xmm4, xmm4, 00110010b
  pand xmm4, xmm7
  paddq xmm3, xmm4 //xmm3=2:0 + 0:2 xmm0=1:1 xmm2=进位
  pshufd xmm0, xmm0, 11011000b
  paddq xmm3, xmm2
  paddq xmm3, xmm0
  movd [ebx+8], xmm3
  psrldq xmm3, 4
  pshufd xmm0, xmm3, 11111001b
  pshufd xmm3, xmm3, 11111100b
  paddq xmm0, xmm3 //xmm0=进位
//1:2 2:1 3:0 0:3
  pshufd xmm1, qword ptr [esi], 00010010b
  pshufd xmm5, qword ptr [edi], 00100001b
  pmuludq xmm1, xmm5
  pshufd xmm2, qword ptr [esi], 00110000b
  pshufd xmm6, qword ptr [edi], 00000011b
  pmuludq xmm2, xmm6 //xmm0进位 xmm1=1:2/2:1 xmm2=3:0/0:3
  pshufd xmm3, xmm1, 00010000b
  pand xmm3, xmm7
  pshufd xmm1, xmm1, 00110010b
  pand xmm1, xmm7
  paddq xmm1, xmm3
  pshufd xmm4, xmm2, 00010000b
  pand xmm4, xmm7
  pshufd xmm2, xmm2, 00110010b
  pand xmm2, xmm7
  paddq xmm2, xmm4
  paddq xmm0, xmm1
  paddq xmm0, xmm2
  movd [ebx+12], xmm0
  psrldq xmm0, 4
  pshufd xmm3, xmm0, 11111001b
  pshufd xmm0, xmm0, 11111100b
  paddq xmm0, xmm3 //xmm0进位
//1:3 2:2 3:1
  pshufd xmm6, [esi], 00110010b
  pshufd xmm5, [edi], 00110010b
  pmuludq xmm6, xmm5 //xmm6=3:3/2:2
  pshufd xmm2, [esi], 00110001b
  pshufd xmm4, [edi], 00010011b
  pmuludq xmm2, xmm4 //xmm2=3:1/1:3
  movq xmm1, xmm6 //xmm1=2:2
  psrldq xmm6, 8 //xmm6=3:3
  pshufd xmm3, xmm2, 00110001b
  pand xmm3, xmm7
  pshufd xmm2, xmm2, 00010010b
  pand xmm2, xmm7
  paddq xmm2, xmm3
  pshufd xmm1, xmm1, 11011100b
  paddq xmm1, xmm2
  paddq xmm0, xmm1
  movd [ebx+16], xmm0
  psrldq xmm0, 4
  pshufd xmm4, xmm0, 11111001b
  pshufd xmm0, xmm0, 11111100b
  paddq xmm0, xmm4 //xmm0=进位
//2:3 3:2
  pshufd xmm1, qword ptr [esi], 00100011b
  pshufd xmm2, qword ptr [edi], 00110010b
  pmuludq xmm1, xmm2
  pshufd xmm3, xmm1, 00010000b
  pand xmm3, xmm7
  pshufd xmm1, xmm1, 00110010b
  pand xmm1, xmm7
  paddq xmm1, xmm3
  paddq xmm0, xmm1
  movd [ebx+20], xmm0
  psrldq xmm0, 4
  pshufd xmm2, xmm0, 11111001b
  pshufd xmm0, xmm0, 11111100b
  paddq xmm0, xmm2 //进位
  pshufd xmm6, xmm6, 11011100b
  paddq xmm0, xmm6
  movd [ebx+24], xmm0
  psrldq xmm0, 4
  pshufd xmm1, xmm0, 11111001b
  pshufd xmm0, xmm0, 11111100b
  paddq xmm0, xmm1
  movd [ebx+28], xmm0