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

无心人 发表于 2008-3-31 11:37:31

那你遇到的第一个问题就是寄存器不足
:lol

liangbch 发表于 2008-3-31 13:15:02

寄存器不足就得使用RAM了，其实，如果这个问题改为 256bit × 256bit，所有的版本都会遇到寄存器不足的问题。
下面是一个完全不使用任何MMX,SSE,SSE2指令的版本,共102条指令。在我的电脑上进行测试，结果表明，除了76楼（fixed bug后为96楼）的程序，比这个程序快2倍以上以外，其他的版本的加速程度都无法达到2倍。
先贴出测试结果。
Test function: UInt128x128To256_ANSI_C32(..) 10000000 times...
Elapsed time: 1905.637 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_SSE2_40F(..) 10000000 times...
Elapsed time: 706.421 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_SSE2_42F(..) 10000000 times...
Elapsed time: 729.128 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_SSE2_54F(..) 10000000 times...
Elapsed time: 695.204 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_SSE2_56F(..) 10000000 times...
Elapsed time: 580.515 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_SSE2_58F(..) 10000000 times...
Elapsed time: 880.662 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_SSE2_69F(..) 10000000 times...
Elapsed time: 834.156 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_SSE2_94F(..) 10000000 times...
Elapsed time: 864.124 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_SSE2_96F(..) 10000000 times...
Elapsed time: 403.787 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

Test function: UInt128x128To256_ALU_102F(..) 10000000 times...
Elapsed time: 945.073 ms
EAAC506E 17E8573A 5A8A540B ACAE5D9E * C30C395E 93C433B8 4AC06050 7E9A50E8
= B2CC75EE B201F29B 1CDA40BB 312C399D A214AF3C 4BD8278E D1BC04E3 10523730

再贴出完整的代码._declspec(naked)
void UInt128x128To256_ALU_102F( UINT32 * const result,
                           const UINT32 * const left,
                           const UINT32 * const right )

{
#undef LEFT_REG
#undef RIGHT_REG
#undef RESULT_REG

#define LEFT_REG esi
#define RIGHT_REG edi
#define RESULT_REG ebp

__asm
{
   push esi
   push edi
   push ebx
   push ebp

   mov RIGHT_REG, dword ptr ;right
   mov LEFT_REG,dword ptr ; left
   mov RESULT_REG, dword ptr ; result
   xor ebx, ebx
   xor ecx, ecx

   // result，一次乘法
   mov eax,dword ptr
   mul dword ptr
   mov dword ptr , eax
   mov dword ptr , edx

   // result，2次乘法
   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov dword ptr , ecx
   mov ecx, ebx

   // result, 3 次乘法
   xor ebx, ebx
   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx,edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov dword ptr , ecx
   mov ecx, ebx

   // result，4次乘法
   xor ebx, ebx
   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov dword ptr , ecx
   mov ecx, ebx

   // result，3次乘法
   xor ebx, ebx
   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0
   mov dword ptr , ecx
   mov ecx, ebx

   // result，2次乘法
   xor ebx, ebx
   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   mov eax,dword ptr
   mul dword ptr
   add dword ptr , eax
   adc ecx, edx
   adc ebx,0

   // result，1次乘法
   mov eax,dword ptr
   mul dword ptr
   add ecx, eax
   adc ebx,edx

   mov dword ptr , ecx
   mov dword ptr , ebx

   pop ebp
   pop ebx
   pop edi
   pop esi

   ret
}
}

无心人 发表于 2008-3-31 13:48:49

可以做到不用额外单元就不错了
===========================
256X256似乎没必要写了

只要用MMX寄存器或者SSE寄存器
可以做到不用内存保存临时变量的
除非超过2^31个双字乘
=============================
之所以加速比达不到2，关键是乘法本身速度慢
进位加法复杂

liangbch 发表于 2008-3-31 13:56:04

我的代码只用了四个push 语句保存一些寄存器的旧的值。除此以外，并没有使用额外的临时变量。

无心人 发表于 2008-3-31 21:42:51

折腾来折腾去
还是Knuth算法最好

无心人 发表于 2008-3-31 21:49:09

voidAsmMulLL(unsigned long*pL, unsigned long*pR, unsigned long*pA, unsigned long tL, unsigned long tR)
{
if ((tL == 0) || (tR == 0))
return;
__asm       {
         mov ecx, tL
         movesi, dword ptr
         movedi, dword ptr
         movebx, dword ptr
      pxor mm3, mm3

mbinmul2:
      mov edx, ecx
      mov eax, ebx
         pxor mm0, mm0
      mov ecx, tR
      movd mm1, dword ptr
      movd mm4, edi
mbinmul3:
         movd mm2, dword ptr
         lea edi,
      movd mm3, dword ptr
         pmuludq mm2, mm1
      paddq mm0, mm3
         paddq mm0, mm2
      movd dword ptr , mm0
         psrlq mm0, 32
         lea ebx,
         loop mbinmul3
      movd edi, mm4
      movd dword ptr , mm0
         mov ebx, eax
         lea esi,
         lea ebx,
         mov ecx, edx
         loop mbinmul2
      emms
         }
}
做个备份
到学校去测试

无心人 发表于 2008-4-1 08:12:36

void UInt128x128To256_SSE2_107F( UINT32 * const result,
                           const UINT32 * const left,
                           const UINT32 * const right )
{
__asm       {
         mov ecx, 4
         movesi, dword ptr
         movedi, dword ptr
         movebx, dword ptr
         pxor xmm0, xmm0
         pxor xmm1, xmm1
         movdqa xmmword ptr , xmm0
         movdqa xmmword ptr , xmm1
         pxor mm3, mm3
mbinmul2:
         mov edx, ecx
         mov eax, ebx
         pxor mm0, mm0
         mov ecx,4
         movd mm1, dword ptr
         movd mm4, edi
mbinmul3:
         movd mm2, dword ptr
         lea edi,
         movd mm3, dword ptr //这里存在问题
         pmuludq mm2, mm1
         paddq mm0, mm3
         paddq mm0, mm2
         movd dword ptr , mm0
         psrlq mm0, 32
         lea ebx,
         loop mbinmul3
         movd edi, mm4
         movd dword ptr , mm0
         mov ebx, eax
         lea esi,
         lea ebx,
         mov ecx, edx
         loop mbinmul2
         emms
         }
}

改成本题要求看下时间：）
//在我这里时间和40差不多

无心人 发表于 2008-4-1 10:38:18

//各种优化均尝试了，lea改add，loop改双指令，均无效
//edi压栈保存造成时间增加
但如果改普遍的，需额外执行一个结果清零的函数
目前问题
1、是否能把edi保存在常规寄存器里，就是再节约一个寄存器
2、能否在不增加很多代码情况下，去掉结果预清零过程
解决这两个问题，俺的B计划之长乘将提前出笼

假设m个双字乘以n个双子
该算法的执行指令数约等于13m + 10mn

liangbch 发表于 2008-4-1 11:26:47

楼主可能是想找出一个解决长乘算法的完美算法。这基本是徒劳的。

最根本的解决之道是：大数的表示法不能采用232进制，而转为采用 2^31或者 2^30进制，即采用每个DWORD存储31位或者30位2进制数，虽然浪费了一点存储空间，也带来乘法次数的增多（增加了大约1/15多一点），但其优点可完全抵消此缺点。并对算法的设计带来实质性的帮助。
对于使用ALU指令来说，可大大降低ADC指令的运行次数。
对于使用SSE2指令来说，可在ALU指令的基础上真正提速2倍，充分发挥SSE2指令 1次可计算2个32bit×32bit和 64bbit + 64bit的能力。平均而言，每次30 bit × 30bit的运算仅需要使用3条指令。

举例来说，以那个128bit*128bit 的函数为例，需要 4×4次 32bit × 32bit 的乘法，平均每次 32bit × 32bit 的乘需要使用 6条以上的指令，以此推算，一个15DWORD( 15×32bit） × 15DWORD 的乘法需要 15×15×6= 1350次指令。而依照本方法，采用2^30进制，一个15×32bit的数需要 16个DWORD 存储，故采用此方式存储，需要16×16次 30bit* 30bit的乘法，所以需要 16*16*3=768条指令，相比前者的1350条指令，单指令数而言，减少至先前的57%，若以指令数和运行时间成反比（即假定每条指令的时间均相同），则此法可提速75%。

无心人 发表于 2008-4-1 11:42:04

那你并没看到解包运算的麻烦吧
你怎么解释不在字节倍数上断开
而造成结果后处理的工作量？
===========================
而且，不光是乘法
加法也是常见运算
这么表示并不能在加上取得任何效果

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

数学研发论坛's Archiver