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

gxqcn

原帖由 无心人 于 2008-3-25 21:50 发表 总感觉你这么安排指令会存在冲突

我在写汇编时，尽量让相同指令聚集，且尽量让相邻语句无关，即减少依赖性； 但确实有许多靠经验、凭感觉的成分在组织安排指令顺序。 算法是经过仔细推敲论证过的。 不太明白你说的“冲突”是指什么？

无心人

比如连续4个乘 每个乘均是5延迟的 最多同时执行两条 那这4个就产生执行单元冲突和指令延时的耽误

gxqcn

ALU 指令讲究 U/V 通道的排布， 不知 SIMD 指令是否也有双通道问题？

无心人

当然啊 你没仔细读你推荐的那些PDF么？ 存在执行单元的限制啊 P4是5个执行单元，但并不是每个都能执行全部指令 每个都有固定的执行指令集合 能否并行是和执行单元数量有关的 指令在进入流水线解码就绪后是在一个地方集合等待 根据不同指令，要等不同执行单元空闲且操作数就绪后 才进入执行单元执行，最后在一个地方集合退出 每次执行单元的延迟数是不同的，虽然大部分都标1周期指令 但执行周期是不同的啊

无心人

一个函数如果指令数少 也许在解码后，可填满指令的cache 就是能保证执行前并不存在延时 但并不能保证执行冲突 1、执行单元限制 2、指令依赖 3、某些具体限制，比如AGI等

liangbch

写了一个 使用 MMX 寄存器，SSE2 乘 和 加 指令的版本，速度大幅领先于 现有的所有版本。 先贴上测试结果（PIV 2.6G 768M RAM），为了避免将函数调用放在后面容易占优势的嫌疑，将对该函数的调用放在前面。 Test function: UInt128x128To256_ANSI_C32(..) 10000000 times... Elapsed time: 2035.286 ms 92F84935 2DD0353C A7166445 8B52280D * ADBE65FA 33EC4DD6 F39E291C 6C1A072B = 63BF184A DE964DAC 8CD871A5 3FB4A4E7 D031F0A2 D4767BD7 A7C40427 3337152F Test function: UInt128x128To256_SSE2_76F(..) 10000000 times... Elapsed time: 395.657 ms 92F84935 2DD0353C A7166445 8B52280D * ADBE65FA 33EC4DD6 F39E291C 6C1A072B = 63BF184A DE964DAC 8CD871A5 3FB4A4E7 D031F0A2 D4767BD7 A7C40427 3337152F Test function: UInt128x128To256_SSE2_40F(..) 10000000 times... Elapsed time: 1027.907 ms 92F84935 2DD0353C A7166445 8B52280D * ADBE65FA 33EC4DD6 F39E291C 6C1A072B = 63BF184A DE964DAC 8CD871A5 3FB4A4E7 D031F0A2 D4767BD7 A7C40427 3337152F Test function: UInt128x128To256_SSE2_42F(..) 10000000 times... Elapsed time: 762.286 ms 92F84935 2DD0353C A7166445 8B52280D * ADBE65FA 33EC4DD6 F39E291C 6C1A072B = 63BF184A DE964DAC 8CD871A5 3FB4A4E7 D031F0A2 D4767BD7 A7C40427 3337152F Test function: UInt128x128To256_SSE2_54F(..) 10000000 times... Elapsed time: 734.089 ms 92F84935 2DD0353C A7166445 8B52280D * ADBE65FA 33EC4DD6 F39E291C 6C1A072B = 63BF184A DE964DAC 8CD871A5 3FB4A4E7 D031F0A2 D4767BD7 A7C40427 3337152F Test function: UInt128x128To256_SSE2_56F(..) 10000000 times... Elapsed time: 607.176 ms 92F84935 2DD0353C A7166445 8B52280D * ADBE65FA 33EC4DD6 F39E291C 6C1A072B = 63BF184A DE964DAC 8CD871A5 3FB4A4E7 D031F0A2 D4767BD7 A7C40427 3337152F Test function: UInt128x128To256_SSE2_58F(..) 10000000 times... Elapsed time: 923.039 ms 92F84935 2DD0353C A7166445 8B52280D * ADBE65FA 33EC4DD6 F39E291C 6C1A072B = 63BF184A DE964DAC 8CD871A5 3FB4A4E7 D031F0A2 D4767BD7 A7C40427 3337152F Test function: UInt128x128To256_SSE2_69F(..) 10000000 times... Elapsed time: 529.028 ms 92F84935 2DD0353C A7166445 8B52280D * ADBE65FA 33EC4DD6 F39E291C 6C1A072B = 63BF184A DE964DAC 8CD871A5 3FB4A4E7 D031F0A2 D4767BD7 A7C40427 3337152F 再贴上所有源代码（注意：该版有bug，作者已在 96# 修正！--gxqcn）

1. _declspec(naked) /* 注意：该版有bug，作者已在 96# 修正！*/
2. void UInt128x128To256_SSE2_76F( UINT32 * const result,
3. const UINT32 * const left,
4. const UINT32 * const right )
5. // 使用MMX寄存器， 但用SSE2指令
6. // 指令的使用，乘法 使用 32bit * 32bit 的MMX指令，加法使用 64 bit + 64 bit 的SSE2指令
8. // 16次乘法，分4趟， 操作数一律在mmx 寄存器中进行
9. // 第1趟计算 left[] * right[0]
10. // 第2趟计算 left[] * right[1]
11. // 第3趟计算 left[] * right[2]
12. // 第4趟计算 left[] * right[3]
14. // 每趟的计算结果为 5 DWORD，最低位DWORD 直接存储到result, 前3趟的计算结果的高4个DWORD（中间结果)存储到一组寄存器
16. // 寄存器的使用
18. // 中间结果 放入 mm0,mm1,mm2,mm3
20. // 掩码寄存器 mm7, 用来请64位寄存器的高32bit
21. // 进位 寄存器 和 当前乘积寄存器 为 mm4 和 mm5,轮换使用
22. // 乘数 寄存器 mm6
23. // 可以看到8个mmx寄存器全部派上用场
25. {
26. #define LEFT_REG eax
27. #define RIGHT_REG edx
28. #define RESULT_REG ecx
29. #define MASK_REG mm7
30. #define MUL_REG mm6
32. __asm
33. {
34. mov eax, 0xffffffff
35. mov RESULT_REG, dword ptr[esp + 04h] ; result
36. movd MASK_REG, eax
38. mov RIGHT_REG, dword ptr[esp + 0Ch] ; right
39. movq MUL_REG, [RIGHT_REG]
40. mov LEFT_REG, dword ptr [esp + 08h] ; left
43. //--------第1次乘
44. movd mm4, [LEFT_REG+0]
45. pmuludq mm4, MUL_REG
46. movd [RESULT_REG], mm4
47. psrlq mm4, 32 //carry
49. //------第2次乘
50. movd mm5, [LEFT_REG+4]
51. pmuludq mm5, MUL_REG
52. paddq mm5, mm4 //64bit + 32 bit carry
53. movq mm0, mm5
54. psrlq mm5, 32 //carry
55. pand mm0, MASK_REG //mm0,中间结果
57. //------第3次乘
58. movd mm4, [LEFT_REG+8]
59. pmuludq mm4, MUL_REG
60. paddq mm4, mm5 //64bit + 32 bit carry
61. movq mm1, mm4
62. psrlq mm4, 32 //carry
63. pand mm1, MASK_REG //mm1,中间结果
65. //------第4次乘
66. movd mm5, [LEFT_REG+12]
67. pmuludq mm5, MUL_REG
68. paddq mm5, mm4 //64bit + 32 bit carry
69. movq mm2, mm5
70. psrlq mm5, 32 //carry
71. pand mm2, MASK_REG //mm2,中间结果
72. movq mm3, mm5 //mm3,中间结果
74. //----得到乘数
75. psrlq MUL_REG, 32 //得到right[1]
78. //---第二趟乘---------
80. //--------第5次乘
81. movd mm4, [LEFT_REG+0]
82. pmuludq mm4, MUL_REG
83. paddq mm4, mm0 // 64bit + 32 bit 中间结果
84. movd [RESULT_REG+4], mm4
85. psrlq mm4, 32 //carry
87. //------第6次乘
88. movd mm5, [LEFT_REG+4]
89. pmuludq mm5, MUL_REG
90. paddq mm5, mm4 //64bit + 32 bit carry
91. paddq mm5, mm1 //64bit + 32 bit 中间结果
92. movq mm0, mm5
93. psrlq mm5, 32 //carry
94. pand mm0, MASK_REG //mm0,中间结果
96. //------第7次乘
97. movd mm4, [LEFT_REG+8]
98. pmuludq mm4, MUL_REG
99. paddq mm4, mm5 //64bit + 32 bit carry
100. paddq mm4, mm2 //64bit + 32 bit 中间结果
101. movq mm1, mm4
102. psrlq mm4, 32 //carry
103. pand mm1, MASK_REG //mm1,中间结果
105. //------第8次乘
106. movd mm5, [LEFT_REG+12]
107. pmuludq mm5, MUL_REG
108. paddq mm5, mm4 //64bit + 32 bit carry
109. paddq mm5, mm3 //64bit + 32 bit 中间结果
110. movq mm2, mm5
111. psrlq mm5, 32 //carry
112. pand mm2, MASK_REG //mm2,中间结果
113. movq mm3, mm5 //mm3,中间结果
115. //----得到乘数
116. movq MUL_REG, [RIGHT_REG+8]
118. //--------第9次乘
119. movd mm4, [LEFT_REG+0]
120. pmuludq mm4, MUL_REG
121. paddq mm4, mm0 // 64bit + 32 bit 中间结果
122. movd [RESULT_REG+8], mm4
123. psrlq mm4, 32 //carry
125. //------第10次乘
126. movd mm5, [LEFT_REG+4]
127. pmuludq mm5, MUL_REG
128. paddq mm5, mm4 //64bit + 32 bit carry
129. paddq mm5, mm1 //64bit + 32 bit 中间结果
130. movq mm0, mm5
131. psrlq mm5, 32 //carry
132. pand mm0, MASK_REG //mm0,中间结果
134. //------第11次乘
135. movd mm4, [LEFT_REG+8]
136. pmuludq mm4, MUL_REG
137. paddq mm4, mm5 //64bit + 32 bit carry
138. paddq mm4, mm2 //64bit + 32 bit 中间结果
139. movq mm1, mm4
140. psrlq mm4, 32 //carry
141. pand mm1, MASK_REG //mm1,中间结果
143. //------第12次乘
144. movd mm5, [LEFT_REG+12]
145. pmuludq mm5, MUL_REG
146. paddq mm5, mm4 //64bit + 32 bit carry
147. paddq mm5, mm3 //64bit + 32 bit 中间结果
148. movq mm2, mm5
149. psrlq mm5, 32 //carry
150. pand mm2, MASK_REG //mm2,中间结果
151. movq mm3, mm5 //mm3,中间结果
154. //----得到乘数
155. psrlq MUL_REG, 32 //得到right[1]
157. //--------第13次乘
158. movd mm4, [LEFT_REG+0]
159. pmuludq mm4, MUL_REG
160. paddq mm4, mm0 // 64bit + 32 bit 中间结果
161. movd [RESULT_REG+12], mm4
162. psrlq mm4, 32 //carry
164. //------第14次乘
165. movd mm5, [LEFT_REG+4]
166. pmuludq mm5, MUL_REG
167. paddq mm5, mm4 //64bit + 32 bit carry
168. paddq mm5, mm1 //64bit + 32 bit 中间结果
169. movd [RESULT_REG+16], mm5
170. psrlq mm5, 32 //carry
173. //------第15次乘
174. movd mm4, [LEFT_REG+8]
175. pmuludq mm4, MUL_REG
176. paddq mm4, mm5 //64bit + 32 bit carry
177. paddq mm4, mm2 //64bit + 32 bit 中间结果
178. movd [RESULT_REG+20], mm4
179. psrlq mm4, 32 //carry
182. //------第16次乘
183. movd mm5, [LEFT_REG+12]
184. pmuludq mm5, MUL_REG
185. paddq mm5, mm4 //64bit + 32 bit carry
186. paddq mm5, mm3 //64bit + 32 bit 中间结果
187. movq [RESULT_REG+24], mm5
188. psrlq mm5, 32 //carry
190. //------保存最高DWORD
191. movq [RESULT_REG+28], mm5 //mm3,中间结果
193. //emms
194. ret
195. }
196. }

复制代码

我的指令数为111条，指令数可能是所有版本中最小的。

gxqcn

在我的 56# 中的指令数为 110 行（包含最后的那个 ret 指令），比ls的更少。 因为 mmx 寄存器与 fpu 寄存器物理上是共享的， 所以 ls 若将 emms 指令屏蔽掉，将导致程序后面无法正常切换到浮点单元，只不过测试程序通篇全整型，不存在这个问题罢了。 不知我这款 AMD CPU 是 MMX 优化得太好，还是对 SSE2 支持得太烂，测试结果如下：

UInt128x128To256_ANSI_C32(..): 640.631 ms UInt128x128To256_SSE2_40F(..): 516.867 ms UInt128x128To256_SSE2_42F(..): 810.261 ms UInt128x128To256_SSE2_54F(..): 755.762 ms UInt128x128To256_SSE2_56F(..): 646.375 ms UInt128x128To256_SSE2_58F(..): 639.292 ms UInt128x128To256_SSE2_69F(..): 633.678 ms UInt128x128To256_SSE2_76F(..): 379.350 ms

在我这台机器上，76F 替代了 40F 坐上了头把交椅 现在，非常期待 Core 2 系列的测试结果，以了解CPU对指令集的优化程度和发展趋势。。。

无心人

你厉害 谁帖Core2测试结果啊？

无心人

P4 2.0测试结果 P4E P4D估计不会出现其他结果 验证了我说的MMX寄存器在P4上快的结论 现在差Core结果 40 1794 42 1907 54 1894 56 894 58 756 76 642 76带不带emms均超其他函数

liangbch

准备再写一个完全的SSE2版的，主要思路： 1.一次SSE2乘法指令计算 2个 32bit 乘以 32bit 的积 2.累加使用SSE2 128bit 加指令,但高64bit不使用 3.和76楼一样,使用4个XMM寄存器保存中间结果 4.依然需要使用掩码寄存器清除高bit,但使用 掩码清除高位的指令的使用次数大大减少,只需要3次. 5.依然是8大XMM寄存器其上阵,最大限度利用系统资源. 6.预计指令数将更少,但不知速度能否超越76楼.