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

liangbch · 发表于 2008-3-24 12:22:34

在PIV上运行， UInt128x128To256_SSE2_58F 依然是最快的版本，下面是在PIV 2.6G上运行附件程序的测试结果。

Test function: UInt128x128To256_ANSI_C32(..) 10000000 times...
Elapsed time: 1922.926 ms
  61E4210A 30360E4B 277E205E 76686B79 * B17E7E8A 5DAC7050 2D821227 A31462E1
= 43DF1983 7F67BDEE F4F9FCCA FF5AE6BD 6ACA3B72 D0384E14 E4ED3A5F 7B5EC759

Test function: UInt128x128To256_SSE2_40F(..) 10000000 times...
Elapsed time: 963.448 ms
  61E4210A 30360E4B 277E205E 76686B79 * B17E7E8A 5DAC7050 2D821227 A31462E1
= 43DF1983 7F67BDEE F4F9FCCA FF5AE6BD 6ACA3B72 D0384E14 E4ED3A5F 7B5EC759

Test function: UInt128x128To256_SSE2_42F(..) 10000000 times...
Elapsed time: 723.966 ms
  61E4210A 30360E4B 277E205E 76686B79 * B17E7E8A 5DAC7050 2D821227 A31462E1
= 43DF1983 7F67BDEE F4F9FCCA FF5AE6BD 6ACA3B72 D0384E14 E4ED3A5F 7B5EC759

Test function: UInt128x128To256_SSE2_54F(..) 10000000 times...
Elapsed time: 688.926 ms
  61E4210A 30360E4B 277E205E 76686B79 * B17E7E8A 5DAC7050 2D821227 A31462E1
= 43DF1983 7F67BDEE F4F9FCCA FF5AE6BD 6ACA3B72 D0384E14 E4ED3A5F 7B5EC759

Test function: UInt128x128To256_SSE2_56F(..) 10000000 times...
Elapsed time: 578.551 ms
  61E4210A 30360E4B 277E205E 76686B79 * B17E7E8A 5DAC7050 2D821227 A31462E1
= 43DF1983 7F67BDEE F4F9FCCA FF5AE6BD 6ACA3B72 D0384E14 E4ED3A5F 7B5EC759

Test function: UInt128x128To256_SSE2_58F(..) 10000000 times...
Elapsed time: 496.951 ms
  61E4210A 30360E4B 277E205E 76686B79 * B17E7E8A 5DAC7050 2D821227 A31462E1
= 43DF1983 7F67BDEE F4F9FCCA FF5AE6BD 6ACA3B72 D0384E14 E4ED3A5F 7B5EC759

gxqcn · 发表于 2008-3-24 12:47:36

谢谢楼上的测试！

看来不同机器对指令的处理周期相差很大，
也许在 Core2 上，56F 反超 58F 也是可能的。。。

无心人 · 发表于 2008-3-24 13:28:34

你的P4是31流水线的P4E
俺的是P4
郭的是Althon
差Core 2版
：）

无心人 · 发表于 2008-3-24 13:33:30

事情复杂了
我执行GxQcn程序也得到liangbch的结果
但我测试却不是
现在
问题在什么地方？

ICC编译器难道对汇编也动？

无心人 · 发表于 2008-3-24 13:46:12

GxQcn你能不能把存储临时值操作取消掉
那个高延时啊

gxqcn · 发表于 2008-3-24 13:57:49

Re: 64#
反正有完整的测试源代码，你不妨将不同编译器编译的程序运行一下，并把结果公布出来，让大家分析分析。。。

Re: 65#
我也不想啊，寄存器那么少，这是权衡之后的下下策啊。

无心人 · 发表于 2008-3-24 14:22:17

疯了疯了
在我VC2008Express上重新编译你代码
结果和你程序一致
也是依次减小

为什么啊

无心人 · 发表于 2008-3-24 16:51:38

这么做有优势么？
先乘四个3:3 2:2 1:1 0:0
正好存在结果里
然后依次乘剩下的12个
加在结果上

gxqcn · 发表于 2008-3-24 19:02:21

对 58F 再进行指令混排（寄存器名局部有轮换），得到新的版本如下：

_declspec(naked)
void UInt128x128To256_SSE2_69F( UINT32 * const result,
const UINT32 * const left,
const UINT32 * const right )
{
__asm
{
mov ecx, dword ptr[esp + 04h] ; result
mov eax, dword ptr[esp + 08h] ; left
mov edx, dword ptr[esp + 0Ch] ; right
movdqa xmm5, xmmword ptr[eax] ; load left
movdqa xmm4, xmmword ptr[edx] ; load right
mov eax, esp ;
sub esp, 0x4F ;
and esp, -16 ;
pshufd xmm0, xmm5, 00000000b ; xmm0 = L0:L0:L0:L0
pshufd xmm1, xmm5, 01010101b ; xmm1 = L1:L1:L1:L1
pshufd xmm2, xmm5, 10101010b ; xmm2 = L2:L2:L2:L2
pshufd xmm3, xmm5, 11111111b ; xmm3 = L3:L3:L3:L3
movdqa xmmword ptr[ecx+0x00], xmm0 ; xmm0 = L0:L0:L0:L0
movdqa xmmword ptr[ecx+0x10], xmm1 ; xmm1 = L1:L1:L1:L1
movdqa xmm6, xmm2 ; xmm6 = L2:L2:L2:L2
movdqa xmm7, xmm3 ; xmm7 = L3:L3:L3:L3
pmuludq xmm0, xmm4 ; xmm0 = L0*R2:L0*R0
pmuludq xmm1, xmm4 ; xmm1 = L1*R2:L1*R0
pmuludq xmm2, xmm4 ; xmm2 = L2*R2:L2*R0
pmuludq xmm3, xmm4 ; xmm3 = L3*R2:L3*R0
movdqa xmmword ptr[esp+0x00], xmm0 ;
movdqa xmmword ptr[esp+0x10], xmm1 ;
movdqa xmmword ptr[esp+0x20], xmm2 ;
movdqa xmmword ptr[esp+0x30], xmm3 ;
psrldq xmm4, 4 ; xmm4 = 00:R3:R2:R1
movdqa xmm1, xmmword ptr[ecx+0x10] ; xmm1 = L1:L1:L1:L1
movdqa xmm0, xmmword ptr[ecx+0x00] ; xmm0 = L0:L0:L0:L0
pmuludq xmm7, xmm4 ; xmm3 = L3*R3:L3*R1 ->V7-4
pmuludq xmm6, xmm4 ; xmm2 = L2*R3:L2*R1 ->V6-3
pmuludq xmm1, xmm4 ; xmm1 = L1*R3:L1*R1
pmuludq xmm0, xmm4 ; xmm0 = L0*R3:L0*R1
movdqa xmm5, xmm7 ;
movdqa xmm4, xmm6 ;
movdqa xmm3, xmm6 ;
pslldq xmm5, 4 ; xmm5 <<= 32
psrldq xmm4, 4 ; xmm4 >>= 32
pslldq xmm3, 12 ; L3*R3:L3*R1 ->V3-0
paddq xmm7, xmm4 ; L2*R3:L2*R1 ->V7-4
paddq xmm6, xmm5 ; L3*R3:L3*R1 ->V6-3
movdqa xmm2, xmm0 ;
movdqa xmm4, xmm0 ;
movdqa xmm5, xmm1 ;
psrldq xmm0, 8 ; xmm0 >>= 64
psrldq xmm1, 4 ; xmm1 >>= 32
pslldq xmm2, 8 ; L0*R3:L0*R1 ->V2-1H
pslldq xmm4, 4 ; xmm4 <<= 32
pslldq xmm5, 8 ; xmm5 <<= 64
paddq xmm6, xmm0 ; L0*R3:L0*R1 ->V6-3
paddq xmm3, xmm4 ; L0*R3:L0*R1 ->V3-0
paddq xmm6, xmm1 ; L1*R3:L1*R1 ->V6-3
paddq xmm3, xmm5 ; L1*R3:L1*R1 ->V3-0
psrldq xmm0, 4 ; xmm0 >>= 32
psrldq xmm1, 4 ; xmm1 >>= 32
pslldq xmm5, 4 ; xmm5 <<= 32
paddq xmm7, xmm0 ; L0*R3:L0*R1 ->V7-4
paddq xmm2, xmm5 ; L1*R3:L1*R1 ->V2-1H
paddq xmm7, xmm1 ; L1*R3:L1*R1 ->V7-4
movdqa xmm5, xmmword ptr[esp+0x30] ; xmm5 = L3*R2:L3*R0
movdqa xmm4, xmmword ptr[esp+0x20] ; xmm4 = L2*R2:L2*R0
movdqa xmm1, xmm5 ;
movdqa xmm0, xmm4 ;
paddq xmm6, xmm5 ; L3*R2:L3*R0 ->V6-3
psrldq xmm0, 4 ; xmm0 >>= 32
psrldq xmm1, 4 ; xmm1 >>= 32
pslldq xmm4, 8 ; xmm4 <<= 64
pslldq xmm5, 12 ; xmm5 <<= 96
paddq xmm6, xmm0 ; L2*R2:L2*R0 ->V6-3
paddq xmm3, xmm4 ; L2*R2:L2*R0 ->V3-0
paddq xmm7, xmm1 ; L3*R2:L3*R0 ->V7-4
paddq xmm3, xmm5 ; L3*R2:L3*R0 ->V3-0
psrldq xmm0, 4 ; xmm0 >>= 32
pslldq xmm4, 4 ; xmm4 <<= 32
movdqa xmm1, xmmword ptr[esp+0x00] ; xmm1 = L0*R2:L0*R0
movdqa xmm5, xmmword ptr[esp+0x10] ; xmm5 = L1*R2:L1*R0
paddq xmm7, xmm0 ; L2*R2:L2*R0 ->V7-4
paddq xmm2, xmm4 ; L2*R2:L2*R0 ->V2-1H
paddq xmm3, xmm1 ; L0*R2:L0*R0 ->V3-0
movdqa xmm0, xmm1 ;
movdqa xmm4, xmm5 ;
pslldq xmm1, 4 ; xmm1 <<= 32
pslldq xmm5, 4 ; xmm5 <<= 32
psrldq xmm0, 12 ; xmm0 >>= 96
psrldq xmm4, 8 ; xmm4 >>= 64
paddq xmm2, xmm1 ; L0*R2:L0*R0 ->V2-1H
paddq xmm6, xmm0 ; L0*R2:L0*R0 ->V6-3
paddq xmm3, xmm5 ; L1*R2:L1*R0 ->V3-0
paddq xmm6, xmm4 ; L1*R2:L1*R0 ->V6-3
pslldq xmm5, 4 ; xmm5 <<= 32
psrldq xmm4, 4 ; xmm4 >>= 32
paddq xmm2, xmm5 ; L1*R2:L1*R0 ->V2-1H
paddq xmm7, xmm4 ; L1*R2:L1*R0 ->V7-4
mov esp, eax ; now, V[1:0]=xmm3[63:0] ->OK
psrldq xmm2, 4 ; xmm2 >>= 32
movdqa xmm4, xmm3 ;
pcmpeqd xmm0, xmm0 ; xmm0 = FF:FF:FF:FF
pandn xmm4, xmm2 ;
psubd xmm2, xmm3 ; xmm2 -= xmm3
psrldq xmm0, 12 ; xmm0 = 00:00:00:FF
pslldq xmm4, 4 ; xmm4 <<= 32
psrld xmm4, 31 ; CF
pshufd xmm1, xmm0, 01001110b ; xmm1 = 00:FF:00:00
paddd xmm2, xmm4 ;
pand xmm2, xmm1 ;
paddq xmm3, xmm2 ; V[3:0]=xmm3 ->OK
movdqa xmmword ptr[ecx+0x00], xmm3 ; store V[3:0]
psrldq xmm3, 12 ; xmm3 >>= 96
psubd xmm3, xmm6 ; xmm3 -= xmm6
pand xmm3, xmm0 ;
paddq xmm6, xmm3 ; V[4:3]=xmm6[63:0] ->OK
pxor xmm1, xmm0 ; xmm1 = 00:FF:00:FF
psrldq xmm6, 4 ; V[4]=xmm6[31:0] ->OK
movdqa xmm2, xmm6 ; xmm2 = xmm6
psubd xmm6, xmm7 ; xmm6 -= xmm0
pand xmm6, xmm1 ;
paddq xmm7, xmm6 ; now, xmm7[63-0] bits OK!
psubusb xmm2, xmm7 ;
psrlq xmm2, 63 ; CF
pslldq xmm2, 8 ;
paddq xmm7, xmm2 ; xmm7[127-64] += CF
movdqa xmmword ptr[ecx+0x10], xmm7 ; store V[7:4]
ret;
}
}

复制代码

当前所有版本测试源代码及编译好的程序压缩包：

UInt128x128To256.zip (32.21 KB, 下载次数: 15)

无心人 · 发表于 2008-3-25 21:50:30

总感觉你这么安排指令会存在冲突

账号		自动登录	找回密码
密码			欢迎注册

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

对 58# 再改进