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

无心人

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

liangbch

　　寄存器不足就得使用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

再贴出完整的代码.

1. _declspec(naked)
   
2. void UInt128x128To256_ALU_102F( UINT32 * const result,
   
3.                                 const UINT32 * const left,
   
4.                                 const UINT32 * const right )
   
5. 
   
6. {
   
7.      #undef LEFT_REG
   
8.      #undef RIGHT_REG
   
9.      #undef RESULT_REG
   
10. 
    
11.      #define LEFT_REG        esi
    
12.      #define RIGHT_REG        edi
    
13.      #define RESULT_REG        ebp
    
14. 
    
15.      __asm
    
16.     {
    
17.         push        esi
    
18.         push        edi
    
19.         push        ebx
    
20.         push        ebp
    
21. 
    
22. 
    
23.         mov        RIGHT_REG, dword ptr[esp +  0Ch+16]   ;  right
    
24.         mov        LEFT_REG,  dword ptr [esp + 08h+16]   ; left
    
25.         mov        RESULT_REG, dword ptr[esp + 04h+16]   ; result
    
26.         xor        ebx, ebx
    
27.         xor        ecx, ecx
    
28. 
    
29.         // result[0]，一次乘法
    
30.         mov        eax,dword ptr [LEFT_REG]
    
31.         mul        dword ptr [RIGHT_REG]
    
32.         mov        dword ptr [RESULT_REG], eax
    
33.         mov        dword ptr [RESULT_REG+4], edx
    
34. 
    
35. 
    
36.         // result[1]，2次乘法
    
37.         mov        eax,dword ptr [LEFT_REG]
    
38.         mul        dword ptr [RIGHT_REG+4]
    
39.         add        dword ptr [RESULT_REG+4], eax
    
40.         adc        ecx, edx
    
41.         adc        ebx,0
    
42. 
    
43.         mov        eax,dword ptr [LEFT_REG+4]
    
44.         mul        dword ptr [RIGHT_REG]
    
45.         add        dword ptr [RESULT_REG+4], eax
    
46.         adc        ecx, edx
    
47.         adc        ebx,0
    
48. 
    
49.         mov        dword ptr [RESULT_REG+8], ecx
    
50.         mov        ecx, ebx
    
51. 
    
52. 
    
53.         // result[2], 3 次乘法
    
54.         xor        ebx, ebx
    
55.         mov        eax,dword ptr [LEFT_REG]
    
56.         mul        dword ptr [RIGHT_REG+8]
    
57.         add        dword ptr [RESULT_REG+8], eax
    
58.         adc        ecx,edx
    
59.         adc        ebx,0
    
60. 
    
61.         mov        eax,dword ptr [LEFT_REG+4]
    
62.         mul        dword ptr [RIGHT_REG+4]
    
63.         add        dword ptr [RESULT_REG+8], eax
    
64.         adc        ecx, edx
    
65.         adc        ebx,0
    
66. 
    
67.         mov        eax,dword ptr [LEFT_REG+8]
    
68.         mul        dword ptr [RIGHT_REG+0]
    
69.         add        dword ptr [RESULT_REG+8], eax
    
70.         adc        ecx, edx
    
71.         adc        ebx,0
    
72. 
    
73.         mov        dword ptr [RESULT_REG+12], ecx
    
74.         mov        ecx, ebx
    
75. 
    
76. 
    
77.         // result[3]，4次乘法
    
78.         xor        ebx, ebx
    
79.         mov        eax,dword ptr [LEFT_REG+0]
    
80.         mul        dword ptr [RIGHT_REG+12]
    
81.         add        dword ptr [RESULT_REG+12], eax
    
82.         adc        ecx, edx
    
83.         adc        ebx,0
    
84. 
    
85.         mov        eax,dword ptr [LEFT_REG+4]
    
86.         mul        dword ptr [RIGHT_REG+8]
    
87.         add        dword ptr [RESULT_REG+12], eax
    
88.         adc        ecx, edx
    
89.         adc        ebx,0
    
90. 
    
91.         mov        eax,dword ptr [LEFT_REG+8]
    
92.         mul        dword ptr [RIGHT_REG+4]
    
93.         add        dword ptr [RESULT_REG+12], eax
    
94.         adc        ecx, edx
    
95.         adc        ebx,0
    
96. 
    
97. 
    
98.         mov        eax,dword ptr [LEFT_REG+12]
    
99.         mul        dword ptr [RIGHT_REG+0]
    
100.         add        dword ptr [RESULT_REG+12], eax
     
101.         adc        ecx, edx
     
102.         adc        ebx,0
     
103. 
     
104.         mov        dword ptr [RESULT_REG+16], ecx
     
105.         mov        ecx, ebx
     
106. 
     
107. 
     
108.         // result[4]，3次乘法
     
109.         xor        ebx, ebx
     
110.         mov        eax,dword ptr [LEFT_REG+4]
     
111.         mul        dword ptr [RIGHT_REG+12]
     
112.         add        dword ptr [RESULT_REG+16], eax
     
113.         adc        ecx, edx
     
114.         adc        ebx,0
     
115. 
     
116.         mov        eax,dword ptr [LEFT_REG+8]
     
117.         mul        dword ptr [RIGHT_REG+8]
     
118.         add        dword ptr [RESULT_REG+16], eax
     
119.         adc        ecx, edx
     
120.         adc        ebx,0
     
121. 
     
122.         mov        eax,dword ptr [LEFT_REG+12]
     
123.         mul        dword ptr [RIGHT_REG+4]
     
124.         add        dword ptr [RESULT_REG+16], eax
     
125.         adc        ecx, edx
     
126.         adc        ebx,0
     
127.         mov        dword ptr [RESULT_REG+20], ecx
     
128.         mov        ecx, ebx
     
129. 
     
130. 
     
131.         // result[5]，2次乘法
     
132.         xor        ebx, ebx
     
133.         mov        eax,dword ptr [LEFT_REG+8]
     
134.         mul        dword ptr [RIGHT_REG+12]
     
135.         add        dword ptr [RESULT_REG+20], eax
     
136.         adc        ecx, edx
     
137.         adc        ebx,0
     
138. 
     
139.         mov        eax,dword ptr [LEFT_REG+12]
     
140.         mul        dword ptr [RIGHT_REG+8]
     
141.         add        dword ptr [RESULT_REG+20], eax
     
142.         adc        ecx, edx
     
143.         adc        ebx,0
     
144. 
     
145. 
     
146.         // result[6，7]，1次乘法
     
147.         mov        eax,dword ptr [LEFT_REG+12]
     
148.         mul        dword ptr [RIGHT_REG+12]
     
149.         add        ecx, eax
     
150.         adc        ebx,edx
     
151. 
     
152.         mov        dword ptr [RESULT_REG+24], ecx
     
153.         mov        dword ptr [RESULT_REG+28], ebx
     
154. 
     
155. 
     
156.         pop        ebp
     
157.         pop        ebx
     
158.         pop        edi
     
159.         pop        esi
     
160. 
     
161.         ret
     
162.      }
     
163. }

复制代码

无心人

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

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

liangbch

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

无心人

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

无心人

void  AsmMulLL(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
             mov  esi, dword ptr [pL]
             mov  edi, dword ptr [pR]
             mov  ebx, dword ptr [pA]
             pxor mm3, mm3

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

无心人

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

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

无心人

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

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

liangbch

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

   最根本的解决之道是：大数的表示法不能采用2³²进制，而转为采用 $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%。

无心人

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