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. */
   
6. 
   
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;
    
15. 
    
16.   n2 = n >> 1;
    
17.   ASSERT (n2 > 0);
    
18. 
    
19.   if ((n & 1) != 0)
    
20.     {
    
21.       /* Odd length. */
    
22.       mp_size_t n1, n3, nm1;
    
23. 
    
24.       n3 = n - n2;
    
25. 
    
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;
    
54. 
    
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;
    
82. 
    
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.         }
     
104. 
     
105.       if (sign)
     
106.         mpn_add_n (ws, p, ws, n1);
     
107.       else
     
108.         mpn_sub_n (ws, p, ws, n1);
     
109. 
     
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);
     
147. 
     
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);
     
168. 
     
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.         }
     
182. 
     
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;
   
5. 
   
6.    if (*a < *b) { long *t = a; a = b; b = t; }
   
7. 
   
8.    sa = *a;
   
9.    sb = *b;
   
10.    sc = sa + sb;
    
11. 
    
12.    if (sb < KARX) {
    
13.       /* classic algorithm */
    
14. 
    
15.       long *pc, i, *pb;
    
16. 
    
17.       pc = c;
    
18.       for (i = sc; i; i--) {
    
19.          pc++;
    
20.          *pc = 0;
    
21.       }
    
22.    
    
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;
    
33. 
    
34.       if (hsa < sb) {
    
35.          /* normal case */
    
36. 
    
37.          long *T1, *T2, *T3;
    
38. 
    
39.          /* allocate space */
    
40. 
    
41.          T1 = stk;  stk += hsa + 2;
    
42.          T2 = stk;  stk += hsa + 2;
    
43.          T3 = stk;  stk += (hsa << 1) + 3;
    
44. 
    
45.          if (stk-kmem > max_kmem) zhalt("internal error: kmem overflow");
    
46. 
    
47.          /* compute T1 = a_lo + a_hi */
    
48. 
    
49.          kar_fold(T1, a, hsa);
    
50. 
    
51.          /* compute T2 = b_lo + b_hi */
    
52. 
    
53.          kar_fold(T2, b, hsa);
    
54.          
    
55.          /* recursively compute T3 = T1 * T2 */
    
56. 
    
57.          kar_mul(T3, T1, T2, stk);
    
58. 
    
59.          /* recursively compute a_hi * b_hi into high part of c */
    
60.          /* and subtract from T3 */
    
61. 
    
62.          {
    
63.             long olda, oldb;
    
64. 
    
65.             olda = a[hsa];  a[hsa] = sa-hsa;
    
66.             oldb = b[hsa];  b[hsa] = sb-hsa;
    
67. 
    
68.             kar_mul(c + (hsa << 1), a + hsa, b + hsa, stk);
    
69.             kar_sub(T3, c + (hsa << 1));
    
70. 
    
71.             a[hsa] = olda;
    
72.             b[hsa] = oldb;
    
73.          }
    
74. 
    
75.          /* recursively compute a_lo*b_lo into low part of c */
    
76.          /* and subtract from T3 */
    
77. 
    
78.          *a = hsa;
    
79.          *b = hsa;
    
80. 
    
81.          kar_mul(c, a, b, stk);
    
82.          kar_sub(T3, c);
    
83. 
    
84.          *a = sa;
    
85.          *b = sb;
    
86. 
    
87.          /* finally, add T3 * NTL_RADIX^{hsa} to c */
    
88. 
    
89.          kar_add(c, T3, hsa);
    
90.       }
    
91.       else {
    
92.          /* degenerate case */
    
93. 
    
94.          long *T;
    
95.          
    
96.          T = stk;  stk += sb + hsa + 1;
    
97. 
    
98.          if (stk-kmem > max_kmem) zhalt("internal error: kmem overflow");
    
99. 
    
100.          /* recursively compute b*a_hi into high part of c */
     
101.          {
     
102.             long olda;
     
103. 
     
104.             olda = a[hsa];  a[hsa] = sa-hsa;
     
105.             kar_mul(c + hsa, a + hsa, b, stk);
     
106.             a[hsa] = olda;
     
107.          }
     
108. 
     
109.          /* recursively compute b*a_lo into T */
     
110. 
     
111.          *a = hsa;
     
112.          kar_mul(T, a, b, stk);
     
113.          *a = sa;
     
114. 
     
115.          /* fix-up result */
     
116. 
     
117.          kar_fix(c, T, hsa);
     
118.       }
     
119.    }
     
120. 
     
121.    /* normalize result */
     
122. 
     
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 才行。
而有许多老爷机是不具备的。