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
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
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
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