用于计算一类朴素卷积问题。
笔者学习的是这份博客
内容中可能有很多相同之处,敬请谅解。
现在要计算两个一元\(n\)次多项式\(F(x)\)与\(G(x)\)的乘积,如何计算?
前置知识:多项式的表示方法
一. 系数表示法
对于一个\(n\)次多项式\(F(x)\),它可以被表示成
\[F(x) = a_nx^n+a_{n-1}x^{n-1}+...+a_1x^1+a_0x^0.\]
更加形式化的来说,它可以表示成
\[F(x) = \sum_{i=0}^{n} a_ix^i.\]
举个例子,2次多项式,其中\(a_0=1,a_1=2,a_2=1\)
那么\(F(x)=1x^2+2x+1.\)这样即可通俗的表示出一个\(n\)次多项式。
二. 点值表示法众所周知,两个点确定一个一次函数,三个点确定一个二次函数。
所以,\(n+1\)个点确定一个一元\(n\)次多项式。
所以我们可以通过\(n+1\)个点来表示它。
那么相乘之后的点值如何计算?
比如说两个2次多项式\(F(x)=x^2+2x+1\)(红色)与\(G(x)=3x^2-4x-2\)(蓝色),它们的图像如图所示:
那么观察图中\(x=1\)时的情况。此时\(F(1)=4,G(1)=-3.\)
所以,显然,\(F(1)\times G(1)=-12\).也就是说在\(Z=F*G\)这一多项式内,带入\(1\),得到的结果是\(-12\).
等等,好像有哪里不对。如果说\(Z=F*G\)的话,那么Z的次数应该是\(2n\).
那\(Z\)需要\(2n+1\)个点来确定。但是原来只需要\(n+1\)个点,咋办?
很简单,在原来的多项式里每个都多加\(n\)个点即可。反正多项式已知。
这样就可以用点值来进行操作。也就是说先转成点值,再一乘,再转回来,就是计算流程。
但是好像还是很慢。那么如何优化呢?
复数部分
复数,即形如\(a+bi\)的数,其中\(i^2=-1.\) \(a\)称为实部,\(bi\)称为虚部。
或者说:在一个数轴上(只有x轴),我们可以表示出任何实数。
那么,多加一维(y轴),也就是类似于平面直角坐标系一样,我们就可以表示出任意一个复数。
所以我们把这个坐标系叫做复平面,其中x轴称为实轴,y轴称为虚轴。
复数运算
复数相加:实部相加,虚部相加,例如
\[(a+bi)+(c+di)=(a+c)+(b+d)i.\]
复数相减:同理。
\[(a+bi)-(c+di)=(a-c)+(b-d)i.\]
复数相乘:像一次多项式一样相乘。 注意\(i^2=-1\).
\[(a+bi)(c+di)=ac+(ad+bc)i-bd=(ac-bd)+(ad+bc)i.\]
复数相除:
相信大家都学过共轭根式。同样的,复数也有共轭。
即:\(a+bi\)的共轭为\(a-bi\)。
这两个复数卡乘在一起一定是个实数。即
\[(a+bi)(a-bi)=a^2-(bi)^2=a^2+b^2.\]
所以再除的时候,将分子分母同乘分母的共轭,就可以将分母有理化。
即
\[\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{c^2+d^2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i.\]
复数逆元:
\[\frac{1}{a+bi}=\frac{a}{a^2+b^2}-\frac{b}{a^2+b^2}i.\]
同样的,复数运算在复平面内也有一定规律可循。
考虑在复平面内的两个复数:(借张图)
表示的是复数\((1+4i)\)与复数\((3+2i)\)相乘所得的结果:\((-5+14i)\)。
设其中\((5,0)\)点为位置\(P\),则\(\angle{POC}=\angle{BOA}.\)
还有:\(\overline{OB}\times\overline{OC}=\overline{OA}.\)
第二个证明:勾股定理。先把\(i\)消掉。
\(\overline{OB}^2=a^2+b^2.\) \(\overline{OC}^2=c^2+d^2.\)
\(\overline{OB}^2\times \overline{OC}^2=(a^2+b^2)(c^2+d^2)=a^2c^2+a^2d^2+b^2c^2+b^2d^2.\)
\(\overline{OA}^2=(ac-bd)^2+(ad+bc)^2=a^2c^2+a^2d^2+b^2c^2+b^2d^2.\)得证。
我们将复数中,复数向量的长度称为模长,向量与x轴正方向的夹角称为幅角。
根据上面的东西:复数相乘时,模长相乘,幅角相加。
单位根
一个n次的单位根即为方程\(x^n=1\)的复数解。
考虑这样一个图:
其中圆上的所有复数模长都是1,这个圆称为单位圆。
考虑\(|x|\)的取值范围:
如果\(|x|<1\),那么\(|x^n|<1\)(模长*模长,越乘越短)
如果\(|x|>1\),那么\(|x^n|>1\)(模长*模长,越乘越长)
所以一定有\(|x|=1\),即点一定在这个单位圆上。
那么还有一条,就是一个n次的单位根共有n个,并且从(1,0)开始,每次逆时针走\(\frac{1}{n}\)个圆周。分别称为\(\omega_{n}^{0},\omega_{n}^{1},...\omega_{n}^{n-1}.\)(逆时针编号!!!)
其实还有\(k\geq n\)和\(k<0\)的单位根,本质上就模个n就可以了,这里不再多加考虑。
有关单位根也有很多性质:
- \(\omega_{n}^{0}=1\)
- \(\omega_{n}^{k}=(\omega_{n}^{1})^k\),很好理解,就是每次转\(\frac{1}{n}\)圆周。
- \(\omega_{n}^{j} \times \omega_{n}^{k}=\omega_{n}^{j+k}.\)
证明:很简单,代入性质2即可。
根据3,有:
\[\omega_{n}^{j+k}=\omega_{n}^{j} \times \omega_{n}^{k}.\]
\[(\omega_{n}^{k})^{-1}=\frac{1}{\omega_{n}^{k}}=\frac{\omega_{n}^{0}}{\omega_{n}^{k}}=\omega_{n}^{-k}=\omega_{n}^{n-k}.\] (根据上面模个n的说法)
\[\omega_{pn}^{pk}=\omega_{n}^{k}.\](将蛋糕分成pn份取pk份=将蛋糕分成n份取k份)
\[\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}.\](相当于绕了半圈\(\rightarrow\)取了个相反数)
\[(\omega_{n}^{k})^j=(\omega_{n}^{j})^k.\](大小为k,共有j个=大小为j,共有k个)
...差不多就这些了
---
\(FFT\)流程:\(DFT,IDFT\)
\(DFT\)加速版本:分治
首先,你需要控制得当将整个n-1次多项式(保证n是2的正整数次幂)按系数奇偶拍成两部分:
\[F(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+...+a_{n-1}x^{n-1}).\]
再设两个n/2项多项式:
\[FL(x)=a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1}\]
\[FR(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}\]
那么有\(F(x)=FL(x^2)+x*FR(x^2).\)(不懂可以自己设一下代入进去)
下面,我们来代入单位根,看看会有什么神奇的事情发生。
比如说我们代一个\(\omega_{n}^{k}(k<\frac{n}{2}):\)
\[F(\omega_{n}^{k})=FL((\omega_{n}^{k})^2)+\omega_{n}^{k}\times FR((\omega_{n}^{k})^2)\]
其中,\((\omega_{n}^{k})^2=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k}.\)然后代入:
\[F(\omega_{n}^{k})=FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}\times FR(\omega_{\frac{n}{2}}^{k}).\]
然后呢?有啥用?
不难发现,如果我们知道了\(FL(x),FR(x)\)在x取\(\omega_{\frac{n}{2}}^{0},\omega_{\frac{n}{2}}^{1}...\omega_{\frac{n}{2}}^{\frac{n}{2}-1}\)下的所有值,我们就可以知道\(F(x)\)在x取\(\omega_{n}^{0},\omega_{n}^{1}...\omega_{n}^{\frac{n}{2}-1}\)下的所有值(代回原式)。
但是这还不是全部啊,剩下那个部分怎么办?
没关系,考虑\(k<\frac{n}{2}\),代入\(\omega_{n}^{k+\frac{n}{2}}\)时的情况。
\[F(\omega_{n}^{k+\frac{n}{2}})=FL((\omega_{n}^{k+\frac{n}{2}})^2)+\omega_{n}^{k+\frac{n}{2}}\times FR((\omega_{n}^{k+\frac{n}{2}})^2).\]
其中,\((\omega_{n}^{k+\frac{n}{2}})^2=\omega_{n}^{2k+n}.\)
故
\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{n}^{2k+n}).\]
\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{n}^{2k})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{n}^{2k}).\]
(模性质)
\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{\frac{n}{2}}^{k}).\]
(性质3 延展)
\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}\times FR(\omega_{\frac{n}{2}}^{k}).\]
(取反性质)
然后观察这个式子和前面式子的区别,发现就是把+号改成了-号。也就是说我们可以通过那组数据来算出当\(k>\frac{n}{2}\)时的情况。这就是加速DFT的原理。
IDFT
至此,我们已经可以将其转化为点值,那么如何还原回来呢?使用IDFT插值。
考虑一个向量\(c_0,c_1...c_{n-1}\),满足\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i\times y_i.\]
其中\(y_i\)表示乘出来的点的y坐标。
那么,这个式子也就代表一个以\(y\)为系数的方程在取\(x=\omega_{n}^{-k}\)时的复数解。
我们尝试着来化简一波。
\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i\times y_i\]
\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i(\sum_{j=0}^{n-1}a_j(\omega_{n}^{i})^j)\]
\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j})^i)\]
\[c_k=\sum _{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j})^i(\omega_{n}^{-k})^i)\]
\[c_k=\sum _{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j-k})^i)\]
\[c_k=\sum _{j=0}^{n-1}a_j\color{red}{(\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i)}\]
发现后面标红的东西就是一个等比数列,那么我们可以直接计算。即
\[\sum_{i=0}^{n-1}X^i=\frac{X^n-1}{X-1}.\]
我们设这个东西为公式\(Z.\)
于是就有
\[c_k=\sum _{j=0}^{n-1}a_jZ(\omega_{n}^{j-k}).\]
观察函数\(Z\)的相关性质。
赫然发现:我们带进去的是单位根!单位根啥性质来着?\(X^n=1!\)
也就是说当\(X≠1\),也就是\(k≠0\)时,原式\(=0.\) (分子为0)
那么\(X=1\)时,也就是\(k=0\)时呢?这时候这个求和公式就不管用了,带回原式算,发现由于\(1^i=1\)(i是自然数),那么\(Z(x)=n.\)
所以说,什么情况下带进去的值\(k=0\)呢?显然是\(j-k\)的时候!此时单项\(=na_k\).
否则的话会乘一个0,答案固然也是0.
所以有:\[c_k=na_k\](飞跃性的一步!!!)
回顾整个式子,发现我们算\(c\)的时候连\(a\)的影子都没有,现在我们竟然推出了一个和\(a\)有关的关系式!太神奇了!这就是\(IDFT\)的原理。
至此,\(\text{FFT}\)的重要部分已经全部介绍完毕,下面进入代码实现与细节处理部分。
代码实现
朴素DFT:
complex < double > omega(const int n , const int k) { // 单位根
complex < double > ret = {cos(2 * Pi / n * k) , sin(2 * Pi / n * k)};
if(!inv) return ret;
else return conj(ret);
}
void DFT(complex < double > a , const int n) {
if(n == 1) return ;
static complex < double > buf[MAXN];
const int m = n / 2;
for(int i = 0; i < m; i++) {
buf[i] = a[(i << 1)];
buf[i + m] = a[(i << 1 | 1)];
}
copy(buf , buf + n , a);
complex < double > *a1 = a , *a2 = a + m;
DFT(a1 , m); DFT(a2 , m);
for(int i = 0; i < m; i++) {
complex < double > w = omega(n , i);
buf[i] = a1[i] + w * a2[i];
buf[i + m] = a1[i] - w * a2[i];
}
copy(buf , buf + n , a);
}
发现这样的效率显然不行,那么考虑迭代形式。
首先观察分治到最后的二进制位表示:
000 001 010 011 100 101 110 111
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
然后...我们惊奇的发现最后的二进制位就是之前的二进制位倒过来!TQL!
接着考虑空间优化。
不难发现,由于DFT的两个式子十分相近,所以我们可以只算一次,同时通过技巧性的操作直接覆盖原数组,而不再新开一个。这个操作被称为蝴蝶操作。
蝴蝶操作的代码:
for(int l = 2; l <= n; l <<= 1) {
int m = l / 2;
for(complex < double > *p = a; p != a + n; p += l) {
for(int i = 0; i < m; i++) {
complex < double > t = omega[n / l * i] * p[m + i];
p[m + i] = p[i] - t;
p[i] += t;
}
}
}
然后穿起来,与以前的FFT模板结合即可。
#include <bits/stdc++.h>
#define dbg(x) cerr << #x " = " << x << "\n"
#define INF 0x3f3f3f3f
/* do not give up ! */
/* try your best! */
/* Read the meaning clearly! */
typedef long long LL;
typedef long double ld;
typedef unsigned long long ULL;
using namespace std;
template < typename T > inline void inp(T& t) {
char c = getchar(); T x = 1; t = 0; while(!isdigit(c)) {if(c == '-') x = -1; c = getchar();}
while(isdigit(c)) t = t * 10 + c - '0' , c = getchar();t *= x;
}
template < typename T , typename... Args > inline void inp(T& t , Args&... args) {inp(t); inp(args...);}
template < typename T > inline void outp(T t) {
if(t < 0) putchar('-') , t = -t; T y = 10 , len = 1;
while(y <= t) y *= 10 , len++; while(len--) y /= 10 , putchar(t / y + 48) , t %= y;
}
template < typename T > inline T mul(T x , T y , T MOD) {x=x%MOD,y=y%MOD;return ((x*y-(T)(((ld)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;}
const int MAXN = 2097152 + 10;
const double Pi = acos(-1.0);
int ans[MAXN];
struct FastFourierTransform {
complex < double > omega[MAXN] , omegaInverse[MAXN];
void init(const int n) {
for(int i = 0; i < n; i++) {
omega[i] = complex < double > (cos(2 * Pi / n * i) , sin(2 * Pi / n * i));
omegaInverse[i] = conj(omega[i]);
}
}
void Transform(complex < double > *a , const int n , const complex < double > *omega) {
int k = 0;
while((1 << k) < n) ++k;
for(int i = 0; i < n; i++) {
int t = 0;
for(int j = 0; j < k; j++) if(i & (1 << j)) t |= (1 << (k - j - 1));
if(i < t) swap(a[i] , a[t]);
}
for(int l = 2; l <= n; l <<= 1) {
int m = l / 2;
for(complex < double > *p = a; p != a + n; p += l) {
for(int i = 0; i < m; i++) {
complex < double > t = omega[n / l * i] * p[m + i];
p[m + i] = p[i] - t;
p[i] += t;
}
}
}
}
void DFT(complex < double > *a , const int n) {
Transform(a , n , omega);
}
void IDFT(complex < double > *a , const int n) {
Transform(a , n , omegaInverse);
for(int i = 0; i < n; i++) a[i] /= n;
}
}FFT;
int main() {
int nn , mm , x; inp(nn , mm);
nn++ , mm++;
int n = 1; while(n < nn + mm) n <<= 1;
complex < double > a[MAXN] , b[MAXN];
for(int i = 0; i < nn; i++) {
inp(x); a[i].real(x);
}
for(int i = 0; i < mm; i++) {
inp(x); b[i].real(x);
}
FFT.init(n);
FFT.DFT(a , n);
FFT.DFT(b , n);
for(int i = 0; i < n; i++) a[i] *= b[i];
FFT.IDFT(a , n);
for(int i = 0; i < nn + mm - 1; i++) {
ans[i] = static_cast < int > (floor(a[i].real() + 0.5));
cout << ans[i] << " ";
}
return 0;
}
至此,该代码已经可以效率很高地通过P3803一题!完结撒花!
FFT例题将会另外放在一个帖子~