14.复杂度理论与分治

发表于 分类于 ，

复杂度理论与分治

本章探讨了算法分析的核心概念，通过求解最大连续和问题展示了不同算法设计策略的效率差异，并引入了复杂性理论、分治法以及它们在优化算法时间复杂度上的应用。从简单的枚举方法到高效的动态规划与分治策略，逐步揭示了解决问题时算法选择的重要性及其对计算资源消耗的深远影响。

感知复杂度

给一个数字序列，求解一段连续的数字最大的和

输入：序列长度\(n\)，序列\(\{A_1, A_2, \ldots, A_n\}\)，\(A_i\) 为任意整数

输出：找到\(1 \leq i \leq j \leq n\)，使得\(A_i + A_{i+1} + \ldots + A_{j-1} + A_j\)和最大

方法1：枚举所有的连续序列，找和最大的序列

1
2
3
4
5
6
7
8
int best = A[1];                        // 初始化答案
for(int i = 1; i <= n; i ++)            // 枚举起点
    for(int j = i; j <= n; j ++) {      // 枚举终点
        int sum = 0;
        for(int k = i; k <= j; k ++)    // 序列求和
            sum += A[k];
        if(sum > best) best = sum;      // 更新答案
    }

对任意输入个数𝑛，加法（ sum += A[k]; ）的执行次数\(T(n)\)是多少？

  • 第一层for循环 \(\sum_{i=1}^n \ldots\)
  • 第二层for循环\(\sum_{j=i}^n \ldots\)
  • 第三层for循环\(\sum_{k=i}^j 1 = j - i + 1\)

\(T(n) = \sum_{i=1}^n \sum_{j=i}^n (j - i + 1) = \frac{n(n+1)(n+2)}{6}\)

我们取个比较大的\(n=10,000\)，此时\(T(n)=166,716,670,000\)。

方法2：先计算前缀和，前\(j\)个数的和减去前\(i-1\)个数的和可以方便地得到\(i~j\)这一连续序列的和\(A_i + A_{i+1} + \ldots + A_j\)

1
2
3
4
5
6
pre[0] = 0;
for(int i = 1; i <= n; i ++) pre[i] = pre[i - 1] + A[i];
for(int i = 1; i <= n; i ++) {
    for(int j = i; j <= n; j ++)
        best = max(best, pre[j] - pre[i - 1]);
}

此方法\(T(n) = n + \sum_{i=1}^n (n - i + 1) = n + \frac{n(n+1)}{2}\)

\(n=10,000\)时\(T(n)=50,015,000\)

方法3：递归地划分问题为左右两半，看左半边的解（最大连续和）大，还是右半边的解大，还是从中间往左右延申的一段连续和大。

此方法\(T(n) = 2T(n/2) + n\)，中间往左右延申最坏情况就是走到头

\(n=10,000\)时\(T(n) \approx n\log n \approx 130,000\)

方法4：考虑这样一个情况，当连续地累加时，只要和不为负，前面累加的值肯定对后面继续累加有贡献，反之若累加到某个位置为负，则可以抛弃前面的累加值。

此方法\(T(n) = n\)，只需要从前往后每个数考虑一次

\(n=10,000\)时\(T(n)=10,000\)

算法不同，处理同一个问题的基础运算次数可能有非常大的差距

复杂性理论

算法更关心随着\(n\)无限增长，\(T(n)\)的增长速度——渐进效率

通常称渐进效率的高低为“阶”的高低，阶更高的函数增长率更高。

表达渐进效率有 5 个符号：

  • 大写\(O\)，\(f(n) = O(g(n))\)，表示\(f(n)\)的阶不大于\(g(n)\)的阶；
  • 大写\(\Omega\)，\(f(n) = \Omega(g(n))\)，表示\(f(n)\)的阶不小于\(g(n)\)的阶；
  • 小写\(o\)，\(f(n) = o(g(n))\)，表示\(f(n)\)的阶严格小于\(g(n)\)的阶；
  • 小写\(\omega\)，\(f(n) = \omega(g(n))\)，表示\(f(n)\)的阶严格大于\(g(n)\)的阶；
  • 大写\(\Theta\)，\(f(n) = \Theta(g(n))\)，表示\(f(n)\)的阶与\(g(n)\)的阶相同。

最常用的是大\(O\)表示，注意不大于这个描述，即 \(n=O(n^2)\)、\(n^2+n=O(n^2)\) 这类表述都是正确的，大 \(O\) 关心的是一个算法复杂度的“上界”，至于是不是严格等于这个上界就不那么重要了，这样更方便评估一个算法的效率是否“可用”。

算法复杂度\(T(n)\)是关于输入规模\(n\)的函数，如果这个函数包含多项，比如 \(5^{n}+2n^{2}+3n\dots\)，大\(O\)表示法只保留阶最高的项，且不关心常数系数，这个例子的大\(O\)表示应该简化为 \(O(5^{n})\)

空间复杂度

算法启动后需要的额外空间。

通常先指定一个常数级别的基本单位，比如 int（4个字节） 大小的内存、长度为10的char等等，根据具体的算法指定一个合乎逻辑的基本单位，统计基本单位的使用数量，得到一个关于输入规模 \(n\) 的函数。

时间复杂度

算法基本运算执行的次数，基本运算是一个或一组有限时间内结束的、算法中执行频率最高的运算，比如加法、乘法、比较、常数级别的矩阵乘法等等。

分治

递归地调用自身解决紧密相关的若干子问题，步骤如下：

  1. 分解原问题为若干互不相交的子问题
    这些子问题是原问题规模较小的实例，与原问题属性相同
  2. 解决这些子问题，递归进行
    如果子问题规模足够小则直接求解
  3. 合并这些子问题的解为原问题的解
    递归的各子问题解决完后回到上一层进行合并

以归并排序为例

  1. 分解：待排序的\(n\)个元素序列分解为两个\(n/2\)个元素的序列
  2. 解决：递归地对两个子序列执行归并排序，长度为1则直接返回
  3. 合并：“回升”阶段将两个已排序的子序列合并为一个已排序的结果
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void MergeSort(int a[], int l, int r) {
    if(l >= r - 1) return;     // 【解决】最小规模
    int mid = l + r >> 1;
    MergeSort(a, l, mid);      // 【分解】并【解决】
    MergeSort(a, mid, r);      // 【分解】并【解决】
    // 【合并】阶段开始
    int i = l, j = mid, tp = l;
    while(i < mid && j < r) {
        if(a[j] < a[i]) tmp[tp ++] = a[j ++];
        else tmp[tp ++] = a[i ++];
    }
    while(i < mid) tmp[tp ++] = a[i ++];
    while(j < r) tmp[tp ++] = a[j ++];
    memcpy(a + l, tmp + l, sizeof(int) * (r - l));
    // 【合并】阶段结束，【解决】了当前子问题[l,r)
}

分治的时间复杂度分析

归并排序要解决 \(2\) 个 \(n/2\) 规模的子问题，把 \(2\) 个子问题的结果合并为一个有序数组有 \(O(n)\) 的代价，复杂度函数表示为

\(T(n) = 2T( n/2 ) + O(n)\)

这种原始问题（规模为\(n\)）的复杂度\(T(n)\)由拆分的若干个（\(a\)个）、相同方法解决的子问题（\(n\)减小的某个规模）复杂度函数\(T(n/b)\)和合并子问题的复杂度\(f(n)\)组成的式子称为递归式，求得由 \(n\) 的公式直接表达 \(T(n)\) 的结果的过程称为解递归式，递归式的解就是这个分治算法的时间复杂度。典型形式为：

多个相同规模子问题：\(T(n)=aT(n/b)+f(n)\)， 多个不同规模子问题：\(T(n)=T(n/b)+T(n/c)+\dots+f(n)\)

最直观的方法是通过递归树求解：

每个结点表示单一子问题合并工作的代价（\(f(n)\)），树中每层代价求和，所有层的代价求和

对于\(T(n) = 2T(n/2) + n\)，

此处我们仍假设\(n\)是2的幂

例：线性（平均情况）第k小

\(n\)个未排序的数，选择其中第\(k\)小的数。线性算法指算法复杂度 \(T(n)\) 是关于 \(n\) 的一次函数。

假设数据规模不允许基数排序等特殊排序方法，则快速排序、归并排序等方式需要 \(nlogn\)，不满足线性。

设数组为\(S\)，以“快速排序”的思想进行划分

  1. 分解：选一个数\(m^{*}\)，按大小划分为\(S_1\)和\(S_2\)（\(m^{*}\)在两集合之间），即数组所有比 \(m^{*}\) 小的放它左边，比 \(m^{*}\) 大的放它右边。
  2. 解决：
    • 若\(k \leq |S_1|\)，则在\(S_1\)中找第\(k\)小
    • 若\(k > |S_1| + 1\)，则在\(S_2\)中找第\(k - |S_1|\)小
    • 若\(k = |S_1| + 1\)，则第\(k\)小为\(m^{*}\)
  3. 合并：第\(k\)小的数只会存在于\(S_1\)、\(S_2\)、\(m^{*}\)其中一个，在这一边递归找

分析：平均情况下，子问题规模（左半边或右半边）\(n/2\)，因为只会在一边找，所以只需要解决 \(1\) 个子问题，合并子问题的工作集中在划分这件事上，即“比 \(m^{*}\) 小的放它左边，比 \(m^{*}\) 大的放它右边”，这个操作是 \(O(n)\) 的，所以递归式为：

\(T(n)=T(n/2)+n\)

解递归式得 \(T(n)=O(n)\) ，在平均情况下是线性的。

为什么要说“平均情况”呢？因为选\(m^{*}\)划分的时候，有可能大多数都分到了一侧，这种不均衡会导致\(T(n)\)不再是线性的，但选择\(m^{*}\)的随机性很高，可以认为较大概率不会很偏，或从统计角度讲，多次执行的平均运算量是线性的。

其实这个问题有一个理论上更强大的解法，在最坏情况下也是\(O(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
36
37
38
39
40
41
42
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<cmath>
#include<complex>
#include<vector>
#include<algorithm>

const int mod = 1e9 + 7;
const int maxn = 1e8 + 10;
int n, m, k;
int a[maxn];

int GetK(int left, int right, int k)
{
    if(left == right - 1) return a[left];
    int low = left, high = right - 1, center = a[low];
    while(low < high)
    {
        while(low < high && a[high] >= center) high --;
        a[low] = a[high];
        while(low < high && a[low] <= center) low ++;
        a[high] = a[low];
    }
    a[low] = center;
    if(low - left >= k) return GetK(left, low, k);
    else if(low - left + 1 == k) return a[low];
    else return GetK(low + 1, right, k - (low - left) - 1);
}
int main()
{

    while(scanf("%d%d%d", &n, &m, &k) != EOF)
    {
        a[0] = m;
        for(int i = 1; i < n; i ++)
            a[i] = 1LL * a[i - 1] * m % mod;
        printf("%d\n", GetK(0, n, k));
    }

    return 0;
}

例：平面最近点对

平面有\(n\)个点，位置以\((x, y)\)坐标表示，求最近的两个点及其距离。

设两个子问题最优解距离为\(d\)

每个小矩形内至多1个点

  • 考虑小矩形内两点最远距离（即对角线\(r\)）
  • \(r=\sqrt{((d/2)^2+(2d/3)^2)}=5d/6<d\)，
  • 如有多个点则与子问题最优解矛盾

合并时，左侧任一点\(P\)在右侧至多考虑6个点

\(T(n)=2T(n/2)+O(n \log n)\)，\(T(n)=O(n \log^2 n)\)

虽然左边每个点考查右边 \(6\) 个点应当是 \(O(n)\) 的，但是要快速找到这 \(6\) 个点，就不得不排序，所以合并子问题的工作量稍大了些。

但可以进一步优化，尝试一次排序，在分治的时候同时划分排序信息，避免每次都排序，自行尝试。

参考代码

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
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<vector>
const int maxn = 1e5 + 10;
const int inf = 0x3f3f3f3f;     // 定义一个极大值
struct Point {
    int x, y, idx;              // 坐标及序号
    bool operator<(const Point &b){return x < b.x;} // 重载“<”用于排序
};
Point p[maxn], csh[maxn];
inline int Sqr(int x){return x * x;}
int Dis(Point &a, Point &b){
    return Sqr(a.x - b.x) + Sqr(a.y - b.y);
}
int n, ansa, ansb, ansdis;
int Update(Point &a, Point &b) {
    // 由于答案要输出点对编号，而不仅仅是距离，用此函数更新答案
    int na = a.idx, nb = b.idx, dis = Dis(a, b);
    if(na > nb) std::swap(na, nb);
    if(ansdis == -1 || dis < ansdis || dis == ansdis && (na < ansa || na == ansa && nb < ansb))
        ansa = na, ansb = nb, ansdis = dis;
    return dis;
}
void Merge(int l, int r, int mid) {
    int i = l, j = mid, k = l;
    while(i < mid && j < r)
    {
        if(p[i].y < p[j].y) csh[k ++] = p[i ++];
        else csh[k ++] = p[j ++];
    }
    while(i < mid) csh[k ++] = p[i ++];
    while(j < r) csh[k ++] = p[j ++];
    memcpy(p + l, csh + l, sizeof(Point) * (r - l));
}
void JudgeNearby(int l, int mid, int r, Point &c) {
    // 注意到本算法采用了一个 trick：回溯的同时对 y 坐标做归并排序
    // 说明处理 [l, r) 区间时，[l, mid)和[mid, r)两个子问题已经 y 坐标有序
    // 提取两个子问题到中间点的 x 距离不超过子问题最优解的那些点 pl[] 和 pr[]
    // 只需左右“异步前进”，就能很容易定位左边 pl[] 每个点在右边上下距离不超过当前最优解的那些点
    std::vector<Point> pl, pr;
    for(int i = l; i < mid; i ++)
        if(Sqr(c.x - p[i].x) <= ansdis) pl.push_back(p[i]);
    for(int i = mid; i < r; i ++)
        if(Sqr(c.x - p[i].x) <= ansdis) pr.push_back(p[i]);
    for(int i = 0, j = 0; i < pl.size(); i ++) {    // 对于左侧 pl[] 的每个点
        // 找到右侧 pr[] 的点中距离 pl[i] 小于目前最优解的序号下界 j
        while(j < pr.size() && pl[i].y > pr[j].y && Sqr(pl[i].y - pr[j].y) > ansdis) j ++;
        // 从 j 开始枚举 pr[] 中的点，直到 y 坐标距离超过当前最优解
        for(int k = j; k < pr.size() && Sqr(pl[i].y - pr[k].y) <= ansdis; k ++)
            Update(pl[i], pr[k]);
    }
}
int MinDis(int l, int r) {
    if(r - l <= 1) return inf;          // 递归终点：子问题中没有点，返回极大距离
    int mid = l + r >> 1;
    Point midp = p[mid];
    MinDis(l, mid);                     // 处理子问题最近点对
    MinDis(mid, r);                     // 处理子问题最近点对
    JudgeNearby(l, mid, r, midp);       // 两侧子问题取相近点求距离更新答案
    Merge(l, r, mid);                   // 对 y 坐标进行归并排序的 Merge 操作
    return ansdis;
}
int main() {
    while(scanf("%d", &n) != EOF) {
        ansa = ansb = ansdis = inf;
        for(int i = 0; i < n; i ++) {
            scanf("%d%d", &p[i].x, &p[i].y);
            p[i].idx = i;
        }
        std::sort(p, p + n);    // 按重载的“<”，以x从小到大排序
        MinDis(0, n);
        printf("%d %d\n", ansa, ansb);
    }
    return 0;
}

分治优化策略

用一部分子问题的解表达另一部分子问题，减少计算子问题的个数

例：快速幂

\(a\)为给定实数，\(n\)为自然数，求$a^n

\(a\)为给定实数，\(n\)为自然数，求\(a^n\)

  1. 分解：两份\(n/2\)个\(a\)相乘
  2. 解决：对于奇数个的子问题，剩下的1个保留，乘在最终结果上
  3. 合并：两个子问题的解相乘

但两个子问题结果相同，不必递归计算两遍

\(T(n) = T(n/2) + \Theta(1)\)，\(T(n) = O(\log n)\)

奇数个数会有“零头”，这个零头令开一个变量“屯”起来。

问题比较简单，可以“自底向上”地合并，不需要递归了。

1
2
3
4
5
6
int PowMod(int a, int b, int mod) {
    int ret = 1;
    for(;b; b >>= 1, a = 1LL * a * a % mod)
        if(b & 1) ret = 1LL * ret * a % mod;
    return ret;
}

例：快速幂计算斐波那契数列第 \(n\) 项

设\(f_0=0\), \(f_1=1\), \(f_n=f_{n-1}+f_{n-2}\)，放入2×2的矩阵，发现：

\[ \begin{pmatrix} f_{n+1} & f_n \\ f_n & f_{n-1} \end{pmatrix} \times \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} f_{n+1} + f_n & f_{n+1} \\ f_n + f_{n-1} & f_n \end{pmatrix} = \begin{pmatrix} f_{n+2} & f_{n+1} \\ f_{n+1} & f_n \end{pmatrix} \]

即矩阵每一项都变成了斐波那契数列对应位置的下一项。

从而，

\[ \begin{pmatrix} f_{n+1} & f_n \\ f_n & f_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n \]

把矩阵类比实数，同样可以用快速幂方法求斐波那契数列第\(n\)项

\(T(n) = O(\log n)\)

例：Strassen矩阵乘法

求\(n\)阶矩阵\(A\)与\(B\)相乘的结果\(AB\)，其中假设\(n=2^k\)

直接相乘方法：\(T(n)=O(n^3)\)

直接分治： 分解： \[ \begin{bmatrix} A_{11} & A_{12}\\ A_{21} & A_{22} \end{bmatrix} \begin{bmatrix} B_{11} & B_{12}\\ B_{21} & B_{22} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12}\\ C_{21} & C_{22} \end{bmatrix} \]

求解：递归按矩阵乘规则求各子矩阵相乘结果

合并： \[C_{11}=A_{11}B_{11}+A_{12}B_{21}，C_{12}=A_{11}B_{12}+A_{12}B_{22}，\] \[C_{21}=A_{21}B_{11}+A_{22}B_{21}，C_{22}=A_{21}B_{12}+A_{22}B_{22}.\]

\(T(n)=8T(n/2)+O(n^2)，T(n)=O(n^3)\)——这样不行

子问题定义：两个\(n/2\)规模的子矩阵相乘

优化合并：将\(A_{11}B_{11}\)等8个子问题通过相互表示优化为7个

设\(M_1 \sim M_7\)为7个子问题

\[M_1 = A_{11}(B_{12} - B_{22})\] \[M_2 = (A_{11} + A_{12})B_{22}\] \[M_3 = (A_{21} + A_{22})B_{11}\] \[M_4 = A_{22}(B_{21} - B_{11})\] \[M_5 = (A_{11} + A_{22})(B_{11} + B_{22})\] \[M_6 = (A_{12} - A_{22})(B_{21} + B_{22})\] \[M_7 = (A_{11} - A_{21})(B_{11} + B_{12})\]

\[C_{11} = M_5 + M_4 - M_2 + M_6,\] \[C_{12} = M_1 + M_2,\] \[C_{21} = M_3 + M_4,\] \[C_{22} = M_5 + M_1 - M_3 - M_7.\]

即Strassen矩阵乘法用7个子问题完成了\(C\)矩阵的计算

\[T(n) = 7T(n/2) + O(n^2)，T(n) = O(n^{\log 7}) \approx O(n^{2.81})\]

Coppersmith-Winograd算法为\(O(n^{2.3727})\)

例：快速傅里叶变换

使用FFT解决多项式乘法问题

设两个多项式\(F(x)\)和\(G(x)\)如下：

\[ F(x) = a_0 + a_1x + a_2x^2 + \ldots \]

\[ G(x) = b_0 + b_1x + b_2x^2 + \ldots \]

它们的乘积可以表示为：

\[ F(x)G(x) = \sum_{i=0} \sum_{j=0}^{i} a_j b_{i-j} x^i \]

应用到大整数乘法

如果把\(x\)换成10，就解决了大整数乘法问题。例如：

\[ 314159 = 9 + 5 \times 10 + 1 \times 10^2 + 4 \times 10^3 + 1 \times 10^4 + 3 \times 10^5 \]

对多项式 \(y = a_0 + a_1x + a_2x^2 + \ldots + a_{n-1}x^{n-1}\) 做一个反向思考：

A同学有多项式，知道所有系数\(a_0, a_1, \ldots, a_{n-1}\)（都是已知常数）

B同学只知道多项式的形式，但不知道多项式具体是什么

A同学向多项式代入\(n\)个不同的\(x_i\)，算出\(n\)个不同的\(y_i\)

得到\((x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n)\)，可以视作“点”

思考：B同学拿到\(n\)“对”\((x_i, y_i)\)，他能否反向算出\(a_0, a_1, \ldots, a_{n-1}\)的值？

\(a_0, a_1, \ldots a_{(n-1)}\)与\((x_1,y_1),(x_2,y_2),\ldots(x_n,y_n)\)都唯一确定一个多项式

系数表示法： \(a_0,a_1,\ldots a_{(n-1)}\)

点值表示法： \((x_1,y_1),(x_2,y_2),\ldots(x_n,y_n)\)

系数表示法的多项式乘法是 \[ F(x)G(x)=\sum_{i=0} \sum_{j=0}^i a_j b_{(i-j)} x^i \]

思考：点值表示法的多项式乘法如何进行？

对每个特定 \(x\) ，两个点表示法的多项式相乘，就是 \(y\) 的值相乘。\(F(x)G(x)\)的点值表示法是\((x_1, (y_{a1} * y_{b1})), (x_2, (y_{a2} * y_{b2})), \ldots\)

如何利用点值表示法加速多项式乘法？直接代入\(n\)个不同的\(x\)并不能实质降低复杂度，仍然要\(O(n^2)\)得到所有点值：

  • \(x_1\) 代入\(F(x)=a_0+a_1 x+\ldots a_{(n-1)} x^{(n-1)}\)得到\(y_1\)需要\(O(n)\)时间
  • \(n\)个不同的\(x_i\)分别代入，总共就要\(O(n^2)\)的时间

设计特殊的\(n\)个不同的\(x\)，加速每个点的计算。

复数域下，\(x^n=1\)的\(x\)是\(n\)次单位根 \[x_k=e^{i 2k\pi/n}, \quad k=0,1,2,\ldots n-1\]

用它作为\(n\)个不同的\(x\)，相当于等分了圆周角 \[x_k=e^{i 2k\pi/n}=e^{i 2\pi/n k}=\omega_n^k\]

欧拉公式\(e^{i2k\pi}=\cos 2k\pi + i \sin 2k\pi = 1\)

\(e^{i2k\pi}\)开\(n\)次根\(e^{i 2k\pi/n}\)就是一个解

\(F(x) = F(\omega_n) = a_0 + a_1 \omega_n + \ldots + a_{(n-1)} \omega_n^{(n-1)}\)

奇偶系数分开定义新多项式

\[A(\omega_n) = a_0 + a_2 \omega_n^1 + a_4 \omega_n^2 + \ldots + a_{(n-2)} \omega_n^{(n/2-1)}\] \[B(\omega_n) = a_1 + a_3 \omega_n^1 + a_5 \omega_n^2 + \ldots + a_{(n-1)} \omega_n^{(n/2-1)}\] \[F(\omega_n) = A(\omega_n^2) + \omega_n B(\omega_n^2)\]

\(F(\omega_n) = A(\omega_n^2) + \omega_n B(\omega_n^2)\)，分别代入\(\omega_n^k\)和\(\omega_n^{(k+n/2)}\)

\(F(\omega_n^k) = A(\omega_n^{2k}) + \omega_n^k B(\omega_n^{2k})\) \(= A(\omega_{(n/2)}^k) + \omega_n^k B(\omega_{(n/2)}^k)\)

\(F(\omega_n^{(k+n/2)}) = A(\omega_n^{(2k+n)}) + \omega_n^{(k+n/2)} B(\omega_n^{(2k+n)})\) \(= A(\omega_{(n/2)}^k) - \omega_n^k B(\omega_{(n/2)}^k)\)

\(F(\omega_n^k) = A(\omega_{(n/2)}^k) + \omega_n^k B(\omega_{(n/2)}^k)\)

\(F(\omega_n^{(k+n/2)}) = A(\omega_{(n/2)}^k) - \omega_n^k B(\omega_{(n/2)}^k)\)

分治！

解决“多项式代入求值”的子问题\(A(\omega_{(n/2)}^k)\)和\(B(\omega_{(n/2)}^k)\) 就能求出\((x_k, y_k)\)，\((x_{(k+n/2)}, y_{(k+n/2)})\)两个点

\(A(\omega_{(n/2)}^k)\)和\(B(\omega_{(n/2)}^k)\)的规模（多项式系数个数）都是原问题的一半

蝴蝶变换

\[A(\omega_n) = a_0 + a_2 \omega_n^1 + a_4 \omega_n^2 + \ldots + a_{(n-2)} \omega_n^{(n/2-1)}\] \[B(\omega_n) = a_1 + a_3 \omega_n^1 + a_5 \omega_n^2 + \ldots + a_{(n-1)} \omega_n^{(n/2-1)}\]

观察系数序号与原始多项式的关系，在继续分治过程中，序号变化：

  • 初始序号：0, 1, 2, 3, 4, 5, 6, 7
  • 第一次分组：0, 2, 4, 6 | 1, 3, 5, 7
  • 第二次分组：0, 4 | 2, 6 | 1, 5 | 3, 7

能否直接得到递归终点的序号序列
0,4 | 2,6 | 1,5 |3,7 的二进制：
000, 100, 010, 110, 001, 101, 011, 111
所有二进制反向：
000, 001, 010, 011, 100, 101, 110, 111

蝴蝶变换代码参考

1
2
3
4
5
6
7
8
9
10
void BitRevChange(Complex y[], int len) {
    // len should be 2^k
    std::vector<int> rev(len, 0);
    for (int i = 0; i < len; i ++) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) rev[i] |= len >> 1;
    }
    for (int i = 0; i < len; ++i) 
        if (i < rev[i]) std::swap(y[i], y[rev[i]]);
}

大整数乘法参考代码

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
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<cmath>
#include<complex>
#include<vector>
#include<algorithm>
const double pi = acos(-1);
const double eps = 1e-6;
const int maxn = 4e6 + 10;
typedef std::complex<double> Complex;
char a[maxn], b[maxn];
Complex xa[maxn], xb[maxn];
int lena, lenb, len2;
void BitRevChange(Complex y[], int len) 
{
    std::vector<int> rev(len, 0);
    for (int i = 0; i < len; i ++) 
    {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) rev[i] |= len >> 1;
    }
    for (int i = 0; i < len; ++i) 
        if (i < rev[i]) std::swap(y[i], y[rev[i]]);
    return;
}
void FFT(Complex y[], int len, int on=1)
{
    BitRevChange(y, len);
    for(int h = 2; h <= len; h <<= 1)
    {
        Complex wn(cos(2 * pi / h), on * sin(2 * pi / h));
        for(int j = 0; j < len; j += h)
        {
            Complex w(1, 0);
            for(int k = j; k < j + h / 2; k ++)
            {
                Complex u = y[k];
                Complex t = w * y[k + h / 2];
                y[k] = u + t;
                y[k + h / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if(on != -1) return;
    for(int i = 0; i < len; i ++) y[i].real(y[i].real() / len);
}
int main()
{
    while(scanf("%s%s", a, b) != EOF)
    {
        lena = strlen(a); lenb = strlen(b);
        for(int i = 0; a[i]; i ++) xa[i] = a[lena - i - 1] - '0';
        for(int i = 0; b[i]; i ++) xb[i] = b[lenb - i - 1] - '0';
        for(len2 = 1; len2 < lena + lenb; len2 <<= 1);
        for(int i = lena; i < len2; i ++) xa[i] = 0;
        for(int i = lenb; i < len2; i ++) xb[i] = 0;
        FFT(xa, len2); FFT(xb, len2);
        for(int i = 0; i < len2; i ++) xa[i] *= xb[i];
        FFT(xa, len2, -1);
        for(int i = 0; i < len2; i ++)
        {
            a[i] = (int)(xa[i].real() + 0.5) % 10 + '0';
            xa[i + 1].real(xa[i + 1].real() + (int)(xa[i].real() + 0.5) / 10);
        }
        for(lena = len2 - 1; a[lena] == '0' && lena > 0; lena --);
        for(; lena >= 0 && printf("%c", a[lena]); lena --);
        printf("\n");
    }
    return 0;
}