两个n位的数相加，结果可能出现进位，这时可用一个变量carry 表示进位(carry=0 或者 carry=1)，一个长为n的数组表示结果。
两个n位的数相减，结果可能出现借位，如果出现借位,carry=-1， 结果采用补码表示，数组的各个元素使用补码表示。 下面贴出GMP kara乘法 的代码：来自GMP 4.22 中的 mpn\generic\mul_n.c ，GMP 4.22 可从 http://swox.com/gmp/ 下载
顺便说一下，我的代码比他的简单，而且是基于2^30进制的。/* Multiplies using 3 half-sized mults and so on recursively.
* p := product of a and b.
* No overlap of p[...] with a[...] or b[...].
* ws is workspace.
*/
void
mpn_kara_mul_n (mp_ptr p, mp_srcptr a, mp_srcptr b, mp_size_t n, mp_ptr ws)
{
mp_limb_t w, w0, w1;
mp_size_t n2;
mp_srcptr x, y;
mp_size_t i;
int sign;
n2 = n >> 1;
ASSERT (n2 > 0);
if ((n & 1) != 0)
{
/* Odd length. */
mp_size_t n1, n3, nm1;
n3 = n - n2;
sign = 0;
w = a;
if (w != 0)
w -= mpn_sub_n (p, a, a + n3, n2);
else
{
i = n2;
do
{
--i;
w0 = a[ i ];
w1 = a;
}
while (w0 == w1 && i != 0);
if (w0 < w1)
{
x = a + n3;
y = a;
sign = ~0;
}
else
{
x = a;
y = a + n3;
}
mpn_sub_n (p, x, y, n2);
}
p = w;
w = b;
if (w != 0)
w -= mpn_sub_n (p + n3, b, b + n3, n2);
else
{
i = n2;
do
{
--i;
w0 = b[ i ];
w1 = b;
}
while (w0 == w1 && i != 0);
if (w0 < w1)
{
x = b + n3;
y = b;
sign = ~sign;
}
else
{
x = b;
y = b + n3;
}
mpn_sub_n (p + n3, x, y, n2);
}
p = w;
n1 = n + 1;
if (n2 < MUL_KARATSUBA_THRESHOLD)
{
if (n3 < MUL_KARATSUBA_THRESHOLD)
{
mpn_mul_basecase (ws, p, n3, p + n3, n3);
mpn_mul_basecase (p, a, n3, b, n3);
}
else
{
mpn_kara_mul_n (ws, p, p + n3, n3, ws + n1);
mpn_kara_mul_n (p, a, b, n3, ws + n1);
}
mpn_mul_basecase (p + n1, a + n3, n2, b + n3, n2);
}
else
{
mpn_kara_mul_n (ws, p, p + n3, n3, ws + n1);
mpn_kara_mul_n (p, a, b, n3, ws + n1);
mpn_kara_mul_n (p + n1, a + n3, b + n3, n2, ws + n1);
}
if (sign)
mpn_add_n (ws, p, ws, n1);
else
mpn_sub_n (ws, p, ws, n1);
nm1 = n - 1;
if (mpn_add_n (ws, p + n1, ws, nm1))
{
mp_limb_t x = (ws + 1) & GMP_NUMB_MASK;
ws = x;
if (x == 0)
ws = (ws + 1) & GMP_NUMB_MASK;
}
if (mpn_add_n (p + n3, p + n3, ws, n1))
{
mpn_incr_u (p + n1 + n3, 1);
}
}
else
{
/* Even length. */
i = n2;
do
{
--i;
w0 = a[ i ];
w1 = a;
}
while (w0 == w1 && i != 0);
sign = 0;
if (w0 < w1)
{
x = a + n2;
y = a;
sign = ~0;
}
else
{
x = a;
y = a + n2;
}
mpn_sub_n (p, x, y, n2);
i = n2;
do
{
--i;
w0 = b[ i ];
w1 = b;
}
while (w0 == w1 && i != 0);
if (w0 < w1)
{
x = b + n2;
y = b;
sign = ~sign;
}
else
{
x = b;
y = b + n2;
}
mpn_sub_n (p + n2, x, y, n2);
/* Pointwise products. */
if (n2 < MUL_KARATSUBA_THRESHOLD)
{
mpn_mul_basecase (ws, p, n2, p + n2, n2);
mpn_mul_basecase (p, a, n2, b, n2);
mpn_mul_basecase (p + n, a + n2, n2, b + n2, n2);
}
else
{
mpn_kara_mul_n (ws, p, p + n2, n2, ws + n);
mpn_kara_mul_n (p, a, b, n2, ws + n);
mpn_kara_mul_n (p + n, a + n2, b + n2, n2, ws + n);
}
/* Interpolate. */
if (sign)
w = mpn_add_n (ws, p, ws, n);
else
w = -mpn_sub_n (ws, p, ws, n);
w += mpn_add_n (ws, p + n, ws, n);
w += mpn_add_n (p + n2, p + n2, ws, n);
MPN_INCR_U (p + n2 + n, 2 * n - (n2 + n), w);
}
} 下面贴出NTL 中的kara乘法 源代码，kar_mul函数来自 WinNTL-5_4_2\src\c_lip_impl.h，WinNTL可从 http://www.shoup.net下载static
void kar_mul(long *c, long *a, long *b, long *stk)
{
long sa, sb, sc;
if (*a < *b) { long *t = a; a = b; b = t; }
sa = *a;
sb = *b;
sc = sa + sb;
if (sb < KARX) {
/* classic algorithm */
long *pc, i, *pb;
pc = c;
for (i = sc; i; i--) {
pc++;
*pc = 0;
}
pc = c;
pb = b;
for (i = sb; i; i--) {
pb++;
pc++;
zaddmul(*pb, pc, a);
}
}
else {
long hsa = (sa + 1) >> 1;
if (hsa < sb) {
/* normal case */
long *T1, *T2, *T3;
/* allocate space */
T1 = stk;stk += hsa + 2;
T2 = stk;stk += hsa + 2;
T3 = stk;stk += (hsa << 1) + 3;
if (stk-kmem > max_kmem) zhalt("internal error: kmem overflow");
/* compute T1 = a_lo + a_hi */
kar_fold(T1, a, hsa);
/* compute T2 = b_lo + b_hi */
kar_fold(T2, b, hsa);
/* recursively compute T3 = T1 * T2 */
kar_mul(T3, T1, T2, stk);
/* recursively compute a_hi * b_hi into high part of c */
/* and subtract from T3 */
{
long olda, oldb;
olda = a;a = sa-hsa;
oldb = b;b = sb-hsa;
kar_mul(c + (hsa << 1), a + hsa, b + hsa, stk);
kar_sub(T3, c + (hsa << 1));
a = olda;
b = oldb;
}
/* recursively compute a_lo*b_lo into low part of c */
/* and subtract from T3 */
*a = hsa;
*b = hsa;
kar_mul(c, a, b, stk);
kar_sub(T3, c);
*a = sa;
*b = sb;
/* finally, add T3 * NTL_RADIX^{hsa} to c */
kar_add(c, T3, hsa);
}
else {
/* degenerate case */
long *T;
T = stk;stk += sb + hsa + 1;
if (stk-kmem > max_kmem) zhalt("internal error: kmem overflow");
/* recursively compute b*a_hi into high part of c */
{
long olda;
olda = a;a = sa-hsa;
kar_mul(c + hsa, a + hsa, b, stk);
a = olda;
}
/* recursively compute b*a_lo into T */
*a = hsa;
kar_mul(T, a, b, stk);
*a = sa;
/* fix-up result */
kar_fix(c, T, hsa);
}
}
/* normalize result */
while (c == 0 && sc > 1) sc--;
*c = sc;
} :)
两个采取了不同的策略啊 原帖由 无心人 于 2008-4-28 09:12 发表 http://images.5d6d.net/dz60/common/back.gif
是否存在比较的必要？
考虑工作量
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 ...
其实，前加法也并非一定会出现进位。 应该说，直接采用楼主的公式还不是最优化的，
更好的形式可以进一步减少一次 大数的加/减法 运算。
下面的公式是从我的 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% 左右。 我公布一下我的kara乘法的 技术数据：
1.因为我的加减法接口是这样的
intDEC_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乘法中，除了乘法外的其他运算的次数。
以第1种形式分析
// ( A1*b + A0 )( B1*b + B0 )// = (( A1*B1*b - A0*B0 ) + ( A0 + A1 )( B0 + B1 ))*b
// - ( A1*B1*b - A0*B0 )
2*len 次整数减法 - A0*B0
3*len 次整数复制 ( A1*B1*b - A0*B0 )
len 次整数加法 + A1
len 次整数加法 + B1
2*len 次整数加法 + ( A0 + A1 )( B0 + B1 )
3*len 次整数减法 - ( A1*B1*b - A0*B0 )
注：选择好适当大小的 b，可以确保上述所有减法无须修正（即结果一定为正）。
合计：加减法 9*len 单元次；复制 3*len 单元次；乘法 3*len*len 单元次；移位 无（在相关运算中直接定位即可） GxQ
你没考虑过你的加减法和求负还有内存拷贝都不是最优的
:)
加减存在一次处理4X32位操作的算法
求负，拷贝存在一次处理4X128位操作的算法
:) 那需要用户CPU至少得支持 SSE2 才行。:(
而有许多老爷机是不具备的。