DNF【题解】A + B Problem 题解
DNF（Do Not Finish）
就不作为题解提交了，不仅可能会被骂，同时作为题解也确实没有意义。
零、前言
参考：
具体学习并实现快速傅里叶变换 - 鹤翔万里
洛谷日报 71：傅里叶变换 (FFT) 学习笔记 - command_block
洛谷日报 4：浅谈线段树 - 皎月半洒花
在这里特别感谢。
这篇文章只是一篇题（che）解（dan），也可以看作是一篇线段树、高精度加法、快速傅里叶变换的融合题解。
如果题解没过我也没有损失，因为还是学到了很多其他的东西，但是还是过一下吧qwq.
或者直接做日报，我看你们已经好久没更新日报了，
由于我只会模板类题目，写这个就只是... 复习基础知识点，加深印象，然后瞎搞，
我 FFT 也就学了一个模板，当然，就模板也忘了差不多了，
可能真的没什么意义，我就是无聊...
我不会 FFT。
由于求解 A + B Problem 较为困难，导致本篇文章篇幅较长，希望大家耐心观看。
在半期考温书假时，我发现了前一天在我们学校 OJ 上面写的 A + B Problem 题解被驳回了，理由：不够详细。
/fn
我非常气愤，决定详解 A + B Problem。
全篇使用 C++ 语言编写代码。
保证代码通过 C++14 编译。
使用万能头 #include。
对于代码中的 ll，为 long long，在最开始有加上 typedef long long ll;。
Pi 的定义：const double Pi = acos(-1.0);。
一、题目大意
题目链接：洛谷 - P1001 A+B Problem。
给定整数 \(a\) 和 \(b\)，求 \(a + b\)，其中 \(|a|,|b| \le 10^9\)。
二、解法
我们可以规定一个长度为 \(2\) 区间，然后把数分别放进去，最后求区间和即可。
线段树可以修改区间以及查询区间。
于是我们可以想到使用线段树来求解这种问题。
三、高精度加法
由于线段树会使用到加法，所以我们先来看看高精度加法。
高精度加法只能做两非负整数相加，于是我们只要在出现负数的时候分类讨论，使用高精度减法，或者 C++ 自己的加法运算。
1. 概念与实现
高精度加法就类似于我们平时手算加法一样，按位相加。
按位相加，每一位使用字符保存，保存进位，在下一次相加时加入，两数长度不一样，可以前补零。
就是一个模拟的过程，比较简单，这里给出代码。
例：洛谷 - P1601 A+B Problem（高精）。
ll add(ll number_1, ll number_2) {
string n1_number = to_string(number_1), n2_number = to_string(number_2);
ll length_of_n1 = n1_number.length();
ll length_of_n2 = n2_number.length();
ll max_of_string_length = max(length_of_n1, length_of_n2);
if(length_of_n1 < max_of_string_length)
for(ll i = length_of_n1; i < max_of_string_length; i++)
n1_number = "0" + n1_number;
else for(ll i = length_of_n2; i < max_of_string_length; i++)
n2_number = "0" + n2_number;
string add_sum = "";
ll last_num = 0, add_sum_1 = 0, add_sum_2 = 0, add_sum_add = 0;
for(ll i = max_of_string_length - 1;i >= 0; i--) {
add_sum_1 = n1_number[i] - '0';
add_sum_2 = n2_number[i] - '0';
add_sum_add = add_sum_1 + add_sum_2 + last_num;
add_sum = to_string(add_sum_add % 10) + add_sum;
last_num = add_sum_add / 10;
}
if(last_num) add_sum = to_string(last_num) + add_sum;
return atoll(add_sum.c_str());
}
四、快速傅里叶变换前言
由于本题线段树中有用到乘法，所以我们使用快速傅里叶变换做乘法。
1. 概念
快速傅里叶变换（Fast Fourier Transform)，即利用计算机计算离散傅里叶变换（DFT）的高效、快速计算方法的统称，简称 FFT。快速傅里叶变换是 1965 年由 J.W.库利 和 T.W.图基 提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数 \(N\) 越多，FFT 算法计算量的节省就越显著。
以上摘自 百度百科 - 快速傅里叶变换。
2. 大数乘法
快速傅里叶变换可以在 \(O(n\log{n})\) 的时间内快速地计算两个一元多项式的卷积。
对于一个整数，可以被分解为 \(a_1+a_2x+a_3x^2+a_4x^3+\cdots\)，其中 \(x = 10\)，这是一个多项式。
若有两个整数，分解成两个一元多项式，计算他们的卷积，最后相加，即为他们的积。
3. 前置知识
为了学习快速傅里叶变换，我们先来学习以下前置知识：
- 复数
- 单位根
- 一元多项式
现在开始讲一下前置知识。
五、复数
1. 概念
复数（复杂的数，complex number），是数的概念扩展。
我们把形如 \(z=a+bi\)（\(a\)、\(b\)均为实数）的数称为复数。其中，\(a\) 称为实部，\(b\) 称为虚部，\(i\) 称为虚数单位。当 \(z\) 的虚部 \(b=0\) 时，则 \(z\) 为实数；当 \(z\) 的虚部 \(b\neq0\) 时，实部 \(a=0\) 时，常称 \(z\) 为纯虚数。复数域是实数域的代数闭包，即任何复系数多项式在复数域中总有根。
以上摘自 百度百科 - 复数。
虚数单位 \(i=\sqrt{-1}\)。
2. 四则运算
复数加法
\[(a+bi)+(c+di)=a+bi+c+di=(a+c)+(bi+di)=(a+c)+(b+d)i \]复数减法
\[(a+bi)-(c+di)=a+bi-c-di=(a-c)+(bi-di)=(a-c)+(b-d)i \]复数乘法
\[(a+bi)(c+di)=ac+adi+bci+bi\cdot di=ac+adi+bci-bd=(ac-bd)+(ad+bc)i \]复数除法
\[\dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{(ac+bd)+(-ad+bc)i}{(c^2+d^2)+(-cd+cd)i}=\dfrac{ac+bd}{c^2+d^2}+\dfrac{bc-ad}{c^2+d^2}i \]3. 欧拉定理
复变函数中，\(e^{ix}=\cos{x}+i\sin{x}\) 称为欧拉公式，\(e\) 是自然对数的底，\(i\) 是虚数单位。
当 \(x=\pi\) 时
\[e^{i\pi}=\cos{\pi}+i\sin{\pi}=-1+0\\e^{i\pi}+1=0 \]即欧拉恒等式。
好的扯远了，我们回来。
4. 三角形式
把 \(x\) 换一下，变成 \(\theta\)，方便一点。
明显地，\(e^{i\theta}=\cos{\theta}+i\sin{\theta}\)，把一个复数写成了 \(e\) 的指数的形式，于是一个复数也可以写成 \(z=re^{i\theta}\)。
像这样的表示形式，被称为三角形式。
其中 \(r\) 是 \(z\) 的模，即 \(r = |z|\)；\(θ\) 是 \(z\) 的辐角，在区间 \([-\pi,\pi]\) 内的辐角称为辐角主值，记作 \(\arg(z)\)。
三角形式下的复数乘法
根据乘方的运算法则，我们可以得到：
\[ae^{i\alpha}\cdot be^{i\beta}=abe^{i(\alpha+\beta)} \]5. 形式互换
\(a+bi=r(\cos{\theta}+i\sin{\theta})=re^{i\theta}\)，其中 \(r=\sqrt{a^2+b^2}\)，\(\cos{\theta}=\dfrac{a}{r}\)，\(\sin{\theta}=\dfrac{b}{r}\)
6. 几何意义
复数的几何意义
由于有实部和虚部，我们可以在一个平面直角坐标系中表示一个复数，\(y\) 轴作为虚轴，\(x\) 作为实轴。而这个平面，叫做复数平面，简称复平面。
若用三角形式表示，复数所对应的向量长度称为复数的模，该向量与实轴正方向的夹角为复数的辐角。
于是我们可以更好地理解形式互换公式了。
复数加法的几何意义
设 \(A=1+i,B=2+i\)，
\(C=A+B=3+2i,D=A+C=4+3i,E=B+D=6+4i\)。
我们把 \(A,B,C,D,E\) 在复平面上表示出来。
发现规律了吗？
同向量加法，我们把每一个复数都看成向量，然后做向量加法。
对于每一个复数，加上另一个复数，相当于以第一个加数为原点重新画复平面，然后在这个新的复平面上画第二个加数，然后回来，即为和。
上述例子是虚部实部都为正数的情况，对于其他的情况也是一样的，大家可以自行探究。
三角形式下的复数乘法的几何意义
设 \(A=1+i,B=2+i\)，
\(C=AB=1+3i,D=AC=-2+4i\)。
我们把 \(A,B,C,D\) 在复平面上表示出来。
我们把复数换成三角形式，然后就会发现，复数乘积的模，即为乘数的模之积，复数乘积的辐角，即为乘数的辐角之和。
诶，你问我为什么不给出上述复数的三角形式？
你可以自己算，根据形式互换公式（其实是因为这玩意儿是无理数）。
7. 代码实现
C++ 有自带的复数头文件 complex，属于 STL，可以使用但是慢，所以还是建议手写复数。
struct Complex_Number {
double real, imag;
Complex_Number() {}
Complex_Number(double num1, double num2) : real(num1) , imag(num2) {}
};
Complex_Number operator + (const Complex_Number &num1, const Complex_Number &num2) {
return Complex_Number(num1.real + num2.real, num1.imag + num2.imag);
}
Complex_Number operator - (const Complex_Number &num1, const Complex_Number &num2) {
return Complex_Number(num1.real - num2.real, num1.imag - num2.imag);
}
Complex_Number operator * (const Complex_Number &num1, const Complex_Number &num2) {
return Complex_Number(num1.real * num2.real - num1.imag * num2.imag, num1.real * num2.imag + num1.imag * num2.real);
}
五、单位根
1. 概念
数学上，\(n\) 次单位根（unit root）是 \(n\) 次幂为 \(1\) 的复数。它们位于复平面的单位圆上，构成正 \(n\) 边形的顶点，其中一个顶点是 \(1\)。
以上内容改编自 百度百科 - 单位根。
简单来说单位根就是方程 \(z^n=1 \ (n=1,2,3,\cdots)\) 在复数范围内的 \(n\) 个根。
这方程的复数根 \(z\) 为 \(n\) 次单位根。
\[z^n=r^ne^{n\theta i}=1 \]因为单位根的 \(n\) 次方为 \(1\)，所以单位根的模一定为 \(1\)。
因为单位根的 \(n\) 次方为 \(1\)，所以单位根的辐角的 \(n\) 倍都是 \(2\pi\) 的倍数。
\[\left\{\begin{matrix} r = 1 \\ n\theta = 2\pi k \end{matrix}\right. \]单位的 \(n\) 次根有 \(n\) 个，我们可以找到这些点：
\[\omega_n^k=e^{i\frac{2\pi k}{n}}=\cos{\frac{2\pi k}{n}}+i\sin{\frac{2\pi k}{n}} \ \ \ (k = 0,1,2,\cdots,n-1) \]每一个单位根都均匀地落在单位圆上，如图（\(8\) 次单位根）：
同时每一个单位根都可以看作 \(\omega_n = e^{i\frac{2\pi}{n}}\) 的幂，所以这个单位根也被称作主 \(n\) 次单位根，记作 \(\omega_n\)。
单位根有它的性质，这里有 \(3\) 个重要的性质
2. 消去引理
基本
\[\omega_{dn}^{dk}=\omega_n^k \]证明：
\[\omega_{dn}^{dk}=(e^{i\frac{2\pi}{dn}})^{dk}=(e^{i\frac{2\pi}{n}})^{k}=\omega_n^k \]几何意义
\(z\) 是 \(4\) 次单位根，\(a\) 是 \(8\) 次单位根。
应该可以看出来吧...
3. 折半引理
基本
\[(\omega_n^{k+\frac{n}{2}})^2=(\omega_n^k)^2=\omega_\frac{n}{2}^k \]证明:
\[\omega_n^{k+\frac{n}{2}}=\omega_n^k\omega_n^\frac{n}{2}=-\omega_n^k,\ (-\omega_n^k)^2=\ (\omega_n^k)^2=\omega_n^{2k}=\omega_{\frac{n}{2}}^k \]几何意义
还是拿这张图：
\(z\) 是 \(4\) 次单位根，\(a\) 是 \(8\) 次单位根。
\((\omega_n^{k+\frac{n}{2}})^2\) 相当于，先绕半圈，再绕半圈，但是会多一位。
\((\omega_n^k)^2\) 相当于 \(\omega_n^k\) 的下一个单位根。
然后使用消去引理。
4. 求和引理
基本
\[\sum_{i = 0}^{n - 1}{(\omega_n^k)^i} = 0 \]证明：
\[\sum_{i = 0}^{n - 1}{(\omega_n^k)^i} = \dfrac{(\omega_n^k)^n - 1}{\omega_n^k - 1} = \dfrac{(\omega_n^n)^k - 1}{\omega_n^k - 1} = \dfrac{1^k - 1}{\omega_n^k - 1} = 0 \]六、一元多项式
1. 概念
由数和字母的积组成的代数式叫做单项式，单独的一个数或一个字母也叫做单项式（例：\(0\) 可看做 \(0\) 乘 \(a\)，\(1\) 可以看做 \(1\) 乘指数为 \(0\) 的字母，\(b\) 可以看做 \(b\) 乘 \(1\)），分数和字母的积的形式也是单项式。
在数学中，由若干个单项式相加组成的代数式叫做多项式（若有减法：减一个数等于加上它的相反数）。多项式中的每个单项式叫做多项式的项，这些单项式中的最高项次数，就是这个多项式的次数。其中多项式中不含字母的项叫做常数项。
以上摘自 百度百科 - 单项式 和 百度百科 - 多项式。
-
单项式次数：所有变数字母的指数之和
-
多项式次数：多项式 \(F\) 的次数为其最高项的次数，记作 \(\operatorname{degree}(F)\)
-
多项式次数界：多项式 \(F\) 的次数界为任意一个严格大于\(F\) 的次数的整数
对于 FFT，我们只研究一元多项式，即只有一个未知数的多项式。
为了方便，对于多项式的项的系数数列编号从 \(0\) 开始。
2. 系数表示
如果用一个 \(n\) 项数列 \(a\) 表示多项式 \(A\) 的每一项数，\(x\) 表示字母部分，则：
\[A(x) = a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \]或者：
\[A(x) = \sum_{i = 0}^{n-1}{a_ix^i} \]这就是多项式的系数表示法。
明显地，计算一个系数表示的多项式的值为 \(O(n)\)。
一元多项式在系数表示下的加法
整式加法，对应次数项加起来即可。
设多项式 \(A(x) = \sum_{i = 0}^{n-1}{a_ix^i}\) 和多项式 \(B(x) = \sum_{i = 0}^{n-1}{b_ix^i}\)。
则：
\[C(x)=A(x)+B(x)=\sum_{i = 0}^{n-1}{a_ix^i}+\sum_{i = 0}^{n-1}{b_ix^i}=\sum_{i = 0}^{n-1}{(a_i + b_i)x^i} \]若我们运用计算机做整式加法，明显地，时间复杂度为 \(O(n)\)。
一元多项式在系数表示下的乘法
整式乘法，用一个多项式的每一项去乘另一个多项式的每一项，将积相加。
设多项式 \(A(x) = \sum_{i = 0}^{n - 1}{a_ix^i}\) 和多项式 \(B(x) = \sum_{i = 0}^{n-1}{b_ix^i}\)。
则：
\[C(x)=A(x)\cdot B(x)=\sum_{i = 0}^{n - 1}{a_ix^i}\cdot \sum_{i = 0}^{n - 1}{b_ix^i}=\sum_{i = 0}^{\operatorname{degree}(C)}{c_ix^i} \]其中 \(c_i = \sum_{j=0}^{i}{a_jb_{i-j}}\)，\(\operatorname{degree}(C)=\operatorname{degree}(A)+\operatorname{degree}(B)=n - 1+n- 1=2n - 2=2(n-1)\)。
\(c\) 为 \(a,b\) 的卷积，记作 \(c=a\otimes b\)。
若我们运用计算机做整式乘法，明显地，时间复杂度为 \(O(n^2)\)，时间复杂度过高，不优秀。
3. 点值表示
平时我们的多项式表示方法一般都是系数表示，因为较为直观，为多项式的概念。
即为单项式之和的形式。
多项式还有一种表示方法，叫做点值表示。
即用至少 \(n\) 个多项式上的点表示。
\[\{(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2),\cdots,(x_n,A(x_n))\} \]求这个多项式的值，被称为插值，可以使用拉格朗日插值公式进行求解。
由于我自己也不会，所以直接给出公式。
\[A(x)=\sum_{i=0}^{n}{A(x_i)\dfrac{{\textstyle \prod_{j\neq i}{(x-x_j)}}}{{\textstyle \prod_{j\neq i}{(x_i-x_j)}}}} \]时间复杂度为 \(O(n^2)\)，不够优秀。
一元多项式在点值表示下的加法
对于每个对应的 \(x\) 的点，直接相加即可
\[C(x_i)=A(x_i)+B(x_i) \]明显地，时间复杂度为 \(O(n)\)。
一元多项式在点值表示下的乘法
对于每个对应的 \(x\) 的点，直接相加即可
\[C(x_i)=A(x_i)B(x_i) \]明显地，时间复杂度为 \(O(n)\)，效率远远高于在系数表示下的乘法。
七、快速傅里叶变换
我们看到，一元多项式在点值表示下的乘法是非常快的，所以我们想把系数表示的多项式转化成点值表示的多项式，这样就可以快速求出两个多项式的卷积了。但是我们转化的过程是很慢的，而快速傅里叶变换（FFT）就是使用单位根优化这个过程的算法。
1. 离散傅里叶变换
离散傅里叶变换（Discrete Fourier Transform，DFT），是这样一个东西：
\[A(x) = \sum_{i = 0}^{n - 1}{a_ix^i} \ \text{在} \ \omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1} \ \text{处的值} \ y_0,y_1,\cdots,y_{n-1} \]\[y_i = A(\omega_n^i) = \sum_{j = 0}^{n - 1}{\omega_n^{ij} a_j} \]记作：\(\boldsymbol{y}=\text{DFT}_n(\boldsymbol{a})\) 或 \(\boldsymbol {y}=\mathcal{F} \boldsymbol{a}\)
看的可能有点迷惑。
简单来说，就是傅里叶非常无聊，把单位根带到一个多项式里面去了...
然后这个东西可以在其他的地方发挥作用。
2. 快速傅里叶变换
离散傅里叶变换并没有降低时间复杂度，霍纳法则（秦九韶算法）求值时间复杂度为 \(O(n^2)\)。
但是我们发现离散傅里叶变换带入了单位根，我们可以使用单位根的性质进行优化。
首先我们有一个 \(n-1\) 次多项式，\(A(x)\)，有系数向量 \(\boldsymbol{a} = [a_0,a_1,a_2,\cdots,a_{n-1}]^T\)。
注意 FFT 的 \(n\) 为 \(2\) 的幂。
我们现在把这个向量拆开，分为偶数项和奇数项。
记为 \(\boldsymbol{a}^{[0]}\) 和 \(\boldsymbol{a}^{[1]}\)。这两个向量所对应的两个多项式我们分别记为 \(A^{[0]}(x)\) 和 \(A^{[1]}(x)\)。
\[\boldsymbol{a}^{[0]} = [a_0,a_2,a_4,\cdots,a_{n-2}]^T \rightarrow A^{[0]}(x) \\ \boldsymbol{a}^{[1]} = [a_1,a_3,a_5,\cdots,a_{n-1}]^T \rightarrow A^{[1]}(x) \]现在把三个多项式都带入 \(x\)。
\[A(x) = a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \\ A^{[0]}(x) = a_0+a_2x+a_4x^2+a_6x^3+\cdots+a_{n-2}x^{\frac{n}{2}-1} \\ A^{[1]}(x) = a_1+a_3x+a_5x^2+a_7x^3+\cdots+a_{n-1}x^{\frac{n}{2}-1} \]把后两个多项式自变量换成 \(x^2\)。
\[A(x) = a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \\ A^{[0]}(x^2) = a_0+a_2x^2+a_4x^4+a_6x^6+\cdots+a_{n-2}x^{n-2} \\ A^{[1]}(x^2) = a_1+a_3x^2+a_5x^4+a_7x^6+\cdots+a_{n-1}x^{n-2} \]把 \(A^{[1]}\) 乘上 \(x\)。
\[A(x) = a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_{n-1}x^{n-1} \\ A^{[0]}(x^2) = a_0+a_2x^2+a_4x^4+a_6x^6+\cdots+a_{n-2}x^{n-2} \\ xA^{[1]}(x^2) = a_1x+a_3x^3+a_5x^5+a_7x^7+\cdots+a_{n-1}x^{n-1} \]然后就会发现：
\[A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2) \]我们现在就将原问题转化成了求次数界为 \(\dfrac{n}{2}\) 的两个多项式在每个单位根的平方上的值然后相加的值。
比较拗口但是还是很好理解的。
我们把两个具体的单位根带入看看。
\[A(\omega_n^k) = A^{[0]}((\omega_n^k)^2) + \omega_n^kA^{[1]}((\omega_n^k)^2) \\ A(\omega_n^{k+\frac{n}{2}}) = A^{[0]}((\omega_n^{k+\frac{n}{2}})^2) + \omega_n^{k+\frac{n}{2}}A^{[1]}((\omega_n^{k+\frac{n}{2}})^2) \]再根据消去引理：\(\omega_{dn}^{dk}=\omega_n^k\)，
和折半引理：\((\omega_n^{k+\frac{n}{2}})^2=(\omega_n^k)^2=\omega_\frac{n}{2}^k\)，
以及折半引理证明中的 \(\omega_n^{k+\frac{n}{2}} = -\omega_n^k\)。
带入化简得到：
\[A(\omega_n^k) = A^{[0]}(\omega_\frac{n}{2}^k) + \omega_n^kA^{[1]}(\omega_\frac{n}{2}^k) \\ A(\omega_n^{k+\frac{n}{2}}) = A^{[0]}(\omega_\frac{n}{2}^k) - \omega_n^kA^{[1]}(\omega_\frac{n}{2}^k) \]发现了吗？除了符号不同之外都是一样的。
所以我们只要知道了前一半，后一半的值我们也能知道了。
这里写一下伪代码。
\[\begin{array}{ll} \\ 1 & \operatorname{FFT}(\boldsymbol{a}, limit) \\ 2 & \qquad \textbf{if } limit = 1 \\ 3 & \qquad \qquad \textbf{return } \\ 4 & \qquad \boldsymbol{a}^{[0]} = (a_0,a_2,a_4,\cdots,a_{n-2}) \\ 5 & \qquad \boldsymbol{a}^{[1]} = (a_1,a_3,a_5,\cdots,a_{n-1}) \\ 6 & \qquad \operatorname{FFT}(\boldsymbol{a}^{[0]}, \frac{limit}{2}) \\ 7 & \qquad \operatorname{FFT}(\boldsymbol{a}^{[1]}, \frac{limit}{2}) \\ 8 & \qquad \omega_n = e^{\frac{2\pi}{n}i} = \cos(\frac{2\pi}{n})+i\sin(\frac{2\pi}{n}) \\ 9 & \qquad \omega = 1 \\ 10 & \qquad \textbf{for } k = 0 \ldots \frac{n}{2} - 1 \\ 11 & \qquad \qquad \boldsymbol{a}_k = \boldsymbol{a}_k^{[0]} + \omega \boldsymbol{a}_k^{[1]} \\ 12 & \qquad \qquad \boldsymbol{a}_{k+\frac{n}{2}} = \boldsymbol{a}_k^{[0]} - \omega \boldsymbol{a}_k^{[1]} \\ 13 & \qquad \qquad \omega = \omega \omega_n \end{array} \]3. 代码实现
例：洛谷 - P3803 【模板】多项式乘法（FFT）。
void FFT(Complex_Number * p, ll len, bool IDFT) {
if(len == 1) return;
for(ll i = 0; i < len; i++) poly[i] = p[i];
Complex_Number * lp = p, * rp = p + (len >> 1);
for(ll i = 0; i < (len >> 1); i++) {
lp[i] = poly[i << 1]; rp[i] = poly[i << 1 | 1];
}
FFT(lp, len >> 1, IDFT);
FFT(rp, len >> 1, IDFT);
Complex_Number unitroot(cos(2.0 * Pi / len), sin(2.0 * Pi / len));
if(IDFT) unitroot.imag *= -1;
Complex_Number w(1.0, 0.0), t;
for(ll i = 0; i < (len >> 1); i++) {
t = w * rp[i];
poly[i] = lp[i] + t;
poly[i + (len >> 1)] = lp[i] - t;
w = w * unitroot;
}
for(ll i = 0; i < len; i++) p[i] = poly[i];
}
int main(){
scanf("%lld%lld",&n,&m);
for(ll i = 0; i <= n; i++) scanf("%lf", &poly1[i].real);
for(ll i = 0; i <= m; i++) scanf("%lf", &poly2[i].real);
for(m += n, n = 1; n <= m; n <<= 1);
FFT(poly1, n, false); FFT(poly2, n, false);
for(ll i = 0; i < n; i++) poly1[i] = poly1[i] * poly2[i];
FFT(poly1, n, true);
for(ll i = 0; i <= m; i++) printf("%lld ", ll(poly1[i].real / n + 0.45));
return 0;
}
八、大数乘法
1. 基本思路
把乘数分解为多项式即 \(a_1+a_2x+a_3x^2+a_4x^3+\cdots\)，其中 \(x = 10\) 的形式。
\[A = a_1+a_2x+a_3x^2+a_4x^3+\cdots \\ B = b_1+b_2x+b_3x^2+b_4x^3+\cdots \]其中 \(x = 10\)。
然后使用 FFT 相乘，最后全加起来求值即可。
2. 代码实现
例：洛谷 - P1919 【模板】A*B Problem升级版（FFT快速傅里叶），洛谷 - P1303 A*B Problem。
ll mul(ll number1, ll number2) {
ll n, m;
fill(poly1, poly1 + N, null);
fill(poly2, poly2 + N, null);
fill(poly, poly + N, null);
string num1 = to_string(number1), num2 = to_string(number2), ans = "";
n = num1.length() - 1; m = num2.length() - 1;
for(ll i = 0; i <= n; i++) poly1[i].real = num1[n - i] - '0';
for(ll i = 0; i <= m; i++) poly2[i].real = num2[m - i] - '0';
for(m += n, n = 1; n <= m; n <<= 1);
FFT(poly1, n, false); FFT(poly2, n, false);
for(ll i = 0; i < n; i++) poly1[i] = poly1[i] * poly2[i];
FFT(poly1, n, true);
ll t = 0;
for(ll i = 0; i <= m; i++) {
ll now = ll(poly1[i].real / n + 0.49);
ans = char(add((now + t) % 10, 48)) + ans;
t = (now + t) / 10;
}
ans = char(add((t % 10), 48)) + ans;
bool firstzero = true;
string ret = "";
for(ll i = 0; i < ans.length(); i++) {
if(ans[i] == '0' && firstzero) continue;
if(firstzero)firstzero = false;
ret += ans[i];
}
if(firstzero) ret += "0";
return atoll(ret.c_str());
}
九、线段树
1. 概念
线段树（Segment Tree），数据结构的一种，可以快速维护某些的序列的修改和查询。
线段树是一颗二叉树，其运用了分块的思想，对于每一个节点的点权是其代表区间之和，从总区间二分下去扩展到每一个项。
线段树由于是平衡的二叉树，树高 \(\log{n}\)，所以可以做到在 \(O(\log{n})\) 的时间内修改节点或求和。
线段树不止可以求 RSQ，还可以求 RMQ 等等奇奇怪怪的东西。
2. 基本思路
在本题中，我们只要求出区间和即可。
要查找一个点，我们可以从根节点走下去，找到一个包含区间即直接返回，直到找到所有的区间，可从定义得证。
要修改一个点，我们可以考虑更普遍的修改区间。若要修改区间，如果一个一个修改，则会导致时间复杂度变高，我们可以使用一个叫做懒标记（lazy tag）的东西，根据线段树的性质，我们可以先得到它的和返回，然后标记，对于有查询到的时候扩展即可。
注意每次都要保证线段树的性质，维护线段树。
如何获得加上一个数之和的区间和？
根据结合律 \((a_1 + x) + (a_2 + x) + \cdots + (a_n + x) = (a_1 + \cdots + a_n) + nx\)。
即可得到区间和。
3. 建树
我们可以使用数组法建造二叉树，如果用其他方法也可以。
在这里，我们，就使用数组法，因为方便。
我习惯用一个叫 Tree 的数组存树，接下来的树上节点，都是 Tree。
则节点 \(p\) 的左孩子：p << 1；
右孩子：p << 1 | 1。
对于每一个节点都要有它的 \(l\) 和 \(r\)（区间左边和区间右边，包含），而这个东西 在本题中 我们递归的时候每一次就可以根据上一个算出，所以不需要记录。
二分递归建树，取 \(mid = \dfrac{l + r}{2}\)，然后两边递归建树。
当 \(l = r\) 时，即为叶子节点，停止递归，让这个点的值为第 \(l\) 或者 \(r\) 项数，因为是区间范围。
可以证明，建树最后时不会漏下一个任何一个点。
建树的时候，要维护线段树。
即根据定义，除叶子节点外的其他点的权为其左右个儿子之和。
void build(ll p, ll l, ll r) {
if(l == r) {
Tree[p] = a[l];
return;
}
ll mid = (l + r) >> 1;
build(mul(p, 2), l, mid);
build(add(mul(p, 2), 1), add(mid, 1), r);
Tree[p] = Tree[mul(p, 2)] + Tree[add(mul(p, 2), 1)];
}
4. lazy tag 的下放
void push_down(ll p, ll l, ll r) {
ll mid = (l + r) >> 1;
tag[mul(p, 2)] += tag[p];
tag[add(mul(p, 2), 1)] += tag[p];
Tree[mul(p, 2)] += tag[p] * add(mid - l, 1);
Tree[add(mul(p, 2), 1)] += tag[p] * (r - mid);
tag[p] = 0;
}
5. 修改
void update(ll p, ll l, ll r, ll L, ll R, ll val) {
if(r < L || R < l) return;
if(L <= l && r <= R) {
tag[p] += val;
Tree[p] += val * add(r - l, 1);
return;
}
ll mid = (l + r) >> 1;
push_down(p, l, r);
if(R <= mid) update(mul(p, 2), l, mid, L, R, val);
else if(mid < L) update(add(mul(p, 2), 1), add(mid, 1), r, L, R, val);
else {
update(mul(p, 2), l, mid, L, R, val);
update(add(mul(p, 2), 1), add(mid, 1), r, L, R, val);
}
Tree[p] = Tree[mul(p, 2)] + Tree[add(mul(p, 2), 1)];
}
6. 查询
ll query(ll p, ll l, ll r, ll L, ll R) {
if(r < L || R < l) return 0;
if(L <= l && r <= R) return Tree[p];
ll mid = (l + r) >> 1;
push_down(p, l, r);
if(R <= mid) return query(mul(p, 2), l, mid, L, R);
else if(mid < L) return query(add(mul(p, 2), 1), add(mid, 1), r, L, R);
else {
return query(mul(p, 2), l, mid, L, R) + query(add(mul(p, 2), 1), add(mid, 1), r, L, R);
}
}
7. 代码实现
结合上述操作，即可得到最终代码。
例：洛谷 - P3372 【模板】线段树 1。
ll a[N] = {0}, Tree[N << 2], tag[N << 2];
void build(ll p, ll l, ll r) {
if(l == r) {
Tree[p] = a[l];
return;
}
ll mid = (l + r) >> 1;
build(mul(p, 2), l, mid);
build(add(mul(p, 2), 1), add(mid, 1), r);
Tree[p] = Tree[mul(p, 2)] + Tree[add(mul(p, 2), 1)];
}
void push_down(ll p, ll l, ll r) {
ll mid = (l + r) >> 1;
tag[mul(p, 2)] += tag[p];
tag[add(mul(p, 2), 1)] += tag[p];
Tree[mul(p, 2)] += tag[p] * add(mid - l, 1);
Tree[add(mul(p, 2), 1)] += tag[p] * (r - mid);
tag[p] = 0;
}
void update(ll p, ll l, ll r, ll L, ll R, ll val) {
if(r < L || R < l) return;
if(L <= l && r <= R) {
tag[p] += val;
Tree[p] += val * add(r - l, 1);
return;
}
ll mid = (l + r) >> 1;
push_down(p, l, r);
if(R <= mid) update(mul(p, 2), l, mid, L, R, val);
else if(mid < L) update(add(mul(p, 2), 1), add(mid, 1), r, L, R, val);
else {
update(mul(p, 2), l, mid, L, R, val);
update(add(mul(p, 2), 1), add(mid, 1), r, L, R, val);
}
Tree[p] = Tree[mul(p, 2)] + Tree[add(mul(p, 2), 1)];
}
ll query(ll p, ll l, ll r, ll L, ll R) {
if(r < L || R < l) return 0;
if(L <= l && r <= R) return Tree[p];
ll mid = (l + r) >> 1;
push_down(p, l, r);
if(R <= mid) return query(mul(p, 2), l, mid, L, R);
else if(mid < L) return query(add(mul(p, 2), 1), add(mid, 1), r, L, R);
else {
return query(mul(p, 2), l, mid, L, R) + query(add(mul(p, 2), 1), add(mid, 1), r, L, R);
}
}
十、完整代码
结合以上知识，应该对于如何求解 A + B Problem 有所了解。
现在我们重新在梳理一遍思路。
我们先建线段树，但是由于初始值都为 \(0\)，于是我们可以省略 build（建树）这一操作。
然后我们向线段树的位置 \(1\) 加上 \(a\)，位置 \(2\) 加上 \(b\)，最后求区间和即可。
然后我们在写线段树时应该注意使用高精度加法，如前面说的。
在本题的线段树中，乘法不会出现负数，其实如果出现负数也非常好除了。
我们可以直接先除里符号位。
由于，乘法不会出现负数，所以我们可以直接套用之前的 FFT 大数乘法模板。
大数乘法模板也有用到加法，注意使用高精度加法。
然后全部结合起来，就完成了本题的程序。
现在给出 A + B Problem 的完整代码。
题目链接：洛谷 - P1001 A+B Problem。
/*Copyright (C) 2013-2022 LZE*/
#include
#define fo(x) freopen(#x".in", "r", stdin); freopen(#x".out", "w", stdout);
#define INF 0x7fffffff
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
const int N = 1010;
const int M = 1010;
const double Pi = acos(-1.0);
struct Complex_Number {
double real, imag;
Complex_Number() {}
Complex_Number(double num1, double num2) : real(num1) , imag(num2) {}
} poly1[N], poly2[N], poly[N], null;
Complex_Number operator + (const Complex_Number &num1, const Complex_Number &num2) {
return Complex_Number(num1.real + num2.real, num1.imag + num2.imag);
}
Complex_Number operator - (const Complex_Number &num1, const Complex_Number &num2) {
return Complex_Number(num1.real - num2.real, num1.imag - num2.imag);
}
Complex_Number operator * (const Complex_Number &num1, const Complex_Number &num2) {
return Complex_Number(num1.real * num2.real - num1.imag * num2.imag, num1.real * num2.imag + num1.imag * num2.real);
}
void FFT(Complex_Number * p, ll len, bool IDFT) {
if(len == 1) return;
for(ll i = 0; i < len; i++) poly[i] = p[i];
Complex_Number * lp = p, * rp = p + (len >> 1);
for(ll i = 0; i < (len >> 1); i++) {
lp[i] = poly[i << 1]; rp[i] = poly[i << 1 | 1];
}
FFT(lp, len >> 1, IDFT);
FFT(rp, len >> 1, IDFT);
Complex_Number unitroot(cos(2.0 * Pi / len), sin(2.0 * Pi / len));
if(IDFT) unitroot.imag *= -1;
Complex_Number w(1.0, 0.0), t;
for(ll i = 0; i < (len >> 1); i++) {
t = w * rp[i];
poly[i] = lp[i] + t;
poly[i + (len >> 1)] = lp[i] - t;
w = w * unitroot;
}
for(ll i = 0; i < len; i++) p[i] = poly[i];
}
ll add(ll number_1, ll number_2) {
string n1_number = to_string(number_1), n2_number = to_string(number_2);
ll length_of_n1 = n1_number.length();
ll length_of_n2 = n2_number.length();
ll max_of_string_length = max(length_of_n1, length_of_n2);
if(length_of_n1 < max_of_string_length)
for(ll i = length_of_n1; i < max_of_string_length; i++)
n1_number = "0" + n1_number;
else for(ll i = length_of_n2; i < max_of_string_length; i++)
n2_number = "0" + n2_number;
string add_sum = "";
ll last_num = 0, add_sum_1 = 0, add_sum_2 = 0, add_sum_add = 0;
for(ll i = max_of_string_length - 1;i >= 0; i--) {
add_sum_1 = n1_number[i] - '0';
add_sum_2 = n2_number[i] - '0';
add_sum_add = add_sum_1 + add_sum_2 + last_num;
add_sum = to_string(add_sum_add % 10) + add_sum;
last_num = add_sum_add / 10;
}
if(last_num) add_sum = to_string(last_num) + add_sum;
return atoll(add_sum.c_str());
}
ll mul(ll number1, ll number2) {
ll n, m;
fill(poly1, poly1 + N, null);
fill(poly2, poly2 + N, null);
fill(poly, poly + N, null);
string num1 = to_string(number1), num2 = to_string(number2), ans = "";
n = num1.length() - 1; m = num2.length() - 1;
for(ll i = 0; i <= n; i++) poly1[i].real = num1[n - i] - '0';
for(ll i = 0; i <= m; i++) poly2[i].real = num2[m - i] - '0';
for(m += n, n = 1; n <= m; n <<= 1);
FFT(poly1, n, false); FFT(poly2, n, false);
for(ll i = 0; i < n; i++) poly1[i] = poly1[i] * poly2[i];
FFT(poly1, n, true);
ll t = 0;
for(ll i = 0; i <= m; i++) {
ll now = ll(poly1[i].real / n + 0.49);
ans = char(add((now + t) % 10, 48)) + ans;
t = (now + t) / 10;
}
ans = char(add((t % 10), 48)) + ans;
bool firstzero = true;
string ret = "";
for(ll i = 0; i < ans.length(); i++) {
if(ans[i] == '0' && firstzero) continue;
if(firstzero)firstzero = false;
ret += ans[i];
}
if(firstzero) ret += "0";
return atoll(ret.c_str());
}
ll a[N] = {0}, Tree[N << 2], tag[N << 2];
void build(ll p, ll l, ll r) {
if(l == r) {
Tree[p] = a[l];
return;
}
ll mid = (l + r) >> 1;
build(mul(p, 2), l, mid);
build(add(mul(p, 2), 1), add(mid, 1), r);
Tree[p] = Tree[mul(p, 2)] + Tree[add(mul(p, 2), 1)];
}
void push_down(ll p, ll l, ll r) {
ll mid = (l + r) >> 1;
tag[mul(p, 2)] += tag[p];
tag[add(mul(p, 2), 1)] += tag[p];
Tree[mul(p, 2)] += tag[p] * add(mid - l, 1);
Tree[add(mul(p, 2), 1)] += tag[p] * (r - mid);
tag[p] = 0;
}
void update(ll p, ll l, ll r, ll L, ll R, ll val) {
if(r < L || R < l) return;
if(L <= l && r <= R) {
tag[p] += val;
Tree[p] += val * add(r - l, 1);
return;
}
ll mid = (l + r) >> 1;
push_down(p, l, r);
if(R <= mid) update(mul(p, 2), l, mid, L, R, val);
else if(mid < L) update(add(mul(p, 2), 1), add(mid, 1), r, L, R, val);
else {
update(mul(p, 2), l, mid, L, R, val);
update(add(mul(p, 2), 1), add(mid, 1), r, L, R, val);
}
Tree[p] = Tree[mul(p, 2)] + Tree[add(mul(p, 2), 1)];
}
ll query(ll p, ll l, ll r, ll L, ll R) {
if(r < L || R < l) return 0;
if(L <= l && r <= R) return Tree[p];
ll mid = (l + r) >> 1;
push_down(p, l, r);
if(R <= mid) return query(mul(p, 2), l, mid, L, R);
else if(mid < L) return query(add(mul(p, 2), 1), add(mid, 1), r, L, R);
else {
return query(mul(p, 2), l, mid, L, R) + query(add(mul(p, 2), 1), add(mid, 1), r, L, R);
}
}
int main() {
null = Complex_Number(0, 0);
ll n1, n2; cin >> n1 >> n2;
build(1, 1, 2);
update(1, 1, 2, 1, 1, n1); update(1, 1, 2, 2, 2, n2);
cout << query(1, 1, 2, 1, 2) << endl;
return 0;
}
好的，恭喜你，你先已经学会了 A + B Problem 了。
我们要学会举一反三，这里给出一道也有用到 A + B 的题目，希望你也可以顺利地 AC 它！
练习题：洛谷 - P6474 [NOI Online #2 入门组] 荆轲刺秦王。