x86上128位二进制乘法最快速算法征解 - 第5页 - 编程擂台 - 数学研发论坛

无心人 发表于 2008-3-18 14:00:51

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

无心人 发表于 2008-3-18 14:12:09

运行时间高于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
movd xmm2, dword ptr
pmuludq xmm0, xmm2
movd , xmm0
psrlq xmm0, 32
// 0:1 1:0
movd xmm1, dword ptr
movd xmm3, dword ptr
pmuludq xmm1, xmm3
movd xmm2, dword ptr
movd xmm4, dword ptr
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 , xmm1 //xmm1低位存储
psrlq xmm1, 32 //xmm1高位
paddq xmm0, xmm1 //得到进位
//1:1 0:2 2:0
movd xmm1, dword ptr
movd xmm4, dword ptr
pmuludq xmm1, xmm4
movd xmm2, dword ptr
movd xmm5, dword ptr
pmuludq xmm2, xmm5
movd xmm3, dword ptr
movd xmm7, dword ptr
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 , xmm6 //xmm6低位存储
psrlq xmm6, 32
paddq xmm0, xmm6 //得到进位
//0:3 1:2 2:1 3:0
movd xmm1, dword ptr
movd xmm5, dword ptr
pmuludq xmm1, xmm5
movd xmm2, dword ptr
movd xmm6, dword ptr
pmuludq xmm2, xmm6
movd xmm3, dword ptr
   movd xmm7, dword ptr
pmuludq xmm3, xmm7
movd xmm4, dword ptr
movd xmm6, dword ptr
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 , xmm0 //xmm0低位存储
psrlq xmm0, 32
paddq xmm0, xmm1 //新进位
//1:3 2:2 3:1
movd xmm2, dword ptr
movd xmm4, dword ptr
pmuludq xmm2, xmm4
movd xmm1, dword ptr
movd xmm5, dword ptr
pmuludq xmm1, xmm5
movd xmm3, dword ptr
movd xmm6, dword ptr
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 , xmm7 //xmm7低位存储
psrlq xmm7, 32
paddq xmm0, xmm7 //得到进位
//2:3 3:2
movd xmm1, dword ptr
movd xmm3, dword ptr
pmuludq xmm1, xmm3
movd xmm2, dword ptr
movd xmm4, dword ptr
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 , xmm7 //xmm7低位存储
psrlq xmm7, 32
paddq xmm0, xmm7 //得到进位
//3:3
movd xmm1, dword ptr
movd xmm2, dword ptr
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 , xmm6 //xmm6低位存储
psrlq xmm6, 32
paddq xmm0, xmm6 //进位
movd dword ptr , xmm0 //存储
// ret
}
}

无心人 发表于 2008-3-18 14:14:25

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

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

gxqcn 发表于 2008-3-18 14:26:44

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

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

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

在我的机器上测试结果如下：Test function: UInt128x128To256_ANSI_C32(..) 10000000 times...
Elapsed time: 696.680ms
FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF * FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF
= FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE 00000000 00000000 00000000 00000001

Test function: UInt128x128To256_SSE2_40F(..) 10000000 times...
Elapsed time: 483.896ms
FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF * FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF
= FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE 00000000 00000000 00000000 00000001

Test function: UInt128x128To256_SSE2_42F(..) 10000000 times...
Elapsed time: 826.723ms
FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF * FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF
= FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFE 00000000 00000000 00000000 00000001

请按任意键继续. . .测试环境：Windows XP SP2，AMD Athlon 64 Processor 3200+，1GB DDR - 200MHz

无心人 发表于 2008-3-18 14:32:49

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

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

无心人 发表于 2008-3-18 14:36:04

我没利用
_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:39:31

原帖由无心人于 2008-3-18 14:32 发表 http://images.5d6d.net/dz60/common/back.gif
嘿嘿
想把你CPU给顺来
比我这里快4倍啊

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

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

gxqcn 发表于 2008-3-18 14:41:25

回复 46# 的帖子

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

//...

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

无心人 发表于 2008-3-18 15:54:24

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

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

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

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

无心人 发表于 2008-3-18 17:52:54

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

mov esi, dword ptr
mov edi, dword ptr
mov ebx, dword ptr
movq xmm0, qword ptr
pmuludq xmm0, qword ptr
pcmpeqq xmm7, xmm7
psrlq xmm7, 32 //00000000FFFFFFFF00000000FFFFFFFF
movd , xmm0 //0:0结果低位保存
pshufd xmm1, qword ptr , 00010000b
pshufd xmm2, qword ptr , 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 , xmm2
psrldq xmm2, 4
pshufd xmm3, xmm2, 11111001b
paddq xmm2, xmm3
//2:0 0:2
pshufd xmm4, qword ptr , 00100000b
pshufd xmm5, qword ptr , 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 , 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 , 00010010b
pshufd xmm5, qword ptr , 00100001b
pmuludq xmm1, xmm5
pshufd xmm2, qword ptr , 00110000b
pshufd xmm6, qword ptr , 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 , xmm0
psrldq xmm0, 4
pshufd xmm3, xmm0, 11111001b
pshufd xmm0, xmm0, 11111100b
paddq xmm0, xmm3 //xmm0进位
//1:3 2:2 3:1
pshufd xmm6, , 00110010b
pshufd xmm5, , 00110010b
pmuludq xmm6, xmm5 //xmm6=3:3/2:2
pshufd xmm2, , 00110001b
pshufd xmm4, , 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 , xmm0
psrldq xmm0, 4
pshufd xmm4, xmm0, 11111001b
pshufd xmm0, xmm0, 11111100b
paddq xmm0, xmm4 //xmm0=进位
//2:3 3:2
pshufd xmm1, qword ptr , 00100011b
pshufd xmm2, qword ptr , 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 , xmm0
psrldq xmm0, 4
pshufd xmm2, xmm0, 11111001b
pshufd xmm0, xmm0, 11111100b
paddq xmm0, xmm2 //进位
pshufd xmm6, xmm6, 11011100b
paddq xmm0, xmm6
movd , xmm0
psrldq xmm0, 4
pshufd xmm1, xmm0, 11111001b
pshufd xmm0, xmm0, 11111100b
paddq xmm0, xmm1
movd , xmm0

页: 1 2 3 4 [5] 6 7 8 9 10 11 12 13 14

数学研发论坛's Archiver

回复 46# 的帖子