任意精度整数计算的实现 (Python 1.5 源码)

在学习 Python 源码中, 我们遇到了 longobject, 它表示任意精度的整数. 要注意的是这里 long 不是 C/Java
等语言中 long 的语义, 一般它们的 long 型是 32-bits (或有 64-bits long long). 然后在 Python 中,
long 可以是任意长度的. 实验一下: 计算 x=2<<64 (这是 unsigned long long 型最大值), 再给它 + 1,
自乘等, 打印出来你会看到很长的数字(还可以更长).

在加密解密领域, 或数论计算等领域, 通常会遇到大的整数的计算, 数字的长度达到十进制 10000 位似乎是
很可能的. 大致 RSA 加密就是利用大的质数的某些数学原理开发的, 如果随着大的数字的计算技术不断发展,
可能破解 1024 位 RSA 会成为可能.

大整数一般用较小的基数(BASE) 的数字的数组表示, 一般的大整数用如下结构表示:
int sign; // 数字的符号, 0 表示正数, 1 表示负数.
int size; // 有多少个有效数字.
int alloc; // 为下面数组实际分配了多少空间.
digit digits[]; // 具有 alloc 个空间的数组, 合起来表示一个大的数字.

这里 digit 是能够表示基数(基底, BASE) 范围的类型, 如在 Python 中 BASE 选择为 2^15, digit 类型
即选择为 unsigned short. (别的语言可能选用的 BASE=2^30, 用 uint 做 digit 类型等).

数组 digits[] 从低位到高位(共 size 位)表示一个数字, 用容易理解的 BASE=10 (即十进制)举例,
数字 1234 在数组 digits[] 中即写作 [4, 3, 2, 1], 即低位 4 在数组的索引 0 位置, 高位 1 在数组
的索引为 3 的位置.

加法和减法

对于加法, 相当于将两个大整数的 digits[] 从低到高逐个位加起来, 如果超过了 BASE 则进位. 这和列
竖式做加法非常相似. 减法也类似, 这些在以前就看过一点, 下面用 python 代码为例列出伪代码:

// longobject.c -- 做无符号整数 a,b 的加法 (不考虑符号 sign 部分).
longobject *x_add(longobject *a, longobject *b) {
int size_a = a->size,   // 数字 a 的位数
size_b = b->size;   // 数字 b 的位数.
// 假设已经保证了 size_a >= size_b. 即数字 a 的位数比 b 要多, 或相同.
// 如果不是, 则交换 a,b 的顺序即可.

// 为结果 z 分配空间. 两数字 a+b 的结果位数可能进位, 所以位数+1.
longobject *z = new longobject(size_a+1);
int carry = 0;   // 表示进位.
// 从低位到高位循环.
for (i = 0; i < size_b; ++i) {
// 将 a,b 的第 i 位加起来, 同时加上前一位的进位.
carry = a->digits[i] + b->digits[i] + carry;
z->digits[i] = carry & MASK;  // 结果, 先取低 digit 部分.
carry >>= SHIFT;   // 新的进位.
}

// 剩下 size_b 到 size_a 部分循环.
for (i = size_b; i < size_a; ++i) {
// ... 类似上面的过程.
}

// 最后的可能进位.
z->digits[i] = carry;
return 最后的规范化处理一下(z);
}

由于基底一般选择为 2^SHIFT, 取和的低部分就变成的了 & MASK 运算; 取和的进位变成右移 SHIFT 位计算.
而 MASK = 2^SHIFT - 1 = BASE-1. (在 Python 中 SHIFT=16, BASE=0x8000, MASK=0x7FFF).

减法类似, 处理一下符号, 以及大小关系. 如果 a-b a>b 则结果符号为正, 否则取 b-a, 结果取负即可.

乘法

列竖式做乘法是小学就学过的, 大整数做乘法的最一般算法也是这种方式. 设有a*b, 则从 b 的低位到高位循环,
每位乘上 a 的从低到高循环, 然后累加到结果中. 简单的伪代码如下:

// longobject.c -- 最简单的竖式乘法.
longobject *long_mul(longobject *a, longobject *b) {
int size_a = a->size,
size_b = b->size;
// 分配结果所需位数, 乘法结果长度一般是两个 size 的和(更细致一点是 <=).
longobject *z = new longobject(size_a+size_b);

// 循环 a 的各位.
for (i = 0; i < size_a; ++i) {
carry = 0;   // 当前计算位置(i+j) 的累加进位.
for (j = 0; j < size_b; ++j) {
// a[i]*b[j] 累加到 z 的第 i+j 位上. (可用10进制辅助思考理解)
carry += z[i+j] + a[i]*b[j];
z[i+j] = carry & MASK;  // 保留此位部分.
carry >>= carry;   // 剩下高部分 进位到高位. 下一个循环 j 即计算高位.
}
// 如果还剩下进位, 要不断累加到 z[] 的高位去.
for (; carry != 0; ++j) {
... 细节略, 类似于上面...
}
}
... 处理符号, 规范化, 返回结果 ...
}

这里的核心是 carry += z[i+j] + a[i]*b[i]; --- 相当于竖式乘法中每位相乘, 累加在结果上, 并进位.

假设数 a,b 都是 n 位的, 则上面的算法重点乘法 a[i]*b[i] 总计执行 n*n 次, 也即这个算法是 O(n^2) 的.
几个世纪以来通用的乘法(列竖式)都是这么做的, 如果没有计算机和超大整数(加密,数论中10000位数字似乎
很正常)计算需求, 似乎没有动力去寻求更快速的算法. 如果是我, 可能能计算乘法就很满意了, 没想过,
也想不到, 最最基本的, 从小学就学过的 --- 乘法, 居然也能找到快速计算的方法.

对于年轻的计算机学科, "近期" 在 1962年发现了 Karatsuba 算法, 这个算法通过将大的 a,b 分解为较小
四个数字的组合, 减少了 3/4 的乘法需要量.

a = a1*BASEN + a2
b = b1*BASEN + b2

这里 BASEN 是 BASE^N (基底的 N 次方). 举例数字 123456, 分解为 123*10^4+456, 就好像将数字
劈成两部分 123 和 456. 这样分解后计算 a*b 变成:

a*b = a1*b1*BASEN^2 + (a1*b2+a2*b1)*BASEN + a2*b2
= W2 * BASEN^2 + W1 * BASEN + W0

其中 W1=a1*b2+a2*b1 = (a1+a2)*(b1+b2) - a1*b1 - a2*b2
= (a1+a2)*(b1+b2) - W2 - W0

现在为计算出 W0,W1,W2 只需3次 n/2*n/2 大小数字的乘法, 即原来的 3/4 比例的乘法需要量了.

递归这样的过程, 对 a1*b1 等也使用这样的分解等, (天才的)计算机科学家们证明了最终计算量大约为
O(n^1.584). 为感性认识一下这个数字有什么变化, 假设有两个 10000 位的数字做乘法, 则原算法需:
n*n = 10000,0000 次基本乘法
Karatsuba 算法大概需要:
n^1.584 = 216,7704 次基本乘法. 只有原来的 2.2%, 多么惊人的数字!
设想黑客们, 能够只用你以为的 2.2% 的时间就破解了加密, 在你还在向老板保证加密是牢不可破的时候...

Python 1.0.1 直到 1.5 还未引入 karastuba 算法, 但在 2.5 的 longobject.c 中看到了此算法.
算法本身不是很复杂, 也有注解, 我想看看注解基本就差不多了:

/* The plan: (for Karatsuba mul)
* 1. Allocate result space (asize + bsize digits:  that's always
*    enough).
*    分配结果所需内存.
* 2. Compute ah*bh, and copy into result at 2*shift.
*    计算 ah*bh (即W0)
* 3. Compute al*bl, and copy into result at 0.  Note that this
*    can't overlap with #2.
*    计算 al*bl (即W2)
* 4. Subtract al*bl from the result, starting at shift.  This may
*    underflow (borrow out of the high digit), but we don't care:
*    we're effectively doing unsigned arithmetic mod
*    BASE**(sizea + sizeb), and so long as the *final* result fits,
*    borrows and carries out of the high digit can be ignored.
*    从结果中减去 W0*BASEN (从 shift 开始即等价于 *BASEN)
* 5. Subtract ah*bh from the result, starting at shift.
*    从结果中减去 W2*BASEN (从 shift 开始等价于 *BASEN)
* 6. Compute (ah+al)*(bh+bl), and add it into the result starting
*    at shift.
*    计算 W1, 加到 shift 开始 (等价于 *BASEN)
*/

一般两个数字 a,b 的位数大到一定程度才适合用 karastuba 算法, 因为有一些附加分解,调用等消耗.
具体 python 选择 KARATSUBA_CUTOFF=70, 即超过 70 个digit 表示的大整数才使用.
到底这个数字多少合适, 需要实际测试, 不同硬件/指令集是有区别的. 例如 intel 芯片
乘法指令:加法指令所需周期为 29:1 (某个地方写的, 芯片在不断发展, 此数字也不断在变化, 甚至不对)
则这个 CUTOFF 就小一点合适. 如果芯片这个比例是 10:1, 15:1, 或别的什么, 那 CUTOFF 就有不同...

对乘法的优化还没完, 在 karastuba 算法之后, Toom-Cook 算法将数字分解为 3 部分, 通过求解
五个 W0,W1,W2,W3,W4 构成的类似系数, 乘法的计算量降到 O(n^1.465).

更进一步, 分解为4,5,...10 部分, n 的指数继续下降到 1.403, 1.365, 1.279. 当然分解本身会带来
额外开销, 使得应用这些复杂算法只对超大整数才有意义. 具体去实现就是"疯子"的任务了, 例如请
参考 GNU MIP 库里面实现的 Toom-Cook 算法等.

这样, 我们只需要知道 python 支持大整数的乘法, 并且还有点优化就够了.

除法

在对大数字的四则运算中, 除法我是很没有弄清楚的, 我想不仔细调试程序, 以及纸上推理(最好
有专家给讲讲... ), 我是弄不太清楚的. 最多将原理大致列出来:

(具体原理请参见 Knuth 计算机程序设计艺术第二卷)

// longobject.c -- 计算大整数的除法.
// 以下我们直接写源码很难懂(我也没彻底弄懂), 所以用伪代码...
longobject *x_divrem(longobject *u1, longobject *v1, longobject **prem) {
// 假设一些简单情况都已经忽略/已处理了, 例如 v1=0, m < n 等.
m = u1->size, n = v1->size;
// D1. [规格化], 此步骤使得 v1 (除数)的高位 v1[n-1] 满足 >= BASE/2.
// 这是 Kunth 算法的要求...
d = BASE/(v1[n-1] + 1)
longobject *u = u1 * d;   // 被除数(numerator)
longobject *v = v1 * d;   // 除数(denumerator)
m = u->size;   // u1 在乘上因子 d 之后可能增加了一位. n 位数不会变.
longobjec *q = new longobject();   // 结果商 (quantity)

// D2. [初始化 j] 准备循环
for (j = m; j >= 0; --j) {
// D3. 计算 q' (试商), 余数 r'
q' = (u[j+n]*BASE + u[j+n-1]) / v[n-1];
r' = (u[j+n]*BASE + u[j+n-1]) % v[n-1];
// 测试 q' = BASE 或 q'*v[n-2] > BASE*r' + u[j+n-2]
// 如果是, 则 q'-1, r'+v[n-1] ...
while (q'* == BASE || q'*v[n-2] > ...)
--q', r' += v[n-1];

// D4. [乘和减] 此步骤计算 u = u - q*v, 可能结果 u < 0, 要保留借位.
u = u - q*v;  // 此步骤同时计算一个乘法和一个减法.
// 实际 python 程序用一个 for 循环完成, 借位为 carry.

// D5. [测试余数] q[j] = q. 如果 D4 有借位,则转到 D6, 否则循环
q[j] = q;
if (carry == 0) 没有借位则 continue 循环.

// D6. [往回加] 由于很小的概率使得 u = u-q*v < 0, 需要 u+q
u = u+q; q[j] = q[j] - 1;

// D7. 即继续循环.
}

// D8. 上述循环完成后, 求出商为 q, 余数在 u 中(但乘上了因子 d).
// 所以 [逆规格化] 得到余数 r(remainder)
r = u / d;
return <q, r>;  // 返回 q,r 对表示 商/余数.
}

神奇的是 Knuth 在他的著作中既然还算出了 D6 步骤的可能概率: 2/BASE.

仅仅为了整数的加减乘除, 计算机科学家/程序员们就已经付出了如此巨大的努力, 为此就绝不能轻视这
一工作的价值, 因为程序员们并不仅仅是简单地作为"翻译"存在. 那些不能够理解计算机器的精微奥秘
之处的, 而试图用它进行计算时, 一旦超出了所用语言的能力时就变得无能为力. 只有计算机专家才能够
深刻地理解为什么, 通过努力研究并克服这一障碍, 使得其它专业的使用者才能够实现他们的计算需求.

知识无涯, 而吾生则有涯, 没有更多的时间去研究/学习这无涯的知识, 只能遗憾地止步于此了.