来学学数学啦。

数学基础

计算机科学与数学紧密相关，而在算法竞赛中尤其强调以数论、排列组合、概率期望、多项式为代表离散、具体的数学：其注重程序实现和现实问题，可以出现在几乎任何类别的题目中。

算法竞赛中的常见符号

以下摘自wiki

整除/同余理论常见符号

  1. 整除符号：\(x\mid y\)，表示 \(x\) 整除 \(y\)，即 \(x\) 是 \(y\) 的因数。
  2. 取模符号：\(x\bmod y\)，表示 \(x\) 除以 \(y\) 得到的余数。
  3. 互质符号：\(x\perp y\)，表示 \(x\)，\(y\) 互质。
  4. 最大公约数：\(\gcd(x,y)\)，在无混淆意义的时侯可以写作 \((x,y)\)。
  5. 最小公倍数：\(\operatorname{lcm}(x,y)\)，在无混淆意义的时侯可以写作 \([x,y]\)。

数论函数常见符号

求和符号：\(\sum\) 符号，表示满足特定条件的数的和。举几个例子：

  • \(\sum_{i=1}^n i\) 表示 \(1+2+\dotsb+n\) 的和。其中 \(i\) 是一个变量，在求和符号的意义下 \(i\) 通常是 正整数或者非负整数（除非特殊说明）。这个式子的含义可以理解为，\(i\) 从 \(1\) 循环到 \(n\)，所有 \(i\) 的和。这个式子用代码的形式很容易表达。当然，学过简单的组合数学的同学都知道 \(\sum_{i=1}^n i=\dfrac{n(n+1)}{2}\)。
  • \(\sum_{S\subseteq T}|S|\) 表示所有被 \(T\) 包含的集合的大小的和。
  • \(\sum_{p\le n,p\perp n}1\) 表示的是 \(n\) 以内有多少个与 \(n\) 互质的数，即 \(\varphi(n)\)，\(\varphi\) 是欧拉函数。

求积符号：\(\prod\) 符号，表示满足特定条件的数的积。举几个例子：

  • \(\prod_{i=1}^ni\) 表示 \(n\) 的阶乘，即 \(n!\)。在组合数学常见符号中会讲到。
  • \(\prod_{i=1}^na_i\) 表示 \(a_1\times a_2\times a_3\times \dotsb\times a_n\)。
  • \(\prod_{x|d}x\) 表示 \(d\) 的所有因数的乘积。

在行间公式中，求和符号与求积符号的上下条件会放到符号的上面和下面，这一点要注意。

其他常见符号

  1. 阶乘符号 \(!\)，\(n!\) 表示 \(1\times 2\times 3\times \dotsb \times n\)。特别地，\(0!=1\)。
  2. 向下取整符号：\(\lfloor x\rfloor\)，表示小于等于 \(x\) 的最大的整数。常用于分数，比如分数的向下取整 \(\left\lfloor\dfrac{x}{y}\right\rfloor\)。
  3. 向上取整符号：\(\lceil x\rceil\)，与向下取整符号相对，表示大于等于 \(x\) 的最小的整数。

位运算基础

位运算就是基于整数的二进制表示进行的运算。由于计算机内部就是以二进制来存储数据，位运算是相当快的。

基本的位运算共 \(6\) 种，分别为按位与、按位或、按位异或、按位取反、左移和右移。

与、或、异或

这三者都是两数间的运算，因此在这里一起讲解。

它们都是将两个整数作为二进制数，对二进制表示中的每一位逐一运算。

运算 运算符 数学符号表示 解释
与 & \(\&\)、\(\operatorname{and}\) 只有两个对应位都为 \(1\) 时才为 \(1\)
或 | \(\mid\)、\(\operatorname{or}\) 只要两个对应位中有一个 \(1\) 时就为 \(1\)
异或 ^ \(\oplus\)、\(\operatorname{xor}\) 只有两个对应位不同时才为 \(1\)

注意区分逻辑与（对应的数学符号为 \(\wedge\)）和按位与、逻辑或（\(\vee\)）和按位或的区别。网络中的资料中使用的符号多有不规范之处，以上下文为准。

异或运算的逆运算是它本身，也就是说两次异或同一个数最后结果不变，即 \(a \oplus b \oplus b = a\) 。

取反

取反是对一个数 \(num\) 进行的位运算，即单目运算。

取反暂无默认的数学符号表示，其对应的运算符为 ~。它的作用是把 \(num\) 的二进制补码中的 \(0\) 和 \(1\) 全部取反（\(0\) 变为 \(1\)，\(1\) 变为 \(0\)）。有符号整数的符号位在 ~ 运算中同样会取反。

补码：在二进制表示下，正数和 \(0\) 的补码为其本身，负数的补码是将其对应正数按位取反后加一

左移和右移

num << i 表示将 \(num\) 的二进制表示向左移动 \(i\) 位所得的值。

num >> i 表示将 \(num\) 的二进制表示向右移动 \(i\) 位所得的值。

移位运算中如果出现如下情况，则其行为未定义：

  1. 右操作数（即移位数）为负值；
  2. 右操作数大于等于左操作数的位数；

例如，对于 int 类型的变量 a ， a<<-1 和 a<<32 都是未定义的。

对于左移操作，需要确保移位后的结果能被原数的类型容纳，否则行为也是未定义的。[^note1]对一个负数执行左移操作也未定义。[^note2]

对于右移操作，右侧多余的位将会被舍弃，而左侧较为复杂：对于无符号数，会在左侧补 \(0\)；而对于有符号数，则会用最高位的数（其实就是符号位，非负数为 \(0\)，负数为 \(1\)）补齐。[^note3]

复合赋值位运算符

和 += , -= 等运算符类似，位运算也有复合赋值运算符： &= , |= , ^= , <<= , >>= 。（取反是单目运算，所以没有。）

关于优先级

位运算的优先级低于算术运算符（除了取反），而按位与、按位或及异或低于比较运算符，所以使用时需多加注意，在必要时添加括号。

比较需要注意的是，取反(~)和逻辑非(!)的优先级特别高，比如 !a=1 ! 会优先结合 a 而不会先运算 a=1。其它位运算的优先级又特别低，比如 l+r>>1 可以表示为 (l+r)/2。

位运算的应用

位运算一般有三种作用：

  1. 高效地进行某些运算，代替其它低效的方式。

  2. 表示集合。（常用于状压 DP 。）

  3. 题目本来就要求进行位运算。

需要注意的是，用位运算代替其它运算方式（即第一种应用）在很多时候并不能带来太大的优化，反而会使代码变得复杂，使用时需要斟酌。（但像“乘 2 的非负整数次幂”和“除以 2 的非负整数次幂”就最好使用位运算，因为此时使用位运算可以优化复杂度。）

有关 2 的幂的应用

由于位运算针对的是变量的二进制位，因此可以推广出许多与 2 的整数次幂有关的应用。

将一个数乘（除） 2 的非负整数次幂：

1
2
3
4
5
6
7
// C++ Version
int mulPowerOfTwo(int n, int m) {  // 计算 n*(2^m)
  return n << m;
}
int divPowerOfTwo(int n, int m) {  // 计算 n/(2^m)
  return n >> m;
}

我们平常写的除法是向 \(0\) 取整，而这里的右移是向下取整（注意这里的区别），即当数大于等于 \(0\) 时两种方法等价，当数小于 \(0\) 时会有区别，如： -1 / 2 的值为 \(0\) ，而 -1 >> 1 的值为 \(-1\) 。

判断一个数是不是 \(2\) 的非负整数次幂：

1
2
// C++ Version
bool isPowerOfTwo(int n) { return n > 0 && (n & (n - 1)) == 0; }

对 \(2\) 的非负整数次幂取模：

1
2
// C++ Version
int modPowerOfTwo(int x, int mod) { return x & (mod - 1); }

模拟集合操作

一个数的二进制表示可以看作是一个集合（\(0\) 表示不在集合中，\(1\) 表示在集合中）。比如集合 \(\{1,3,4,8\}\) ，可以表示成 \((100011010)_2\) 。而对应的位运算也就可以看作是对集合进行的操作。

操作 集合表示 位运算语句
交集 \(a \cap b\) a & b
并集 \(a \cup b\) a|b
补集 \(\bar{a}\) ~a （全集为二进制都是 1）
差集 \(a \setminus b\) a & (~b)
对称差 \(a\triangle b\) a ^ b

子集遍历：

1
2
3
4
// 遍历 u 的非空子集
for (int s = u; s; s = (s - 1) & u) {
  // s 是 u 的一个非空子集
}

用这种方法可以在 \(O(2^{\text{popcount}(u)})\) （ \(\text{popcount}(u)\) 表示 \(u\) 二进制中 1 的个数）的时间复杂度内遍历 \(u\) 的子集，进而可以在 \(O(3^n)\) 的时间复杂度内遍历大小为 \(n\) 的集合的每个子集的子集。（复杂度为 \(O(3^n)\) 是因为每个元素都有 不在大子集中/只在大子集中/同时在大小子集中 三种状态。）

汉明权重

汉明权重是一串符号中不同于（定义在其所使用的字符集上的）零符号（zero-symbol）的个数。对于一个二进制数，它的汉明权重就等于它 \(1\) 的个数（即 popcount）。

求一个数的汉明权重可以循环求解：我们不断地去掉这个数在二进制下的最后一位（即右移 \(1\) 位），维护一个答案变量，在除的过程中根据最低位是否为 \(1\) 更新答案。

代码如下：

1
2
3
4
5
6
7
8
9
// 求 x 的汉明权重
int popcount(int x) {
    int cnt = 0;
    while (x) {
        cnt += x & 1;
        x >>= 1;
    }
    return cnt;
}

求一个数的汉明权重还可以使用 lowbit 操作：我们将这个数不断地减去它的 lowbit[^note4]，直到这个数变为 \(0\)。

代码如下：

1
2
3
4
5
6
7
8
9
// 求 x 的汉明权重
int popcount(int x) {
    int cnt = 0;
    while (x) {
        cnt++; 
        x -= x & -x;
    }
    return cnt;
}

构造汉明权重递增的排列

在 状压 DP 中，按照 popcount 递增的顺序枚举有时可以避免重复枚举状态。这是构造汉明权重递增的排列的一大作用。

下面我们来具体探究如何在 \(O(n)\) 时间内构造汉明权重递增的排列。

我们知道，一个汉明权重为 \(n\) 的最小的整数为 \(2^n-1\)。只要可以在常数时间构造出一个整数汉明权重相等的后继，我们就可以通过枚举汉明权重，从 \(2^n-1\) 开始不断寻找下一个数的方式，在 \(O(n)\) 时间内构造出 \(0\sim n\) 的符合要求的排列。

而找出一个数 \(x\) 汉明权重相等的后继有这样的思路，以 \((10110)_2\) 为例：

  • 把 \((10110)_2\) 最右边的 \(1\) 向左移动，如果不能移动，移动它左边的 \(1\)，以此类推，得到 \((11010)_2\)。

  • 把得到的 \((11010)_2\) 最后移动的 \(1\) 原先的位置一直到最低位的所有 \(1\) 都移到最右边。这里最后移动的 \(1\) 原来在第三位，所以最后三位 \(010\) 要变成 \(001\)，得到 \((11001)_2\)。

这个过程可以用位运算优化：

1
2
int t = x + (x & -x);
x = t | ((((t&-t)/(x&-x))>>1)-1);
  • 第一个步骤中，我们把数 \(x\) 加上它的 lowbit，在二进制表示下，就相当于把 \(x\) 最右边的连续一段 \(1\) 换成它左边的一个 \(1\)。如刚才提到的二进制数 \((10110)_2\)，它在加上它的 lowbit 后是 \((11000)_2\)。这其实得到了我们答案的前半部分。
  • 我们接下来要把答案后面的 \(1\) 补齐，\(t\) 的 lowbit 是 \(x\) 最右边连续一段 \(1\) 最左边的 \(1\) 移动后的位置，而 \(x\) 的 lowbit 则是 \(x\) 最右边连续一段 \(1\) 最左边的位置。还是以 \((10110)_2\) 为例，\(t = (11000)_2\)，\(\operatorname{lowbit}(t) = (01000)_2\)，\(\operatorname{lowbit}(x)=(00010)_2\)。
  • 接下来的除法操作是这种位运算中最难理解的部分，但也是最关键的部分。我们设原数最右边连续一段 \(1\) 最高位的 \(1\) 在第 \(r\) 位上（位数从 \(0\) 开始），最低位的 \(1\) 在第 \(l\) 位，\(t\) 的 lowbit 等于 1 << (r+1) ，\(x\) 的 lowbit 等于 1 << l， (((t&-t)/(x&-x))>>1) 得到的，就是 (1<<(r+1))/(1<<l)/2 = (1<<r)/(1<<l) = 1<<(r-l) ，在二进制表示下就是 \(1\) 后面跟上 \(r-l\) 个零，零的个数正好等于连续 \(1\) 的个数减去 \(1\) 。举我们刚才的数为例，\(\frac{\operatorname{lowbit(t)\div 2}}{\operatorname{lowbit(x)}} = \frac{(00100)_2}{(00010)_2} = (00010)_2\) 。把这个数减去 \(1\) 得到的就是我们要补全的低位，或上原来的数就可以得到答案。

所以枚举 \(0\sim n\) 按汉明权重递增的排列的完整代码为：

1
2
3
4
5
for (int i = 0; (1<<i)-1 <= n; i++) {
    for (int x = (1<<i)-1, t; x <= n; t = x+(x&-x), x = x ? (t|((((t&-t)/(x&-x))>>1)-1)) : (n+1)) {
        // 写下需要完成的操作
    }
}

其中要注意 \(0\) 的特判，因为 \(0\) 没有相同汉明权重的后继。

STL

如果需要更大的位数支持可以使用 std::bitset，它的位运算符号被重载过，可以就当一个数来操作。

快速幂（板子）

1
2
3
4
5
6
7
8
9
10
long long qpow(int base,int exp){
    long long ans=1;
    for(long long k=base;exp;exp>>=1){
        if(exp&1){
            ans=(ans*k)%MOD;
        }
        k=k*k%MOD;
    }
    return ans;
}

高精度运算（板子）

存储

在平常的实现中，高精度数字利用字符串表示，每一个字符表示数字的一个十进制位。因此可以说，高精度数值计算实际上是一种特别的字符串处理。

读入字符串时，数字最高位在字符串首（下标小的位置）。但是习惯上，下标最小的位置存放的是数字的 最低位，即存储反转的字符串。这么做的原因在于，数字的长度可能发生变化，但我们希望同样权值位始终保持对齐（例如，希望所有的个位都在下标 [0]，所有的十位都在下标 [1]……）；同时，加、减、乘的运算一般都从个位开始进行（回想小学的竖式运算～），这都给了「反转存储」以充分的理由。

此后我们将一直沿用这一约定。定义一个常数 LEN = 1004 表示程序所容纳的最大长度。

由此不难写出读入高精度数字的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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\) 时仍希望输出一个字符 0。

1
2
3
4
5
6
7
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');
}

拼起来就是一个完整的复读机程序咯。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#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）乘法和高精度—高精度乘法，最后是高精度—高精度除法。

我们将按这个顺序分别实现所有要求的功能。

加法

高精度加法，其实就是竖式加法啦。

也就是从最低位开始，将两个加数对应位置上的数码相加，并判断是否达到或超过 \(10\)。如果达到，那么处理进位：将更高一位的结果上增加 \(1\)，当前位的结果减少 \(10\)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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;
    }
  }
}

试着和上一部分结合，可以得到一个加法计算器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#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;
}

减法

高精度减法，也就是竖式减法啦。

从个位起逐位相减，遇到负的情况则向上一位借 \(1\)。整体思路与加法完全一致。

1
2
3
4
5
6
7
8
9
10
11
12
13
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()，就有了一个减法计算器。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#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，诶这个怎么给了我一份假的代码啊……

事实上，上面的代码只能处理减数 \(a\) 大于等于被减数 \(b\) 的情况。处理被减数比减数小，即 \(a<b\) 时的情况很简单。

\(a-b=-(b-a)\)

要计算 \(b-a\) 的值，因为有 \(b>a\)，可以调用以上代码中的 sub 函数，写法为 sub(b,a,c)。要得到 \(a-b\) 的值，在得数前加上负号即可。

乘法

高精度—单精度

高精度乘法，也就是竖……等会儿等会儿！

先考虑一个简单的情况：乘数中的一个是普通的 int 类型。有没有简单的处理方法呢？

一个直观的思路是直接将 \(a\) 每一位上的数字乘以 \(b\)。从数值上来说，这个方法是正确的，但它并不符合十进制表示法，因此需要将它重新整理成正常的样子。

重整的方式，也是从个位开始逐位向上处理进位。但是这里的进位可能非常大，甚至远大于 \(9\)，因为每一位被乘上之后都可能达到 \(9b\) 的数量级。所以这里的进位不能再简单地进行 \(-10\) 运算，而是要通过除以 \(10\) 的商以及余数计算。详见代码注释，也可以参考下图展示的一个计算高精度数 \(1337\) 乘以单精度数 \(42\) 的过程。

当然，也是出于这个原因，这个方法需要特别关注乘数 \(b\) 的范围。若它和 \(10^9\)（或相应整型的取值上界）属于同一数量级，那么需要慎用高精度—单精度乘法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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;
        }
    }
}
高精度—高精度

如果两个乘数都是高精度，那么竖式乘法又可以大显身手了。

回想竖式乘法的每一步，实际上是计算了若干 \(a \times b_i \times 10^i\) 的和。例如计算 \(1337 \times 42\)，计算的就是 \(1337 \times 2 \times 10^0 + 1337 \times 4 \times 10^1\)。

于是可以将 \(b\) 分解为它的所有数码，其中每个数码都是单精度数，将它们分别与 \(a\) 相乘，再向左移动到各自的位置上相加即得答案。当然，最后也需要用与上例相同的方式处理进位。

注意这个过程与竖式乘法不尽相同，我们的算法在每一步乘的过程中并不进位，而是将所有的结果保留在对应的位置上，到最后再统一处理进位，但这不会影响结果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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;
        }
    }
}

除法

高精度除法，也就是竖～～～～竖式长除法啦！

竖式长除法实际上可以看作一个逐次减法的过程。例如上图中商数十位的计算可以这样理解：将 \(45\) 减去三次 \(12\) 后变得小于 \(12\)，不能再减，故此位为 \(3\)。

为了减少冗余运算，我们提前得到被除数的长度 \(l_a\) 与除数的长度 \(l_b\)，从下标 \(l_a - l_b\) 开始，从高位到低位来计算商。这和手工计算时将第一次乘法的最高位与被除数最高位对齐的做法是一样的。

参考程序实现了一个函数 greater_eq() 用于判断被除数以下标 last_dg 为最低位，是否可以再减去除数而保持非负。此后对于商的每一位，不断调用 greater_eq()，并在成立的时候用高精度减法从余数中减去除数，也即模拟了竖式除法的过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
// 被除数 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;
            // 返回循环开头，重新检查
        }
    }
}

入门篇完成！

将上面介绍的四则运算的实现结合，即可完成开头提到的计算器程序。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
static const int LEN = 1004;

int a[LEN], b[LEN], c[LEN], d[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;
        }
    }
}

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;
        }
    }
}

void mul(int a[], int b[], int c[]) {
    clear(c);

    for (int i = 0; i < LEN - 1; ++i) {
        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;
        }
    }
}

inline bool greater_eq(int a[], int b[], int last_dg, int len) {
    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;
    }

    for (int i = 0; i < la; ++i) d[i] = a[i];
    for (int i = la - lb; i >= 0; --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;
                }
            }
            c[i] += 1;
        }
    }
}

int main() {
    read(a);

    char op[4];
    scanf("%s", op);

    read(b);

    switch (op[0]) {
        case '+':
            add(a, b, c);
            print(c);
            break;
        case '-':
            sub(a, b, c);
            print(c);
            break;
        case '*':
            mul(a, b, c);
            print(c);
            break;
        case '/':
            div(a, b, c, d);
            print(c);
            print(d);
            break;
        default:
            puts("> <");
    }

    return 0;
}

压位高精度

在一般的高精度加法，减法，乘法运算中，我们都是将参与运算的数拆分成一个个单独的数码进行运算。

例如计算 \(8192\times 42\) 时，如果按照高精度乘高精度的计算方式，我们实际上算的是 \((8000+100+90+2)\times(40+2)\)。

在位数较多的时候，拆分出的数也很多，高精度运算的效率就会下降。

有没有办法作出一些优化呢？

注意到拆分数字的方式并不影响最终的结果，因此我们可以将若干个数码进行合并。

还是以上面这个例子为例，如果我们每两位拆分一个数，我们可以拆分成 \((8100+92)\times 42\)。

这样的拆分不影响最终结果，但是因为拆分出的数字变少了，计算效率也就提升了。

从 进位制 的角度理解这一过程，我们通过在较大的进位制（上面每两位拆分一个数，可以认为是在 \(100\) 进制下进行运算）下进行运算，从而达到减少参与运算的数字的位数，提升运算效率的目的。

这就是 压位高精度 的思想。

下面我们给出压位高精度的加法代码，用于进一步阐述其实现方法：

1
2
3
4
5
6
7
8
9
10
11
12
// 这里的 a,b,c 数组均为 p 进制下的数
// 最终输出答案时需要将数字转为十进制
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] >= p) {  // 在普通高精度运算下，p=10
            c[i + 1] += 1;
            c[i] -= p;
        }
    }
}

压位高精下的高效竖式除法

在使用压位高精时，如果试商时仍然使用上文介绍的方法，由于试商次数会很多，计算常数会非常大。例如在万进制下，平均每个位需要试商 5000 次，这个巨大的常数是不可接受的。因此我们需要一个更高效的试商办法。

我们可以把 double 作为媒介。假设被除数有 4 位，是 \(a_4,a_3,a_2,a_1\)，除数有 3 位，是 \(b_3,b_2,b_1\)，那么我们只要试一位的商：使用 \(base\) 进制，用式子 \(\dfrac{a_4 base + a_3}{b_3 + b_2 base^{-1} + (b_1+1)base^{-2}}\) 来估商。而对于多个位的情况，就是一位的写法加个循环。由于除数使用 3 位的精度来参与估商，能保证估的商 q' 与实际商 q 的关系满足 \(q-1 \le q' \le q\)，这样每个位在最坏的情况下也只需要两次试商。但与此同时要求 \(base^3\) 在 double 的有效精度内，即 \(base^3 < 2^{53}\)，所以在运用这个方法时建议不要超过 32768 进制，否则很容易因精度不足产生误差从而导致错误。

另外，由于估的商总是小于等于实际商，所以还有再进一步优化的空间。绝大多数情况下每个位只估商一次，这样在下一个位估商时，虽然得到的商有可能因为前一位的误差造成试商结果大于等于 base，但这没有关系，只要在最后再最后做统一进位便可。举个例子，假设 base 是 10，求 \(395081/9876\)，试商计算步骤如下：

  1. 首先试商计算得到 \(3950/988=3\)，于是 \(395081-(9876 \times 3 \times 10^1) = 98801\)，这一步出现了误差，但不用管，继续下一步计算。
  2. 对余数 98801 继续试商计算得到 \(9880/988=10\)，于是 \(98801-(9876 \times 10 \times 10^0) = 41\)，这就是最终余数。
  3. 把试商过程的结果加起来并处理进位，即 \(3 \times 10^1 + 10 \times 10^0 = 40\) 便是准确的商。

方法虽然看着简单，但具体实现上很容易进坑，所以以下提供一个经过多番验证确认没有问题的实现供大家参考，要注意的细节也写在注释当中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
// 完整模板和实现 https://baobaobear.github.io/post/20210228-bigint1/
    // 对b乘以mul再左移offset的结果相减，为除法服务
BigIntSimple &sub_mul(const BigIntSimple &b, int mul, int offset) {
    if (mul == 0) return *this;
    int borrow = 0;
    // 与减法不同的是，borrow可能很大，不能使用减法的写法
    for (size_t i = 0; i < b.v.size(); ++i) {
        borrow += v[i + offset] - b.v[i] * mul - BIGINT_BASE + 1;
        v[i + offset] = borrow % BIGINT_BASE + BIGINT_BASE - 1;
        borrow /= BIGINT_BASE;
    }
    // 如果还有借位就继续处理
    for (size_t i = b.v.size(); borrow; ++i) {
        borrow += v[i + offset] - BIGINT_BASE + 1;
        v[i + offset] = borrow % BIGINT_BASE + BIGINT_BASE - 1;
        borrow /= BIGINT_BASE;
    }
    return *this;
}
BigIntSimple div_mod(const BigIntSimple &b, BigIntSimple &r) const {
    BigIntSimple d;
    r = *this;
    if (absless(b)) return d;
    d.v.resize(v.size() - b.v.size() + 1);
    // 提前算好除数的最高三位+1的倒数，若最高三位是a3,a2,a1
    // 那么db是a3+a2/base+(a1+1)/base^2的倒数，最后用乘法估商的每一位
    // 此法在BIGINT_BASE<=32768时可在int32范围内用
    // 但即使使用int64，那么也只有BIGINT_BASE<=131072时可用（受double的精度限制）
    // 能保证估计结果q'与实际结果q的关系满足q'<=q<=q'+1
    // 所以每一位的试商平均只需要一次，只要后面再统一处理进位即可
    // 如果要使用更大的base，那么需要更换其它试商方案
    double t = (b.get((unsigned)b.v.size() - 2) +
                (b.get((unsigned)b.v.size() - 3) + 1.0) / BIGINT_BASE);
    double db = 1.0 / (b.v.back() + t / BIGINT_BASE);
    for (size_t i = v.size() - 1, j = d.v.size() - 1; j <= v.size();) {
        int rm = r.get(i + 1) * BIGINT_BASE + r.get(i);
        int m = std::max((int)(db * rm), r.get(i + 1));
        r.sub_mul(b, m, j);
        d.v[j] += m;
        if (!r.get(i + 1))  // 检查最高位是否已为0，避免极端情况
            --i, --j;
    }
    r.trim();
    // 修正结果的个位
    int carry = 0;
    while (!r.absless(b)) {
        r.subtract(b);
        ++carry;
    }
    // 修正每一位的进位
    for (size_t i = 0; i < d.v.size(); ++i) {
        carry += d.v[i];
        d.v[i] = carry % BIGINT_BASE;
        carry /= BIGINT_BASE;
    }
    d.trim();
    d.sign = sign * b.sign;
    return d;
}

BigIntSimple operator/(const BigIntSimple &b) const {
    BigIntSimple r;
    return div_mod(b, r);
}

BigIntSimple operator%(const BigIntSimple &b) const {
    BigIntSimple r;
    div_mod(b, r);
    return r;
}

Karatsuba 乘法

记高精度数字的位数为 \(n\)，那么高精度—高精度竖式乘法需要花费 \(O(n^2)\) 的时间。本节介绍一个时间复杂度更为优秀的算法，由前苏联（俄罗斯）数学家 Anatoly Karatsuba 提出，是一种分治算法。

考虑两个十进制大整数 \(x\) 和 \(y\)，均包含 \(n\) 个数码（可以有前导零）。任取 \(0 < m < n\)，记

\[ \begin{aligned} x &= x_1 \cdot 10^m + x_0, \\ y &= y_1 \cdot 10^m + y_0, \\ x \cdot y &= z_2 \cdot 10^{2m} + z_1 \cdot 10^m + z_0, \end{aligned} \]

其中 \(x_0, y_0, z_0, z_1 < 10^m\)。可得

\[ \begin{aligned} z_2 &= x_1 \cdot y_1, \\ z_1 &= x_1 \cdot y_0 + x_0 \cdot y_1, \\ z_0 &= x_0 \cdot y_0. \end{aligned} \]

观察知

\[ z_1 = (x_1 + x_0) \cdot (y_1 + y_0) - z_2 - z_0, \]

于是要计算 \(z_1\)，只需计算 \((x_1 + x_0) \cdot (y_1 + y_0)\)，再与 \(z_0\)、\(z_2\) 相减即可。

上式实际上是 Karatsuba 算法的核心，它将长度为 \(n\) 的乘法问题转化为了 \(3\) 个长度更小的子问题。若令 \(m = \left\lceil \dfrac n 2 \right\rceil\)，记 Karatsuba 算法计算两个 \(n\) 位整数乘法的耗时为 \(T(n)\)，则有 \(T(n) = 3 \cdot T \left(\left\lceil \dfrac n 2 \right\rceil\right) + O(n)\)，由主定理可得 \(T(n) = \Theta(n^{\log_2 3}) \approx \Theta(n^{1.585})\)。

整个过程可以递归实现。为清晰起见，下面的代码通过 Karatsuba 算法实现了多项式乘法，最后再处理所有的进位问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
int *karatsuba_polymul(int n, int *a, int *b) {
    if (n <= 32) {
        // 规模较小时直接计算，避免继续递归带来的效率损失
        int *r = new int[n * 2 + 1]();
        for (int i = 0; i <= n; ++i)
            for (int j = 0; j <= n; ++j) r[i + j] += a[i] * b[j];
        return r;
    }
    int m = n / 2 + 1;
    int *r = new int[m * 4 + 1]();
    int *z0, *z1, *z2;

    z0 = karatsuba_polymul(m - 1, a, b);
    z2 = karatsuba_polymul(n - m, a + m, b + m);

    // 计算 z1
    // 临时更改，计算完毕后恢复
    for (int i = 0; i + m <= n; ++i) a[i] += a[i + m];
    for (int i = 0; i + m <= n; ++i) b[i] += b[i + m];
    z1 = karatsuba_polymul(m - 1, a, b);
    for (int i = 0; i + m <= n; ++i) a[i] -= a[i + m];
    for (int i = 0; i + m <= n; ++i) b[i] -= b[i + m];
    for (int i = 0; i <= (m - 1) * 2; ++i) z1[i] -= z0[i];
    for (int i = 0; i <= (n - m) * 2; ++i) z1[i] -= z2[i];

    // 由 z0、z1、z2 组合获得结果
    for (int i = 0; i <= (m - 1) * 2; ++i) r[i] += z0[i];
    for (int i = 0; i <= (m - 1) * 2; ++i) r[i + m] += z1[i];
    for (int i = 0; i <= (n - m) * 2; ++i) r[i + m * 2] += z2[i];

    delete[] z0;
    delete[] z1;
    delete[] z2;
    return r;
}

void karatsuba_mul(int a[], int b[], int c[]) {
    int *r = karatsuba_polymul(LEN - 1, a, b);
    memcpy(c, r, sizeof(int) * LEN);
    for (int i = 0; i < LEN - 1; ++i)
        if (c[i] >= 10) {
            c[i + 1] += c[i] / 10;
            c[i] %= 10;
        }
    delete[] r;
}

但是这样的实现存在一个问题：在 \(b\) 进制下，多项式的每一个系数都有可能达到 \(n \cdot b^2\) 量级，在压位高精度实现中可能造成整数溢出；而若在多项式乘法的过程中处理进位问题，则 \(x_1 + x_0\) 与 \(y_1 + y_0\) 的结果可能达到 \(2 \cdot b^m\)，增加一个位（如果采用 \(x_1 - x_0\) 的计算方式，则不得不特殊处理负数的情况）。因此，需要依照实际的应用场景来决定采用何种实现方式。