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 )
6. {
7. #undef LEFT_REG
8. #undef RIGHT_REG
9. #undef RESULT_REG
11. #define LEFT_REG esi
12. #define RIGHT_REG edi
13. #define RESULT_REG ebp
15. __asm
16. {
17. push esi
18. push edi
19. push ebx
20. push ebp
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
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
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
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
49. mov dword ptr [RESULT_REG+8], ecx
50. mov ecx, ebx
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
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
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
73. mov dword ptr [RESULT_REG+12], ecx
74. mov ecx, ebx
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
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
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
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
104. mov dword ptr [RESULT_REG+16], ecx
105. mov ecx, ebx
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
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
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
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
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
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
152. mov dword ptr [RESULT_REG+24], ecx
153. mov dword ptr [RESULT_REG+28], ebx
156. pop ebp
157. pop ebx
158. pop edi
159. pop esi
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%。

无心人

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