多项式
假设有
n
n
n次多项式:
∑
i
=
0
n
a
i
x
i
\sum_{i=0}^na_ix^i
i=0∑naixi
如何表示这个多项式.
首先想到的是用所有系数表示,即:
{
a
n
,
a
n
−
1
,
a
n
−
2
,
.
.
.
,
a
2
,
a
1
,
a
0
}
\{a_n,a_{n-1},a_{n-2},...,a_2,a_1,a_0\}
{an,an−1,an−2,...,a2,a1,a0}
我们将这样的表示方法为系数表示法.
但是,系数表示法在表示多项式相加和相乘很不方便.
譬如又有
n
n
n次表达式:
∑
i
=
0
n
b
i
x
i
\sum_{i=0}^nb_ix^i
i=0∑nbixi
那么将两式子相乘的时间复杂度是
O
(
n
2
)
O(n^2)
O(n2)的.
如果
∣
a
∣
|a|
∣a∣=
1
0
5
10^5
105,这种算法就不是那么优越.
所幸,我们还有另一种表示方法——点值表示法.
设想,如果我们用足够多的值对(x,y) 表示这个函数:
y
=
∑
i
=
0
n
a
i
x
i
y=\sum_{i=0}^na_ix^i
y=i=0∑naixi
表达式的系数是可以确定的.
更美妙的是,在点值表达式下,多项式加法和乘法是可以在
O
(
n
)
O(n)
O(n)实现的!
比如,有点对
(
a
1
,
X
1
)
(
a
2
,
X
2
)
.
.
.
(
a
n
,
X
n
)
(a_1,X_1)(a_2,X_2)...(a_n,X_n)
(a1,X1)(a2,X2)...(an,Xn)
(
a
1
,
Y
1
)
(
a
2
,
Y
2
)
.
.
.
(
a
n
,
Y
n
)
(a_1,Y_1)(a_2,Y_2)...(a_n,Y_n)
(a1,Y1)(a2,Y2)...(an,Yn)
那么相乘得
(
a
1
,
X
1
Y
1
)
(
a
2
,
X
2
Y
2
)
.
.
.
(
a
n
,
X
n
Y
n
)
(a_1,X_1Y_1)(a_2,X_2Y_2)...(a_n,X_nY_n)
(a1,X1Y1)(a2,X2Y2)...(an,XnYn)
妙啊
转换
但是,我们在实际生活中并不经常使用点值表示法,而是经常使用系数表示法.
于是,便有了两者之间的转换.
系数表示法和点值表示法可以用如下矩阵乘法转换:
[
1
x
1
x
1
2
.
.
.
x
1
n
1
x
2
x
2
2
.
.
.
x
2
n
.
.
.
.
.
1
x
n
x
n
2
.
.
.
x
n
n
]
[
a
0
a
1
.
.
.
a
n
]
=
[
X
1
X
2
.
.
.
X
n
]
\left[\begin{array}{cc} 1\ \ x_1\ \ x_1^2\ ...\ x_1^n\\ \\ 1\ \ x_2\ \ x_2^2\ ...\ x_2^n\\ \\ .\ \ \ .\ \ \ .\ \ \ .\ \ \ .\\ \\ 1\ \ x_n\ \ x_n^2\ ...\ x_n^n\\ \end{array}\right] \left[\begin{array}{cc} a_0\\ \\ a_1\\ \\ ...\\ \\ a_n \end{array}\right]= \left[\begin{array}{cc} X_1\\ \\ X_2\\ \\ ...\\ \\ X_n \\ \end{array}\right]
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡1 x1 x12 ... x1n1 x2 x22 ... x2n. . . . .1 xn xn2 ... xnn⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡a0a1...an⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡X1X2...Xn⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
这种变换,叫作离散傅里叶变换(DFT).
反之,由点值表示法退回系数表示法,叫作离散傅里叶逆变换(IDFT).
花了好久点个赞呗
但是发现,这样转换时间复杂度还是
O
(
n
2
)
O(n^2)
O(n2)的!
所以,就有了快速傅里叶变换(FFT).
快速傅里叶变换
又称(Fast Fourier Transform)(FFT),详情见Mr.Du.
快速傅里叶变换是离散傅里叶变换的快速版本.这还用说吗
它可以以优秀的
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)解决这个问题.
举个例子:
我们代入
a
=
1
a=1
a=1这个特殊值,那么可以得到
X
=
∑
i
=
0
n
a
i
X=\sum_{i=0}^na_i
X=∑i=0nai.
我们代入
a
=
−
1
a=-1
a=−1、
a
=
0
a=0
a=0也可以得到特殊值.
但是,这样的特殊值总是有限的.
所以,我们将眼光转向别处——复数域.
比如说有一个复数坐标系:
这个圆叫作单位圆没错我们见过,上面的点都可以乘方到1.
如果想要乘方
n
n
n次,那么这些值就是这样分布的:
这是
n
=
8
n=8
n=8的情况.最好画了不是吗
//下面所有
n
n
n都是
2
2
2的幂,注意!
我们称这样的数列为
{
ω
1
,
ω
2
,
.
.
.
,
ω
n
}
\{\omega_1,\omega_2,...,\omega_n\}
{ω1,ω2,...,ωn}.(当然是乘方n次后为1的 )
并称这些为
1
1
1的
n
n
n次单位根 .
铺垫
由神奇而美妙的欧拉公式:
e i x = cos x + i sin x e^{ix}=\cos x+i\sin x eix=cosx+isinx
还有复数的极角表达式:
( a , θ 1 ) ( b , θ 2 ) = ( a b , θ 1 + θ 2 ) (a, \theta_1)(b,\theta_2)=(ab,\theta_1+\theta_2) (a,θ1)(b,θ2)=(ab,θ1+θ2)
以及复数运算满足交换律和结合律.
我们就可以以这些为铺垫,进行接下来的研究了 ?
定理
下面的定理不再证明,可以结合奆老博客和Mr.Du证明:
//提示:欧拉定理和复平面(上图)是个好东西.
ω n 1 = cos ( 2 π / n ) + i sin ( 2 π / n ) \omega_n^1=\cos(2\pi/n)+i\sin(2\pi/n) ωn1=cos(2π/n)+isin(2π/n)
对于 ω n \omega_n ωn序列, ω k = ω 1 k \omega_k=\omega_1^k ωk=ω1k(很实用)
ω 2 n 2 k = ω n k \omega_{2n}^{2k}=\omega_n^k ω2n2k=ωnk
ω n k + n 2 = − ω n k \omega_n^{k+\frac{n}{2}}=-\omega_n^k ωnk+2n=−ωnk
ω n k + n = ω n k \omega_n^{k+n}=\omega_n^k ωnk+n=ωnk
算法构建
经过这么一波操作之后,好像什么都没做…
因为直到现在,我们还是没有在这个数据上构造优化算法.
构造那么神奇的数据,不实现又有何用
要不我们先将一个
ω
n
\omega_n
ωn代入试试.
比如我们设定
A
(
n
)
=
∑
i
=
0
n
a
n
ω
n
i
A(n)=\sum_{i=0}^na_n\omega_n^i
A(n)=i=0∑nanωni
看上去没什么变化…
但是我把它的奇数指数和偶数指数分开,分别写成:
A
1
(
n
)
=
a
0
+
a
2
ω
n
+
.
.
.
+
a
n
ω
n
n
2
A_1(n)=a_0+a_2\omega_n+...+a_n\omega_n^\frac{n}{2}
A1(n)=a0+a2ωn+...+anωn2n
和
A
2
(
n
)
=
a
1
ω
n
+
a
3
ω
n
2
+
.
.
.
+
a
n
−
1
ω
n
n
2
A_2(n)=a_1\omega_n+a_3\omega_n^2+...+a_{n-1}\omega_n^\frac{n}{2}
A2(n)=a1ωn+a3ωn2+...+an−1ωn2n
//
n
n
n是
2
2
2的幂.
所以我们震惊的发现,竟然有这样的表达式!
A
(
x
)
=
A
1
(
x
2
)
+
ω
n
A
2
(
x
2
)
A(x)=A_1(x^2)+\omega_nA_2(x^2)
A(x)=A1(x2)+ωnA2(x2)
哇我真棒
但是,这还是没有用…
在仔细研究为什么需要使用复数.
因为,
ω
n
k
+
n
2
=
−
ω
n
k
\omega_n^{k+\frac{n}{2}}=-\omega_n^k
ωnk+2n=−ωnk!
所以,我们取
x
=
ω
n
k
+
n
2
x=\omega_n^{k+\frac{n}{2}}
x=ωnk+2n得:
A
(
x
)
=
A
1
(
x
2
)
−
ω
n
A
2
(
x
2
)
A(x)=A_1(x^2)-\omega_nA_2(x^2)
A(x)=A1(x2)−ωnA2(x2)
整理一下得(
k
<
n
2
k<\frac{n}{2}
k<2n):
A
(
ω
n
k
)
=
A
1
(
ω
n
2
k
)
+
ω
n
k
A
2
(
ω
n
2
k
)
A(\omega_n^k)=A_1(\omega_\frac{n}{2}^k)+\omega_n^kA_2(\omega_\frac{n}{2}^k)
A(ωnk)=A1(ω2nk)+ωnkA2(ω2nk)
A
(
ω
n
k
+
n
2
)
=
A
1
(
ω
n
2
k
)
−
ω
n
k
A
2
(
ω
n
2
k
)
A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_\frac{n}{2}^k)-\omega_n^kA_2(\omega_\frac{n}{2}^k)
A(ωnk+2n)=A1(ω2nk)−ωnkA2(ω2nk)
出现了!FFT!
IFFT
看似好像不那么容易.
但是我们已经会FFT了!
在分支向下的时候,保存每一项的系数,
具体怎么保存呢?
在转换的时候,我们乘的
ω
n
\omega_n
ωn.
转换回来的时候,就乘上它的共轭复数 What?你不会?
于是,我们就很快活
递归版FFT&IFFT
题目传送门
//为什么叫BF很快揭晓
先随便打个真·暴力上去…
#include<bits/stdc++.h>
using namespace std;
const int N=1e6+5;
long long a[N],b[N],c[N];
int n,m;
int main()
{
//True BF
//O(mn)...
scanf("%d%d",&m,&n);
for (int i=0;i<=m;i++) scanf("%d",&a[i]);
for (int i=0;i<=n;i++) scanf("%d",&b[i]);
for (int i=0;i<=m;i++)
for (int j=0;j<=n;j++)
c[i+j]+=a[i]*b[j];
for (int i=0;i<=n+m;i++)
printf("%d ",c[i]);
return 0;
}
Result:
//哇竟然能有55分!
而上面的诸多公式告诉我们,FFT可以完美解决这好像就是一道FFT模板题…
根据我们的FFT算法,很容易想到递归实现FFT:
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const double pi=acos(-1.0);
const int N=5e6+5;
struct complex
{
double x,y;
complex(double _x=0,double _y=0) {x=_x,y=_y;}
}a[N],b[N];
complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void FFT(int n,complex *a,int inv)
{
if (n==1) return;
int k=n>>1;
complex a1[k],a2[k];
for (int i=0;i<n;i+=2) a1[i>>1]=a[i],a2[i>>1]=a[i+1];//A1和A2
FFT(k,a1,inv);//
FFT(k,a2,inv);//分治
complex Wn=complex(cos(2.*pi/n),inv*sin(2.*pi/n)),w=complex(1,0);//omegaN和单位元
for (int i=0;i<k;i++,w=w*Wn)
{
a[i]=a1[i]+w*a2[i];
a[i+k]=a1[i]-w*a2[i];//一次推两个
}
}
int main()
{
int n,m;
scanf("%d%d",&n,&m);
for (int i=0;i<=n;i++) scanf("%lf",&a[i].x);
for (int i=0;i<=m;i++) scanf("%lf",&b[i].x);
int p=1;
while (p<=n+m) p<<=1;
FFT(p,a,1),FFT(p,b,1);//系数转点值
for (int i=0;i<=p;i++) a[i]=a[i]*b[i];
FFT(p,a,-1);//点值转系数
for (int i=0;i<=n+m;i++) printf("%d ",(int)(abs)(a[i].x/p+0.5));
return 0;
}
然而我竟然AC了!这是怎么回事?
Result:
这还怎么提出迭代版FFT!别提了
迭代版FFT&IFFT
接下来,我们将它优化一下:
这里是经过FFT分治后的序列
初始序列(id) | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
---|---|---|---|---|---|---|---|---|
1次操作 | 0 | 2 | 4 | 6 | 1 | 3 | 5 | 7 |
2次操作 | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 |
2进制编码 | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |
2进制逆序编码 | 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 |
我真聪明
可以发现,最后的编码是按逆序升序排列的.
换句话说:每个下标位置处理后的位置是它二进制逆序后的位置!
用这个短小精悍的代码就可以解决一切问题:
//名叫Rader算法
int p=1,L=0;
while (p<=n+m) p<<=1,L++;
for (int i=0;i<p;i++)
order[i]=(order[i>>1]>>1)|((i&1)<<(L-1));
//(i&1)<<(L-1) 是为了考虑后面的1<<(L-1)个数
真有趣 ?
接下来,我们就更快活开心了 ?
蝴蝶效应
属于小优化的范畴.
主要是因为复数乘法过于缓慢,所以可以预先处理好.
而已qwq.不要想到龙卷风那里去了
Code
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=5e6+5;
const double pi=acos(-1.);
struct complex
{
double x,y;
complex(double _x=0,double _y=0) {x=_x,y=_y;}
}a[N],b[N];
complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
int order[N];
int p=1,L;
void FFT(complex *a,int inv)
{
for (int i=0;i<p;i++)
if (i<order[i]) swap(a[i],a[order[i]]);
for (int l=1;l<p;l<<=1)
{
complex Wn(cos(pi/l),inv*sin(pi/l));
for (int R=l<<1,j=0;j<p;j+=R)
{
complex w(1,0);
for (int k=0;k<l;k++,w=w*Wn)
{
complex x=a[j+k],y=w*a[j+l+k];//B
a[j+k]=x+y;
a[j+l+k]=x-y;
}
}
}
}
int main()
{
int n,m;
scanf("%d%d",&n,&m);
for (int i=0;i<=n;i++) scanf("%lf",&a[i].x);
for (int i=0;i<=m;i++) scanf("%lf",&b[i].x);
while (p<=n+m) p<<=1,L++;
for (int i=0;i<p;i++)
order[i]=(order[i>>1]>>1)|((i&1)<<(L-1));
FFT(a,1),FFT(b,1);
for (int i=0;i<=p;i++) a[i]=a[i]*b[i];
FFT(a,-1);
for (int i=0;i<=n+m;i++)
printf("%d ",(int)(a[i].x/p+0.5));
return 0;
}
Result:
感谢奆老关注 qwq ?
后记
FFT,又称Fast-Fast-TLE…
这里为什么不直接从
<
c
o
m
p
l
e
x
>
<complex>
<complex>这个库里调用
c
o
m
p
l
e
x
complex
complex?
因为…太慢了…