高精度计算介绍
定义
高精度计算(Arbitrary-Precision Arithmetic),也被称作大整数(bignum)计算,运用了一些算法结构来支持更大整数间的运算(数字大小超过语言内建整型)。
存储
在平常的实现中,高精度数字利用字符串表示,每一个字符表示数字的一个十进制位。因此可以说,高精度数值计算实际上是一种特别的字符串处理。
读入字符串时,数字最高位在字符串首(下标小的位置)。但是习惯上,下标最小的位置存放的是数字的 最低位,即存储反转的字符串。这么做的原因在于,数字的长度可能发生变化,但我们希望同样权值位始终保持对齐(例如,希望所有的个位都在下标 [0]
,所有的十位都在下标 [1]
……);同时,加、减、乘的运算一般都从个位开始进行,这都给了「反转存储」以充分的理由。
此后我们将一直沿用这一约定。定义一个常数 LEN = 1004
表示程序所容纳的最大长度。
由此不难写出读入高精度数字的代码:
void clear(int a[]) { for (int i = 0; i < LEN; ++i) a[i] = 0; } void read(int a[]) { static char s[LEN + 1]; scanf("%s", s); clear(a); int len = strlen(s); // 如上所述,反转 for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0'; // s[i] - '0' 就是 s[i] 所表示的数码 // 也可用 ord(s[i]) - ord('0') 的方式 }
输出也按照存储的逆序输出。由于不希望输出前导零,故这里从最高位开始向下寻找第一个非零位,从此处开始输出;终止条件 i >= 1
而不是 i >= 0
是因为当整个数字等于 时仍希望输出一个字符 0
。
void print(int a[]) { int i; for (i = LEN - 1; i >= 1; --i) if (a[i] != 0) break; for (; i >= 0; --i) putchar(a[i] + '0'); putchar('\\n'); }
拼起来就是一个完整的复读机程序:
#include <cstdio> #include <cstring> static const int LEN = 1004; int a[LEN], b[LEN]; void clear(int a[]) { for (int i = 0; i < LEN; ++i) a[i] = 0; } void read(int a[]) { static char s[LEN + 1]; scanf("%s", s); clear(a); int len = strlen(s); for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0'; } void print(int a[]) { int i; for (i = LEN - 1; i >= 1; --i) if (a[i] != 0) break; for (; i >= 0; --i) putchar(a[i] + '0'); putchar('\\n'); } int main() { read(a); print(a); return 0; }
四则运算
四则运算中难度也各不相同。最简单的是高精度加减法,其次是高精度—单精度(普通的 int
)乘法和高精度—高精度乘法,最后是高精度—高精度除法。
我们将按这个顺序分别实现所有要求的功能。
加法
高精度加法,其实就是竖式加法啦。
也就是从最低位开始,将两个加数对应位置上的数码相加,并判断是否达到或超过 。如果达到,那么处理进位:将更高一位的结果上增加 ,当前位的结果减少 。
void add(int a[], int b[], int c[]) { clear(c); // 高精度实现中,一般令数组的最大长度 LEN 比可能的输入大一些 // 然后略去末尾的几次循环,这样一来可以省去不少边界情况的处理 // 因为实际输入不会超过 1000 位,故在此循环到 LEN - 1 = 1003 已经足够 for (int i = 0; i < LEN - 1; ++i) { // 将相应位上的数码相加 c[i] += a[i] + b[i]; if (c[i] >= 10) { // 进位 c[i + 1] += 1; c[i] -= 10; } } }
试着和上一部分结合,可以得到一个加法计算器:
#include <cstdio> #include <cstring> static const int LEN = 1004; int a[LEN], b[LEN], c[LEN]; void clear(int a[]) { for (int i = 0; i < LEN; ++i) a[i] = 0; } void read(int a[]) { static char s[LEN + 1]; scanf("%s", s); clear(a); int len = strlen(s); for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0'; } void print(int a[]) { int i; for (i = LEN - 1; i >= 1; --i) if (a[i] != 0) break; for (; i >= 0; --i) putchar(a[i] + '0'); putchar('\\n'); } void add(int a[], int b[], int c[]) { clear(c); for (int i = 0; i < LEN - 1; ++i) { c[i] += a[i] + b[i]; if (c[i] >= 10) { c[i + 1] += 1; c[i] -= 10; } } } int main() { read(a); read(b); add(a, b, c); print(c); return 0; }
减法
高精度减法,也就是竖式减法啦。
从个位起逐位相减,遇到负的情况则向上一位借 。整体思路与加法完全一致。
void sub(int a[], int b[], int c[]) { clear(c); for (int i = 0; i < LEN - 1; ++i) { // 逐位相减 c[i] += a[i] - b[i]; if (c[i] < 0) { // 借位 c[i + 1] -= 1; c[i] += 10; } } }
将上一个程序中的 add()
替换成 sub()
,就有了一个减法计算器。
#include <cstdio> #include <cstring> static const int LEN = 1004; int a[LEN], b[LEN], c[LEN]; void clear(int a[]) { for (int i = 0; i < LEN; ++i) a[i] = 0; } void read(int a[]) { static char s[LEN + 1]; scanf("%s", s); clear(a); int len = strlen(s); for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0'; } void print(int a[]) { int i; for (i = LEN - 1; i >= 1; --i) if (a[i] != 0) break; for (; i >= 0; --i) putchar(a[i] + '0'); putchar('\\n'); } void sub(int a[], int b[], int c[]) { clear(c); for (int i = 0; i < LEN - 1; ++i) { c[i] += a[i] - b[i]; if (c[i] < 0) { c[i + 1] -= 1; c[i] += 10; } } } int main() { read(a); read(b); sub(a, b, c); print(c); return 0; }
试一试,输入 1 2
——输出 /9999999
,怎么程序有问题呢……
事实上,上面的代码只能处理减数 大于等于被减数 的情况。处理被减数比减数小,即 时的情况很简单。
要计算 的值,因为有 ,可以调用以上代码中的 sub
函数,写法为 sub(b,a,c)
。要得到 的值,在得数前加上负号即可。
乘法
高精度—单精度
高精度乘法,也就是竖……等会儿等会儿!
先考虑一个简单的情况:乘数中的一个是普通的 int
类型。有没有简单的处理方法呢?
一个直观的思路是直接将 每一位上的数字乘以 。从数值上来说,这个方法是正确的,但它并不符合十进制表示法,因此需要将它重新整理成正常的样子。
重整的方式,也是从个位开始逐位向上处理进位。但是这里的进位可能非常大,甚至远大于 ,因为每一位被乘上之后都可能达到 的数量级。所以这里的进位不能再简单地进行 运算,而是要通过除以 的商以及余数计算。详见代码注释,也可以参考下图展示的一个计算高精度数 乘以单精度数 的过程。
当然,也是出于这个原因,这个方法需要特别关注乘数 的范围。若它和 (或相应整型的取值上界)属于同一数量级,那么需要慎用高精度—单精度乘法。
void mul_short(int a[], int b, int c[]) { clear(c); for (int i = 0; i < LEN - 1; ++i) { // 直接把 a 的第 i 位数码乘以乘数,加入结果 c[i] += a[i] * b; if (c[i] >= 10) { // 处理进位 // c[i] / 10 即除法的商数成为进位的增量值 c[i + 1] += c[i] / 10; // 而 c[i] % 10 即除法的余数成为在当前位留下的值 c[i] %= 10; } } }
高精度—高精度
如果两个乘数都是高精度,那么竖式乘法又可以大显身手了。
回想竖式乘法的每一步,实际上是计算了若干 的和。例如计算 ,计算的就是 。
于是可以将 分解为它的所有数码,其中每个数码都是单精度数,将它们分别与 相乘,再向左移动到各自的位置上相加即得答案。当然,最后也需要用与上例相同的方式处理进位。
注意这个过程与竖式乘法不尽相同,我们的算法在每一步乘的过程中并不进位,而是将所有的结果保留在对应的位置上,到最后再统一处理进位,但这不会影响结果。
void mul(int a[], int b[], int c[]) { clear(c); for (int i = 0; i < LEN - 1; ++i) { // 这里直接计算结果中的从低到高第 i 位,且一并处理了进位 // 第 i 次循环为 c[i] 加上了所有满足 p + q = i 的 a[p] 与 b[q] 的乘积之和 // 这样做的效果和直接进行上图的运算最后求和是一样的,只是更加简短的一种实现方式 for (int j = 0; j <= i; ++j) c[i] += a[j] * b[i - j]; if (c[i] >= 10) { c[i + 1] += c[i] / 10; c[i] %= 10; } } }
除法
高精度除法,也就是竖~~~~竖式长除法啦!
竖式长除法实际上可以看作一个逐次减法的过程。例如上图中商数十位的计算可以这样理解:将 减去三次 后变得小于 ,不能再减,故此位为 。
为了减少冗余运算,我们提前得到被除数的长度 与除数的长度 ,从下标 开始,从高位到低位来计算商。这和手工计算时将第一次乘法的最高位与被除数最高位对齐的做法是一样的。
参考程序实现了一个函数 greater_eq()
用于判断被除数以下标 last_dg
为最低位,是否可以再减去除数而保持非负。此后对于商的每一位,不断调用 greater_eq()
,并在成立的时候用高精度减法从余数中减去除数,也即模拟了竖式除法的过程。
// 被除数 a 以下标 last_dg 为最低位,是否可以再减去除数 b 而保持非负 // len 是除数 b 的长度,避免反复计算 inline bool greater_eq(int a[], int b[], int last_dg, int len) { // 有可能被除数剩余的部分比除数长,这个情况下最多多出 1 位,故如此判断即可 if (a[last_dg + len] != 0) return true; // 从高位到低位,逐位比较 for (int i = len - 1; i >= 0; --i) { if (a[last_dg + i] > b[i]) return true; if (a[last_dg + i] < b[i]) return false; } // 相等的情形下也是可行的 return true; } void div(int a[], int b[], int c[], int d[]) { clear(c); clear(d); int la, lb; for (la = LEN - 1; la > 0; --la) if (a[la - 1] != 0) break; for (lb = LEN - 1; lb > 0; --lb) if (b[lb - 1] != 0) break; if (lb == 0) { puts("> <"); return; } // 除数不能为零 // c 是商 // d 是被除数的剩余部分,算法结束后自然成为余数 for (int i = 0; i < la; ++i) d[i] = a[i]; for (int i = la - lb; i >= 0; --i) { // 计算商的第 i 位 while (greater_eq(d, b, i, lb)) { // 若可以减,则减 // 这一段是一个高精度减法 for (int j = 0; j < lb; ++j) { d[i + j] -= b[j]; if (d[i + j] < 0) { d[i + j + 1] -= 1; d[i + j] += 10; } } // 使商的这一位增加 1 c[i] += 1; // 返回循环开头,重新检查 } } }