模板题
点值表达
大莉模拟的复杂度是
O
(
n
2
)
O(n^2)
O(n2)。
引入点值表达qwq(这部分摘抄自毒瘤czh的FFT讲稿):
1.
1.
1.例子:
A
(
x
)
=
x
2
+
2
x
−
1
A(x)=x^2+2x-1
A(x)=x2+2x−1可以被表达为
{
(
0
,
−
1
)
,
(
1
,
2
)
,
(
2
,
7
)
}
\{(0 , -1), (1 , 2),(2 , 7)\}
{(0,−1),(1,2),(2,7)}
2. 2. 2.用处:在多项式乘法中,系数表达暴力模拟复杂度 O ( n 2 ) O(n^2) O(n2),但是点值表达式可以 O ( n ) O(n) O(n)求出
更加具体地说,我们有两个点值表达
A
(
x
)
=
{
(
x
0
,
A
(
x
0
)
)
,
(
x
1
,
A
(
x
1
)
)
,
(
x
2
,
A
(
x
2
)
)
,
.
.
.
.
,
(
x
n
,
A
(
x
n
)
)
}
A(x)=\{(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2)),....,(x_n,A(x_n))\}
A(x)={(x0,A(x0)),(x1,A(x1)),(x2,A(x2)),....,(xn,A(xn))}
B
(
x
)
=
{
(
x
0
,
B
(
x
0
)
)
,
(
x
1
,
B
(
x
1
)
)
,
(
x
2
,
B
(
x
2
)
)
,
.
.
.
.
,
(
x
n
,
B
(
x
n
)
)
}
B(x)=\{(x_0,B(x_0)),(x_1,B(x_1)),(x_2,B(x_2)),....,(x_n,B(x_n))\}
B(x)={(x0,B(x0)),(x1,B(x1)),(x2,B(x2)),....,(xn,B(xn))}
乘积
C
(
x
)
=
A
(
x
)
×
B
(
x
)
=
{
(
x
0
,
A
(
x
0
)
×
B
(
x
0
)
)
,
(
x
1
,
A
(
x
1
)
×
B
(
x
1
)
)
,
(
x
2
,
A
(
x
2
)
×
B
(
x
2
)
)
,
.
.
.
.
,
(
x
n
,
A
(
x
n
)
×
B
(
x
n
)
)
}
C(x)=A(x)\times B(x)=\{(x_0,A(x_0)\times B(x_0)),(x_1,A(x_1)\times B(x_1)),(x_2,A(x_2)\times B(x_2)),....,(x_n,A(x_n)\times B(x_n))\}
C(x)=A(x)×B(x)={(x0,A(x0)×B(x0)),(x1,A(x1)×B(x1)),(x2,A(x2)×B(x2)),....,(xn,A(xn)×B(xn))}
所以我们可以通过
O
(
n
)
O(n)
O(n)复杂度,求出两个点值表达式乘积的结果(两个多项式的乘积又称卷积)
知道了点值表达式就珂以确定一个多项式了qwq(珂以列方程用高斯消元爆解一波qwq)
其中DFT算法是暴力 O ( n 2 ) O(n^2) O(n2)把系数表达转点值表达,FFT算法是 O ( n l o g n ) O(nlogn) O(nlogn)把系数表达转点值表达
DFT(离散傅里叶变换)
普通做法:随机取
n
n
n个值代入求值qwq。
DFT是考虑建立一个复平面,在上面怼出 画一个单位圆:
因为是单位圆,所以圆上的每一个点对应的复数模均为1。
把这个圆平均分成
n
n
n份,把每个点对应的复数代入。
具体地说,这个函数的点值表达是
{
(
w
n
0
,
f
(
w
n
0
)
)
,
(
w
n
1
,
f
(
w
n
1
)
)
,
.
.
.
,
(
w
n
n
−
1
,
f
(
w
n
n
−
1
)
)
}
\{(w_n^0,f(w_n^0)),(w_n^1,f(w_n^1)),...,(w_n^{n-1},f(w_n^{n-1}))\}
{(wn0,f(wn0)),(wn1,f(wn1)),...,(wnn−1,f(wnn−1))},
其中
w
n
k
=
c
o
s
(
k
n
2
π
)
+
i
s
i
n
(
k
n
2
π
)
w_n^k=cos(\frac{k}{n}2\pi)+isin(\frac{k}{n}2\pi)
wnk=cos(nk2π)+isin(nk2π)
FFT(快速傅里叶变换)
w w w的性质
证明珂以看巨神mhy的博客qwq
这里只列出性质,就不证明了:
1.
w
n
k
=
w
2
n
2
k
1.w_n^k=w_{2n}^{2k}
1.wnk=w2n2k
2.
w
n
k
+
n
2
=
−
w
n
k
2.w_n^{k+\frac{n}{2}}=-w_n^k
2.wnk+2n=−wnk
3.
w
n
0
=
w
n
n
=
1
3.w_n^0=w_n^n=1
3.wn0=wnn=1
4.
w
n
p
∗
w
n
q
=
w
n
p
+
q
4.w_n^p*w_n^q=w_n^{p+q}
4.wnp∗wnq=wnp+q
表达式的变形
首先把一个函数用系数表达式写出来(假设
n
n
n为偶数,如果不是就补一个
0
∗
x
n
0*x^n
0∗xn):
f
(
x
)
=
a
0
+
a
1
x
1
+
a
2
x
2
+
a
3
x
3
+
a
4
x
4
+
.
.
.
+
a
n
−
1
x
n
−
1
f(x)=a_0+a_1x^1+a_2x^2+a_3x^3+a_4x^4+...+a_{n-1}x^{n-1}
f(x)=a0+a1x1+a2x2+a3x3+a4x4+...+an−1xn−1
按照下标的奇偶性把
f
f
f分成两半:
f
(
x
)
=
(
a
0
+
a
2
x
2
+
.
.
.
+
a
n
−
2
x
n
−
2
)
+
x
(
a
1
+
a
3
x
2
+
.
.
.
+
a
n
−
1
x
n
−
2
)
f(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})
f(x)=(a0+a2x2+...+an−2xn−2)+x(a1+a3x2+...+an−1xn−2)
令
A
1
(
x
)
=
a
0
+
a
2
x
+
a
4
x
2
+
.
.
.
+
a
n
−
2
x
n
−
2
2
A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n-2}{2}}
A1(x)=a0+a2x+a4x2+...+an−2x2n−2,
A
2
(
x
)
=
a
1
+
a
3
x
+
a
5
x
2
+
.
.
.
+
a
n
−
1
x
n
−
2
2
A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n-2}{2}}
A2(x)=a1+a3x+a5x2+...+an−1x2n−2
不难发现
f
(
x
)
=
A
1
(
x
2
)
+
x
A
2
(
x
2
)
f(x)=A_1(x^2)+xA_2(x^2)
f(x)=A1(x2)+xA2(x2)。
以下分别把
w
n
k
w_n^k
wnk和
w
n
k
+
n
2
w_n^{k+\frac{n}{2}}
wnk+2n代入:
把
w
n
k
w_n^k
wnk代入,得到
f
(
w
n
k
)
=
A
1
(
w
n
k
∗
w
n
k
)
+
x
A
2
(
w
n
k
∗
w
n
k
)
=
A
1
(
w
n
2
k
)
+
w
n
k
A
2
(
w
n
2
k
)
f(w_n^k)=A_1(w_n^k*w_n^k)+xA_2(w_n^k*w_n^k)=A_1(w_n^{2k})+w_n^kA_2(w_n^{2k})
f(wnk)=A1(wnk∗wnk)+xA2(wnk∗wnk)=A1(wn2k)+wnkA2(wn2k)
由
w
w
w的性质1,珂以得到:
f
(
w
n
k
)
=
A
1
(
w
n
2
k
)
+
w
n
k
A
2
(
w
n
2
k
)
f(w_n^k)=A_1(w_{\frac{n}{2}}^k)+w_n^kA_2(w_{\frac{n}{2}}^k)
f(wnk)=A1(w2nk)+wnkA2(w2nk)
把
w
n
k
+
n
2
w_n^{k+\frac{n}{2}}
wnk+2n代入,得到:
f
(
w
n
k
+
n
2
)
=
A
1
(
w
n
2
k
+
n
)
+
w
n
k
+
n
2
A
2
(
w
n
2
k
+
n
)
=
A
1
(
w
n
2
k
∗
w
n
n
)
−
w
n
k
A
2
(
w
n
2
k
∗
w
n
n
)
=
A
1
(
w
n
2
k
)
−
w
n
k
A
2
(
w
n
2
k
)
f(w_n^{k+\frac{n}{2}})=A_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}A_2(w_n^{2k+n})=A_1(w_n^{2k}*w_n^n)-w_n^kA_2(w_n^{2k}*w_n^n)=A_1(w_n^{2k})-w_n^kA_2(w_n^{2k})
f(wnk+2n)=A1(wn2k+n)+wnk+2nA2(wn2k+n)=A1(wn2k∗wnn)−wnkA2(wn2k∗wnn)=A1(wn2k)−wnkA2(wn2k)
分治写法
认真观察……这
f
(
w
n
k
)
f(w_n^k)
f(wnk)和
f
(
w
n
k
+
n
2
)
f(w_n^{k+\frac{n}{2}})
f(wnk+2n)长得很像?
于是珂以递归(分治)解决,每次求解左半部分的答案,再推出右半部分的答案:
void FFT(complex<double> *ans,int n,int type) {
//type=1表示FFT;type=-1表示IFFT
if(n==1) return;
//偶数存L,奇数存R
complex<double> L[(n>>1)+1],R[(n>>1)+1];
for(re i=0; i<n; i+=2) {
L[i>>1]=ans[i];
R[i>>1]=ans[i+1];
}
FFT(L,n>>1,type);
FFT(R,n>>1,type);
//玄学单位圆乱搞
complex<double> w(1,0);
complex<double> k(cos(2.0*pi/n),type*sin(2.0*pi/n));
int maxn=n>>1;
for(re i=0; i<maxn; i++) {
complex<double> tmp=w*R[i];
ans[i]=L[i]+tmp;
ans[i+maxn]=L[i]-tmp;
w*=k;
}
}
然后交到例题上,TLE??
递归常数蜃大qwq,考虑使用迭代。
迭代写法
观察下面这张图:
发现原序列的每个数按对称轴翻转一下就是后序列的数qwq,称之为位逆序置换。
所以,我们可以预处理出系数位置递归到最底层的位置,自底向上进行递推。
处理出位逆序置换之后的方法(
R
[
i
]
R[i]
R[i]表示
i
i
i进行位逆序置换之后得到的值,
L
L
L表示二进制下的总位数):
for(re i=0; i<MAXN; i++) { //MAXN表示要处理出来的数的总数
//表示先把最后一位去掉,剩下的部分进行反转,然后把最后一位反转到最高位qwq
R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
}
所以把递归强行用循环写出来:
第一维枚举区间长度,即递归版本中的
n
n
n。第二维枚举区间开始位置。第三维枚举当前数是区间中的第几个。
然后珂以大莉写一波qwq(我的版本中,区间长度实际上为
2
i
2i
2i):
void FFT(Complex *ans,int n,int type) {
//type=1表示FFT,type=-1表示IFFT
//把序列交换成后序列的形式
for(re i=0; i<n; i++) if(i<R[i]) swap(ans[i],ans[R[i]]);
for(re i=1; i<n; i<<=1) {
//(2*pi)/(2*i)=pi/i
Complex wn=(Complex){cos(pi/(double)i),type*sin(pi/(double)i)};
for(re j=0; j<n; j+=(i<<1)) {
Complex w=(Complex){1,0};
for(re k=0; k<i; k++) {
//类似递归版的处理
const Complex x=ans[j+k],y=w*ans[i+j+k];
ans[j+k]=x+y;
ans[i+j+k]=x-y;
w=w*wn;
}
}
}
}
IFFT(快速傅里叶逆变换)
毒瘤czh:
美妙结论:一个多项式在分治的过程中乘上单位根的共轭复数,分治完的每一项除以n即为原多项式的每一项系数(具体看程序实现,证明略)
也就是在调用的时候把type设成-1就珂以了
例题完整代码
至此,就珂以写出例题了(注意预处理和边界条件qwq):
#include<stdio.h>
#include<cstring>
#include<algorithm>
#include<math.h>
#define re register int
#define rl register ll
#define pi 3.14159265358979323846264338328
using namespace std;
int read() {
re x=0,f=1;
char ch=getchar();
while(ch<'0' || ch>'9') {
if(ch=='-') f=-1;
ch=getchar();
}
while(ch>='0' && ch<='9') {
x=10*x+ch-'0';
ch=getchar();
}
return x*f;
}
inline void write(const int x) {
if(x>9) write(x/10);
putchar(x%10+'0');
}
const int Size=4000005;
//手写complex类,不过貌似也没快多少
struct Complex {
double x,y;
inline double real() {
return x;
}
};
inline Complex operator + (const Complex a,const Complex b) {
return (Complex){a.x+b.x,a.y+b.y};
}
inline Complex operator - (const Complex a,const Complex b) {
return (Complex){a.x-b.x,a.y-b.y};
}
inline Complex operator * (const Complex a,const Complex b) {
return (Complex){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}
Complex A[Size],B[Size];
int R[Size];
void FFT(Complex *ans,const int n,const int type) {
//type=1表示FFT,type=-1表示IFFT
for(re i=0; i<n; i++) if(i<R[i]) swap(ans[i],ans[R[i]]);
for(re i=1; i<n; i<<=1) {
//(2*pi)/(2*i)=pi/i
const Complex wn=(Complex){cos(pi/i),type*sin(pi/i)};
for(re j=0; j<n; j+=i<<1) {
Complex w=(Complex){1,0};
for(re k=0; k<i; k++) {
Complex x=ans[j+k],y=w*ans[i+j+k];
ans[j+k]=x+y;
ans[i+j+k]=x-y;
w=w*wn;
}
}
}
}
int ans[Size];
int main() {
int n=read();
int m=read();
for(re i=0; i<=n; i++) {
A[i]=(Complex){read(),0};
}
for(re i=0; i<=m; i++) {
B[i]=(Complex){read(),0};
}
//now为大于n+m的第一个2的幂,因为总区间长度必须是2的幂qwq
int now=1,L=0;
while(now<=n+m) {
now<<=1;
L++;
}
for(re i=0; i<now; i++) {
R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
}
FFT(A,now,1);
FFT(B,now,1);
for(re i=0; i<=now; i++) {
A[i]=A[i]*B[i];
}
FFT(A,now,-1);
for(re i=0; i<=n+m; i++) {
write(int(A[i].real()/now+0.5)); //注意四舍五入
putchar(' ');
}
return 0;
}
其他题目
洛谷P1919 A*B Problem升级版
即bzoj2179 (权限题)快速傅里叶变换
bzoj4827 礼物
……
…………
………………
题解?咕咕咕
还是有一篇的qwq