karatsuba乘法两种算法比较测试

liangbch

我的做法： 两个n位的数相加，结果可能出现进位，这时可用一个变量carry 表示进位(carry=0 或者 carry=1)，一个长为n的数组表示结果。 两个n位的数相减，结果可能出现借位，如果出现借位,carry=-1， 结果采用补码表示，数组的各个元素使用补码表示。

liangbch

下面贴出GMP kara乘法 的代码：来自GMP 4.22 中的 mpn\generic\mul_n.c ，GMP 4.22 可从 http://swox.com/gmp/ 下载 顺便说一下，我的代码比他的简单，而且是基于2^30进制的。

1. /* Multiplies using 3 half-sized mults and so on recursively.
2. * p[0..2*n-1] := product of a[0..n-1] and b[0..n-1].
3. * No overlap of p[...] with a[...] or b[...].
4. * ws is workspace.
5. */
7. void
8. mpn_kara_mul_n (mp_ptr p, mp_srcptr a, mp_srcptr b, mp_size_t n, mp_ptr ws)
9. {
10. mp_limb_t w, w0, w1;
11. mp_size_t n2;
12. mp_srcptr x, y;
13. mp_size_t i;
14. int sign;
16. n2 = n >> 1;
17. ASSERT (n2 > 0);
19. if ((n & 1) != 0)
20. {
21. /* Odd length. */
22. mp_size_t n1, n3, nm1;
24. n3 = n - n2;
26. sign = 0;
27. w = a[n2];
28. if (w != 0)
29. w -= mpn_sub_n (p, a, a + n3, n2);
30. else
31. {
32. i = n2;
33. do
34. {
35. --i;
36. w0 = a[ i ];
37. w1 = a[n3 + i];
38. }
39. while (w0 == w1 && i != 0);
40. if (w0 < w1)
41. {
42. x = a + n3;
43. y = a;
44. sign = ~0;
45. }
46. else
47. {
48. x = a;
49. y = a + n3;
50. }
51. mpn_sub_n (p, x, y, n2);
52. }
53. p[n2] = w;
55. w = b[n2];
56. if (w != 0)
57. w -= mpn_sub_n (p + n3, b, b + n3, n2);
58. else
59. {
60. i = n2;
61. do
62. {
63. --i;
64. w0 = b[ i ];
65. w1 = b[n3 + i];
66. }
67. while (w0 == w1 && i != 0);
68. if (w0 < w1)
69. {
70. x = b + n3;
71. y = b;
72. sign = ~sign;
73. }
74. else
75. {
76. x = b;
77. y = b + n3;
78. }
79. mpn_sub_n (p + n3, x, y, n2);
80. }
81. p[n] = w;
83. n1 = n + 1;
84. if (n2 < MUL_KARATSUBA_THRESHOLD)
85. {
86. if (n3 < MUL_KARATSUBA_THRESHOLD)
87. {
88. mpn_mul_basecase (ws, p, n3, p + n3, n3);
89. mpn_mul_basecase (p, a, n3, b, n3);
90. }
91. else
92. {
93. mpn_kara_mul_n (ws, p, p + n3, n3, ws + n1);
94. mpn_kara_mul_n (p, a, b, n3, ws + n1);
95. }
96. mpn_mul_basecase (p + n1, a + n3, n2, b + n3, n2);
97. }
98. else
99. {
100. mpn_kara_mul_n (ws, p, p + n3, n3, ws + n1);
101. mpn_kara_mul_n (p, a, b, n3, ws + n1);
102. mpn_kara_mul_n (p + n1, a + n3, b + n3, n2, ws + n1);
103. }
105. if (sign)
106. mpn_add_n (ws, p, ws, n1);
107. else
108. mpn_sub_n (ws, p, ws, n1);
110. nm1 = n - 1;
111. if (mpn_add_n (ws, p + n1, ws, nm1))
112. {
113. mp_limb_t x = (ws[nm1] + 1) & GMP_NUMB_MASK;
114. ws[nm1] = x;
115. if (x == 0)
116. ws[n] = (ws[n] + 1) & GMP_NUMB_MASK;
117. }
118. if (mpn_add_n (p + n3, p + n3, ws, n1))
119. {
120. mpn_incr_u (p + n1 + n3, 1);
121. }
122. }
123. else
124. {
125. /* Even length. */
126. i = n2;
127. do
128. {
129. --i;
130. w0 = a[ i ];
131. w1 = a[n2 + i];
132. }
133. while (w0 == w1 && i != 0);
134. sign = 0;
135. if (w0 < w1)
136. {
137. x = a + n2;
138. y = a;
139. sign = ~0;
140. }
141. else
142. {
143. x = a;
144. y = a + n2;
145. }
146. mpn_sub_n (p, x, y, n2);
148. i = n2;
149. do
150. {
151. --i;
152. w0 = b[ i ];
153. w1 = b[n2 + i];
154. }
155. while (w0 == w1 && i != 0);
156. if (w0 < w1)
157. {
158. x = b + n2;
159. y = b;
160. sign = ~sign;
161. }
162. else
163. {
164. x = b;
165. y = b + n2;
166. }
167. mpn_sub_n (p + n2, x, y, n2);
169. /* Pointwise products. */
170. if (n2 < MUL_KARATSUBA_THRESHOLD)
171. {
172. mpn_mul_basecase (ws, p, n2, p + n2, n2);
173. mpn_mul_basecase (p, a, n2, b, n2);
174. mpn_mul_basecase (p + n, a + n2, n2, b + n2, n2);
175. }
176. else
177. {
178. mpn_kara_mul_n (ws, p, p + n2, n2, ws + n);
179. mpn_kara_mul_n (p, a, b, n2, ws + n);
180. mpn_kara_mul_n (p + n, a + n2, b + n2, n2, ws + n);
181. }
183. /* Interpolate. */
184. if (sign)
185. w = mpn_add_n (ws, p, ws, n);
186. else
187. w = -mpn_sub_n (ws, p, ws, n);
188. w += mpn_add_n (ws, p + n, ws, n);
189. w += mpn_add_n (p + n2, p + n2, ws, n);
190. MPN_INCR_U (p + n2 + n, 2 * n - (n2 + n), w);
191. }
192. }

复制代码

liangbch

下面贴出NTL 中的kara乘法 源代码，kar_mul函数来自 WinNTL-5_4_2\src\c_lip_impl.h，WinNTL可从 http://www.shoup.net下载

1. static
2. void kar_mul(long *c, long *a, long *b, long *stk)
3. {
4. long sa, sb, sc;
6. if (*a < *b) { long *t = a; a = b; b = t; }
8. sa = *a;
9. sb = *b;
10. sc = sa + sb;
12. if (sb < KARX) {
13. /* classic algorithm */
15. long *pc, i, *pb;
17. pc = c;
18. for (i = sc; i; i--) {
19. pc++;
20. *pc = 0;
21. }
23. pc = c;
24. pb = b;
25. for (i = sb; i; i--) {
26. pb++;
27. pc++;
28. zaddmul(*pb, pc, a);
29. }
30. }
31. else {
32. long hsa = (sa + 1) >> 1;
34. if (hsa < sb) {
35. /* normal case */
37. long *T1, *T2, *T3;
39. /* allocate space */
41. T1 = stk; stk += hsa + 2;
42. T2 = stk; stk += hsa + 2;
43. T3 = stk; stk += (hsa << 1) + 3;
45. if (stk-kmem > max_kmem) zhalt("internal error: kmem overflow");
47. /* compute T1 = a_lo + a_hi */
49. kar_fold(T1, a, hsa);
51. /* compute T2 = b_lo + b_hi */
53. kar_fold(T2, b, hsa);
55. /* recursively compute T3 = T1 * T2 */
57. kar_mul(T3, T1, T2, stk);
59. /* recursively compute a_hi * b_hi into high part of c */
60. /* and subtract from T3 */
62. {
63. long olda, oldb;
65. olda = a[hsa]; a[hsa] = sa-hsa;
66. oldb = b[hsa]; b[hsa] = sb-hsa;
68. kar_mul(c + (hsa << 1), a + hsa, b + hsa, stk);
69. kar_sub(T3, c + (hsa << 1));
71. a[hsa] = olda;
72. b[hsa] = oldb;
73. }
75. /* recursively compute a_lo*b_lo into low part of c */
76. /* and subtract from T3 */
78. *a = hsa;
79. *b = hsa;
81. kar_mul(c, a, b, stk);
82. kar_sub(T3, c);
84. *a = sa;
85. *b = sb;
87. /* finally, add T3 * NTL_RADIX^{hsa} to c */
89. kar_add(c, T3, hsa);
90. }
91. else {
92. /* degenerate case */
94. long *T;
96. T = stk; stk += sb + hsa + 1;
98. if (stk-kmem > max_kmem) zhalt("internal error: kmem overflow");
100. /* recursively compute b*a_hi into high part of c */
101. {
102. long olda;
104. olda = a[hsa]; a[hsa] = sa-hsa;
105. kar_mul(c + hsa, a + hsa, b, stk);
106. a[hsa] = olda;
107. }
109. /* recursively compute b*a_lo into T */
111. *a = hsa;
112. kar_mul(T, a, b, stk);
113. *a = sa;
115. /* fix-up result */
117. kar_fix(c, T, hsa);
118. }
119. }
121. /* normalize result */
123. while (c[sc] == 0 && sc > 1) sc--;
124. *c = sc;
125. }

复制代码

无心人

两个采取了不同的策略啊

gxqcn

原帖由 无心人 于 2008-4-28 09:12 发表 是否存在比较的必要？ 考虑工作量 1、前加法2k 乘法(k+1)^2 后加法2k 总减法 2k+1 2、前减法2k 纠正2k+2 乘法k^2 后加法2k 总加法 2k+1 差异在 1、乘法(k+1)^2 2、纠正2k+2 乘法k^2 (k+1)^2 - k^2 - (2k ...

其实，前加法也并非一定会出现进位。

gxqcn

应该说，直接采用楼主的公式还不是最优化的， 更好的形式可以进一步减少一次 大数的加/减法 运算。 下面的公式是从我的 HugeCalc 源代码注释中抠出的：

// ( A1*b + A0 )( B1*b + B0 ) // = (( A1*B1*b - A0*B0 ) + ( A0 + A1 )( B0 + B1 ))*b // - ( A1*B1*b - A0*B0 ) // ( A1*b + A0 )( B1*b + B0 ) // = (( A1*B1*b + A0*B0 ) - ( A0 - A1 )( B0 - B1 ))*b // + ( A1*B1*b + A0*B0 )

计算400000!（屏蔽Toom算法及卷积算法），前者需 5.919233s，后者需 6.188290s 也就是说，在 HugeCalc 上的表现，用第一种形式可比第二种至少快 4% 左右。

liangbch

我公布一下我的kara乘法的 技术数据： 1.因为我的加减法接口是这样的 int DEC_Add(DWORD *a,const DWORD *b,size_t len) 故在做加减法是可使用更少的寄存器，访问更少的内存。缺点应用到kara_mul时，需要使用内存拷贝函数. 2.若计算2个长度为 2*len的数的乘积，则在每次递归调用kara_mul时，除了乘法外，需要附加的操作有： 1. 2 * len 次整数减法。 2. len + len + 4*len 次整数复制 3. 3 * (len*2) 次整数加法 4. 2 * len * 0.5 次整数求补 (大约有50%的概率需要做求补运算) 共计：2×len 次减法，6×len 次整数加法，6×len次整数复制，len次 求补运算。 gxq可否给出你的kara乘法中，除了乘法外的其他运算的次数。

gxqcn

// ( A1*b + A0 )( B1*b + B0 ) // = (( A1*B1*b -^[1] A0*B0 ) +^[5] ( A0 +^[3] A1 )( B0 +^[4] B1 ))*b // -^[6] ( A1*B1*b - A0*B0 )^[2]

[1] 2*len 次整数减法 - A0*B0 [2] 3*len 次整数复制 ( A1*B1*b - A0*B0 ) [3] len 次整数加法 + A1 [4] len 次整数加法 + B1 [5] 2*len 次整数加法 + ( A0 + A1 )( B0 + B1 ) [6] 3*len 次整数减法 - ( A1*B1*b - A0*B0 ) 注：选择好适当大小的 b，可以确保上述所有减法无须修正（即结果一定为正）。 合计：加减法 9*len 单元次；复制 3*len 单元次；乘法 3*len*len 单元次；移位 无（在相关运算中直接定位即可）

无心人

GxQ 你没考虑过你的加减法和求负还有内存拷贝都不是最优的 加减存在一次处理4X32位操作的算法 求负，拷贝存在一次处理4X128位操作的算法

gxqcn

那需要用户CPU至少得支持 SSE2 才行。 而有许多老爷机是不具备的。