uoj模板题库T1
传说中的快速傅里叶变换
以前在算法导论上看了几眼然后放弃了。。
不过静下心来读。。加上问了李学霸关于\(n\)次单位复数根的事情,大概\(40\)分钟搞定
然后代码实现就是蒯算法导论和hzwer犇的事了。。理解了还是很好码的
————————————-开始讲课————————————-
一、多项式乘法
给出两个多项式\(A=\sum_{i=0}^n a_i x^i,B=\sum_{i=0}^m b_i x^i\)
求\(A\times B\)
朴素的方法就是各位相乘,复杂度\(O(nm)\)。有没有更快的方法呢?
首先介绍一下多项式的点值表示法:
对于一个多项式\(A=\sum_{i=0}^n a_i x^i\)
取\(n\)个不同的\(x\)值,算出\(n\)个不同的\(y\)值,得到\(n\)个点\((x_1,y_1),(x_2,y_2)…\)这些点唯一确定了多项式\(A\)。用霍纳算法可以在\(O(n^2)\)的复杂度内得到一个多项式的点值表示。由系数表示法转为点值表示法的过程称为求值;其逆过程称为插值。
而点值表示法的好处为:取同样的x求值,多项式乘法的复杂度为\(O(n)\):只要将两个多项式的\(y\)相乘即可。
则多项式乘法的过程可以转换为:对\(A\),\(B\)求值->点值乘法->插值
若能将求值和插值的复杂度降低,就能达到我们的目的了!
二、FFT
【【声明:以下规定\(n\)为\(2\)的幂】】【你问我\(n\)不是\(2\)的幂怎么办?扩大一下,高位系数全为\(0\)不就完了
FFT的核心思想:通过恰当选取\(x\)的值,并采用分治策略使得求值和插值的复杂度降为\(O(nlogn)\)。
首先我们要了解的是\(n\)次单位复数根(%数学组李浚鼎大薛霸)
\(n\)次单位复数根:表示为\(ω_n,ω_n^n=1\)
并且\(n\)次单位复数根的个数为\(n\)
算法导论告诉我们,\(n\)个单位复数根均匀的分布在以复平面的原点为圆心的单位半径的圆周上

再举个例子。。\(4\)次单位复数根就是\(±1\)和复数中喜闻乐见的\(±i\)
那么究竟怎么求呢
有公式\(ω_n^k=e^{2πi\frac{k}{n}}\)
然后又根据\(e^{iu}=cos(u)+isin(u)\)
于是\(ω_n^k=cos(2πk/n)+isin(2πk/n)\)
代码就好写了是吧
那么这玩意有什么用呢?
再给出两个引理:
消去引理:\(ω_{ak}^{bk}=e^{2πi\frac{bk}{ak}}=e^{2πi\frac{b}{a}}=ω_a^b\)
折半引理:由消去引理易证,如果\(n\)为大于\(0\)的偶数,那么\(ω_n\)平方的集合为\(ω_{n/2}\)的集合
即,\(ω_n\)平方的集合中元素的个数仅为\(\frac{n}{2}\)
故可以进行分治
将\(A=\sum_{i=0}^n a_i x^i\)拆分为
\(A_0={a_0,a_2x,a_4x^2,…,a_{n-2}^{n/2-1}}\)(所有下标二进制末尾为\(0\)的系数组成的多项式)
\(A_1={a_1,a_3x,a_5x^2,…,a_{n-1}^{n/2-1}}\)(所有下标二进制末尾为\(1\)的系数组成的多项式)
则\(A(x)=A_0(x^2)+xA_1(x^2)\)
那么若我们代入A的各个x值为n个n次单位复数根。。则\(x^2\)的取值仅有\(\frac{n}{2}\)个。。
那么两个子问题\(A_0(x^2),A_1(x^2)\)的规模都只有原问题的一半。。
Congratulations! 复杂度已经被我们降到了\(O((n+m)log(n+m))\)。
至于插值,只要将\(ω_n变成ω_n^{-1}\),再将所得结果都除以\(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 43 44 45 46 47 48 49 50 |
// <FFT.cpp> - 02/23/16 20:45:22 // This file is created by XuYike's black technology automatically. // Copyright (C) 2015 ChangJun High School, Inc. // I don't know what this program is. #include <iostream> #include <vector> #include <algorithm> #include <cstring> #include <cstdio> #include <cmath> #include <complex> using namespace std; typedef long long lol; typedef complex<double> C; int gi(){ int res=0,fh=1;char ch=getchar(); while((ch>'9'||ch<'0')&&ch!='-')ch=getchar(); if(ch=='-')fh=-1,ch=getchar(); while(ch>='0'&&ch<='9')res=res*10+ch-'0',ch=getchar(); return fh*res; } const int N=262145; const double pi=acos(-1); int n,m; C a[N],b[N]; void fft(C *a,int n,int f){ if(n==1)return; C wn(cos(2.0*pi/n),sin(f*2.0*pi/n)),w(1,0),t; C a0[n>>1],a1[n>>1]; for(int i=0;i<n>>1;i++)a0[i]=a[i<<1],a1[i]=a[i<<1|1]; fft(a0,n>>1,f);fft(a1,n>>1,f); for(int i=0;i<n>>1;i++,w*=wn){ t=w*a1[i]; a[i]=a0[i]+t; a[i+(n>>1)]=a0[i]-t;//由于ω_n图像的对称性,想想就明白了 } } int main(){ n=gi();m=gi(); for(int i=0;i<=n;i++)a[i]=gi(); for(int i=0;i<=m;i++)b[i]=gi(); m+=n; for(n=1;n<=m;n<<=1); fft(a,n,1);fft(b,n,1); for(int i=0;i<=n;i++)a[i]*=b[i]; fft(a,n,-1); for(int i=0;i<=m;i++)printf("%d ",int(a[i].real()/n+0.5));//+0.5是为了防止精度丢失造成结果略小,小数部分被舍掉 return 0; } |
其中核心的FFT函数也只有\(9\)行呢。。有点不可思议的感觉
这样我们就能够AC了,总时间将近\(4\)秒的样子
但是!我们要精益求精是不是!
既然是这么simple的递归就肯定能转成非递归版本的,常数更小对不对!
那我们来分析一下这个递归的结构:
我们只要把原数组先排成最底下一排的顺序再两两合并就可以自底向上地实现了
那么观察一下最后一排的顺序:发现把每一个二进制数前后颠倒即为原数组
那么就很简单了(详见代码)
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 |
// <FFT.cpp> - 02/23/16 20:45:22 // This file is created by XuYike's black technology automatically. // Copyright (C) 2015 ChangJun High School, Inc. // I don't know what this program is. #include <iostream> #include <vector> #include <algorithm> #include <cstring> #include <cstdio> #include <cmath> #include <complex> using namespace std; typedef long long lol; typedef complex<double> C; int gi(){ int res=0,fh=1;char ch=getchar(); while((ch>'9'||ch<'0')&&ch!='-')ch=getchar(); if(ch=='-')fh=-1,ch=getchar(); while(ch>='0'&&ch<='9')res=res*10+ch-'0',ch=getchar(); return fh*res; } const int MAXN=262145; const double pi=acos(-1); int n,m,L; int R[MAXN]; C a[MAXN],b[MAXN]; void fft(C *a,int f){ for(int i=0;i<n;i++)if(i<R[i])swap(a[i],a[R[i]]);//交换位置 for(int i=1;i<n;i<<=1){//待合并区间长度 C wn(cos(pi/i),sin(f*pi/i)),x,y;//这里就不用再*2了,因为合并后的区间长度是i的两倍 for(int j=0;j<n;j+=i<<1){//起始位置 C w(1,0); for(int k=0;k<i;k++,w*=wn){//第k个 x=a[j+k];y=w*a[j+i+k]; a[j+k]=x+y; a[j+i+k]=x-y; } } } } int main(){ n=gi();m=gi(); for(int i=0;i<=n;i++)a[i]=gi(); for(int i=0;i<=m;i++)b[i]=gi(); m+=n; for(n=1;n<=m;n<<=1)L++; for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));//递推 fft(a,1);fft(b,1); for(int i=0;i<=n;i++)a[i]*=b[i]; fft(a,-1); for(int i=0;i<=m;i++)printf("%d ",int(a[i].real()/n+0.5)); return 0; } |
泥又在刷火题
不错,不错
劲啊劲啊~~~参看大神的模版
光光好劲啊