如果说出题人给了你一个卷积的话,那么你应该也友好的回赠他一个FFT;
-------------FFT
·一个问题:FFT的作用是:给定两个多项式,在nlgn内求出其乘积多项式。
·多项式的两种表达:
多项式$\sum_{i = 0}^{n-1}a_i*x^i$可以用系数$(a_1,a_2,a_3,…,a_{n-1})$来表示;也可以用点值${(x_0,y_0),(x_1,y_1),…,(x_{n-1},y_{n-1})}$,其中$y_i$为$x_i$处对应的值。可以,一个次数界为n的多项式需要至少n个点值确定。(少了正确性难以保证) 对于多项式的运算,可以先求出各个多项式在某些点的点值,转化成对应点值的运算,求出最后的点值向量,最后插值求出结果。上述是FFT的基本思路,或者说FFT进行了上述两步转化。
·关于复数:
①复数的指数形式定义: $e^{\pm iu}=cos(u) \pm isin(u)$
:
$ e^x = 1 +\frac{x}{1!} + \frac{x^2}{2!} + \frac{x^3}{3!} + … $
$cos(x) = 1 – \frac{x^2}{2!} + \frac{x^4}{4!} - \frac{x^6}{6!} + … $
$sin(x) = x – \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + … $
令x = ±iu ,代换比较一下就OJBK了。
取u=π,,我们于是激动地有了一个很美的公式:$$ e^{\pi i} + 1 = 0 $$
②:单位复数根是满足$\omega ^n = 1$的复数$\omega $。利用上面的结论,n次单位根恰好有n个,这些根为 $$ e^{2k \pi i /n}, k \in {0,1,...,n-1} $$ 。其中k=1时,$ \omega_{n} = e^{2 \pi i /n} $叫做主n次单位根。
我们在上的点是$\omega_{n}$ ( cos( $\frac{2k\pi}{n}$ ) , sin( $\frac{2k\pi}{n}$ ) ),这样我们可以得到如上图的群,这个群再复平面的单位圆上均匀分布,有这样的结论$$ \omega_{n}^{j}\omega_{n}^{k} = \omega_{n}^{j+k} = \omega_{n}^{(j+k) mod n} $$。
n次单位根有一些性质:①消去引理:对任何整数n≥0,k≥0,以及d>0,有$$\omega_{dn}^{dk} = \omega_{n}^{k}$$(证明:大米兔式显然)
②折半引理:如果n为大于零的偶数,则{$(\omega_n^0)^2 , (\omega_n^1)^2 , (\omega_n^3)^2... , (\omega_n^{n-1})^2 $} = {$ \omega_{n/2}^{0} , \omega_{n/2}^{1} , \omega_{n/2}^{2} , ... , \omega_{n/2}^{n/2 - 1} $} , 因为 $(\omega_n^{k + n/2})^2 = \omega_n^{2k + n} = (\omega_n^{k})^2 = \omega_n^{2k} = \omega_{n/2}^{k}$ ,也可以用上面的复平面理解,多滚一遍就好了。
③求和引理:对任意整数n≥1和不能被n整除的非负整数k,有$\sum_{j=0}^{n-1}(\omega_n^k)^j = 0$;因为根据等比数列求和,大米兔发现:$\sum_{j=0}^{n-1}(\omega_{n}^{k})^j = \frac{(\omega_n^k)^n - 1}{\omega_{n}^{k} - 1} = \frac{(\omega_n^n)^k - 1}{\omega_n^k - 1} = 0$(所以你现在应该反映过来为什么$k \nmid n$了)
(有公式恐惧症的朋友坚持啊。。。。)
DFT
对于系数为a=$(a_0,a_1,a_2,…,a_{n-1})$的次数界为n的多项式:$$A(x) = \sum_{j=0}^{n-1}a_{j} x^{j}$$,对求值$(\omega_n^0,\omega_n^1,…,\omega_n^{n-1})$的结果定义为$$y_k = A(\omega_n^k) = \sum_{j=0}^{n-1}\omega_{n}^{jk} ①$$,记向量$y=(y_0,y_1,y_2,...,y_{n-1})$为系数向量a的离散傅里叶变换(DFT),即$y=DFT_{n}(a)$;
插值
如果我们已经得到了,DFT,可以写成如下形式:y = V*a,V即:$$\begin{bmatrix} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2n-2} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3n-3} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{bmatrix}$$ ,
什么!你竟然已经发现了$V_{j1j2} = \omega_{n}^{j1j2}$;那我们就可以尝
试解这个东东,V的逆为: $V^{-1}_{j1j2} = \omega_{n}^{-j1j2} /n$ ,因为$$ VV^{-1} = \sum_{k=0}^{n-1} (\omega_n^{j1k} * \omega_n^{-kj2} /n) = \sum_{k=0}^{n-1} (\omega_n^{(j1-j2)k}/n)$$当j1=j2时值为1,否则n不整除(j1-j2),根据求和引理,值为0,所以$VV^{-1}$为单位矩阵。
$I * y = V * a –> V * y^{-1}$ = a (I为单位矩阵
)所以$$ a_j = \frac{1}{n} \sum_{k=0}^{n-1} y_k \omega_n^{-jk} ② $$
对比①式和②式(你应该还找得到①式吧。。。),如果传说中的FFT可以算出①,那么把 $ \omega_n^{kj} $ 换成 $ \omega_n^{-kj} $ 最后的结果再除n,就可以算出②,或者说DFT的逆。分析到这里,快速傅里叶变换(FFT)已经呼之欲出了!
FFT
让我们再来悠闲地回顾下FFT要做的事情,对于多项式A(x),求:$$y_k = A(\omega_n^k) = \sum_{j=0}^{n-1}\omega_{n}^{jk} ①$$ 。
基本思路是分治使其达到nlgn的复杂度,先将n填充到2的整数次幂,然后假设我们计算$$A(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + … a_{n-1} x^{n-1}$$ 令,$$A_0(x) = a_0 + a_2x + a_4x^2 + … + a_{n-2}x^{n/2-1}$$ $$A_1(x) = a_1 + a_3x + a_5x^2 + … + a_{n-1} x^{n/2-1}$$emm。。。$$A(x) = A_0(x^2) + x*A_1(x^2)$$ 为什么要机关算尽地去平方呢?是这样的,根据折半定理,x^2 才可以保证我们需要计算的规模由n变为n/2,递推地去求就好了。
但是有一个问题,原先的数组顺序需要被改变,但是这不是问题。如下图:
第一层会选取二进制末尾为0的数到左边,第二层会选取二进制倒数第二位为0的到左边,所以,将所有数二进制翻转,从小到大写即可!
具体的,用f表示是计算DFT还是$DFT^{-1}$我们枚举计算的规模的一半i,计算出此时的主n次单位根:
接着枚举每一个区间的开头j,初始化$\omega$,最后枚举具体的位置j+k; 我们知道:$$\omega_n^{k} = – \omega_n^{k + n/2}$$,所以A(j+k)=x+y,A(j+i+k) = x - y;
最后是因为如果f==-1不能忘了除n;
就这样。。。。
因为精疲力尽的大米兔下午就开学了,所以只能附上一些简单的练习和代码:
(请自动忽略每个代码前的注释)
luogu 3803 A * B1 /*数据较强的良心模板题; 2 3 */ 4 #include5 #include 6 #include 7 #include 8 #define N 3000010 9 using namespace std;10 const double pi = acos(-1);11 typedef complex C;12 int n,m,len,R[N];13 C a[N],b[N];14 char gc() {15 static char *p1,*p2,s[1000000];16 if(p1==p2) p2 = (p1=s) + fread(s,1,1000000,stdin);17 return (p1==p2)?EOF:*p1++;18 }19 int read() {20 int x = 0;21 char c = gc();22 while(c < '0' || c > '9') c = gc();23 while(c >= '0' && c <= '9') x = x * 10 + c - '0',c = gc();24 return x;25 }26 void FFT(C *a,int f) {27 for(int i = 0; i < len; i++) if(i < R[i]) swap(a[i],a[R[i]]);28 for(int i = 1; i < len; i<<=1) {29 C wn(cos(pi/i),f * sin(pi/i));30 for(int j = 0; j < len; j+=(i<<1)) {31 C w(1,0);32 for(int k = 0; k < i; k++,w*=wn) {33 C x = a[j + k],y = w * a[j + i + k];34 a[j + k] = x + y,a[j + i + k] = x - y;35 }36 }37 }38 if(f==-1) for(int i = 0; i < len; i++) a[i] /= len;39 }40 int main() {41 freopen("P3803.in","r",stdin);42 freopen("P3803.out","w",stdout);43 n = read();44 m = read();45 for(int i = 0; i <= n; i++) a[i] = read();46 for(int i = 0; i <= m; i++) b[i] = read();47 int t = n + m;48 int L = 0;49 for(len=1; len<=t; len<<=1) L++;50 for(int i = 1; i < len; i++) R[i] = (R[i>>1]>>1) | ((i&1)<<(L-1));51 for(int i = n + 1; i < len; i++) a[i] = 0;52 for(int i = m + 1; i < len; i++) b[i] = 0;53 FFT(a,1);54 FFT(b,1);55 for(int i = 0; i < len; i++) a[i] *= b[i];56 FFT(a,-1);57 for(int i = 0; i <= t; i++) {58 printf("%d ",(int)(a[i].real() + 0.1));59 }60 printf("\n");61 return 0;62 }//by tkys_Austin;
BZOJ 3527 [Zjoi2014]力
1 /*化式子很重要; 2 分开求sigma 3 把第二个化成sum_{j=0}{n-1-i}q_{n-1-i-j}p_{j} 并把原数组翻转; 4 其实就是倒着来就与第一个一样了。 5 C输出时要取real或img,不然自己的电脑跑出来没事,BZOJ上回RE,小心点好; 6 好!决定了,下一个卷积我一定要自己化出来! 7 */ 8 #include9 #include 10 #include 11 #include 12 #include 13 #define N 26214514 using namespace std;15 typedef complex C;16 const double Pi = acos(-1);17 int n,m,L,rev[N];18 C a[N],b[N],c[N];19 void FFT(C *a,int f){20 for(int i = 1;i < n;i++) if(i < rev[i]) swap(a[i],a[rev[i]]);21 for(int i = 1;i < n;i<<=1){22 C wn(cos(Pi/i),f * sin(Pi/i)); 23 for(int j = 0;j < n;j+=(i<<1)){24 C w(1,0);25 for(int k = 0;k < i;k++,w*=wn){26 C x = a[j + k],y = w * a[j + k + i];27 a[j + k] = x + y,a[j + k + i] = x - y;28 }29 }30 }31 if(f==-1) for(int i = 0;i < n;i++) a[i]/=n;32 }33 int main()34 { //freopen("BZOJ3527.in","r",stdin);35 //freopen("BZOJ3527.out","w",stdout);36 scanf("%d",&n);n--;37 for(int i = 0;i <= n;i++){38 double x;scanf("%lf",&x);39 a[i] = b[n - i] = x; 40 }41 for(int i = 1;i <= n;i++) c[i] = 1.0 / i / i;42 int len = 2 * n ; for(n=1;n<=len;n<<=1)L++;43 for(int i = 1;i < n;i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(L-1));44 FFT(a,1); FFT(b,1); FFT(c,1);45 for(int i = 0;i < n;i++) a[i] *= c[i],b[i] *= c[i];46 FFT(a,-1); FFT(b,-1);47 len/=2;48 for(int i = 0;i <= len;i++){49 printf("%.3lf\n",a[i].real() - b[len - i].real());50 }51 return 0;52 }//by tkys_Austin;
BZOJ 3513 [MUTC2013]idiots
1 /*首先,能被我秒掉的题还真的不多; 2 其次,complex我以后再也不用了; 3 因为数值只有1e5左右,所以可以考虑三角不等式其实可以和卷积联系起来; 4 但是记住正难则反,此题直接统计ans比较麻烦(其实也不是很麻烦),计算不合法的不怎么容易写错; 5 L不清零又遭了、、、 6 */ 7 #include8 #include 9 #include 10 #include 11 #include 12 #define N 26214513 #define ll long long14 using namespace std;15 const double Pi = acos(-1);16 struct C {17 double x,y;18 C(double x = 0,double y = 0):x(x),y(y) {}19 // bool operator =(const C &a)const {*this = a;}20 // bool operator =(const int &a)const {*this = C(a,0);}21 C operator +(const C &a) { return C(x + a.x,y + a.y);}22 C operator -(const C &a) { return C(x - a.x,y - a.y);}23 C operator *(const C &a) { return C(x * a.x - y * a.y,x * a.y + y * a.x);}24 C operator /(const int &a) { return C(x/a,y/a);}25 };26 int n,m,L;27 ll a[N],rev[N];28 C b[N];29 char gc(){30 static char *p1,*p2,s[1000000];31 if(p1==p2) p2 = (p1 = s) + fread(s,1,1000000,stdin);32 return (p1==p2)?EOF:*p1++;33 }34 int read(){35 int x = 0; char c = gc();36 while(c < '0' || c > '9') c = gc();37 while(c >='0' && c <='9') x = x * 10 + c - '0',c = gc();38 return x;39 }40 void FFT(C *a,int f){41 for(int i = 1;i < n;i++) if(i < rev[i]) swap(a[i],a[rev[i]]);42 for(int i = 1;i < n;i<<=1){43 C wn(cos(Pi/i),f * sin(Pi/i));44 for(int j = 0;j < n;j+=(i<<1)){45 C w(1,0);46 for(int k = 0;k < i;k++,w=w*wn){47 C x = a[j + k],y = w * a[j + k + i];48 a[j + k] = x + y,a[j + k + i] = x - y;;49 }50 }51 }52 if(f==-1) for(int i = 0;i < n;i++) a[i] =a[i]/n;53 }54 int main()55 { freopen("BZOJ3513.in","r",stdin);56 freopen("BZOJ3513.out","w",stdout);57 int T = read();58 while(T--){59 m = read();60 for(int i = 0;i <= 100000;i++) a[i] = 0;61 int mx = 0;62 for(int i = 1,t;i <= m;i++) a[t = read()]++,mx = max(t,mx);63 L = 0;for(n=1;n<=mx * 2;n<<=1) L++;64 for(int i = 1;i < n;i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(L-1));65 for(int i = 0;i < n;i++) b[i] = C(a[i],0);66 FFT(b,1); for(int i = 0;i < n;i++) b[i] = b[i]*b[i]; FFT(b,-1);67 ll ans = 0,sum = 0,all = 1ll * m * (m - 1) * (m - 2) / 6;68 for(int i = 1;i <= mx;i++){69 ll t = (ll)(b[i].x + 0.1);70 sum += t; if(!(i&1)) sum -= a[i/2];71 ans += sum * a[i];72 }73 ans /= 2;74 printf("%.7lf\n",1.0 * (all - ans) / all);75 }76 return 0;77 }//by tkys_Austin;
BZOJ 3160 万径人踪灭
1 /*FFT好短啊; 2 学到了新的思想构造多项式。 3 脑海中要浮现卷积的样子; 4 */ 5 #include6 #include 7 #include 8 #include 9 #include 10 #define N 26214411 #define mod 100000000712 #define ll long long13 using namespace std;14 typedef complex C; 15 const double pi = acos(-1);16 int n,m,L,len,rev[N],bin[N],p[N],f[N];17 char s[N/2],st[N];18 C a[N],b[N];19 int manacher(){20 st[0] = '$'; st[1] = '#';21 for(int i = 0;i < len;i++) st[i+1<<1] = s[i],st[(i+1<<1)|1] = '#';22 int mx = 0,ret = 0,pos = 0;23 for(int i = 0;i < 2 * len + 1;i++){24 if(mx>i) p[i] = min(p[2*pos-i],mx-i);else p[i] = 0;25 while(st[i+p[i]]==st[i-p[i]]) p[i]++;26 if(i+p[i]>mx) mx = i + p[i],pos = i;27 ret = (ret + p[i]/2) % mod;28 }29 return ret;30 }31 void FFT(C *a,int f){32 for(int i = 0;i < n;i++) if(i < rev[i]) swap(a[i],a[rev[i]]);33 for(int i = 1;i < n;i<<=1){34 C wn(cos(pi/i),f * sin(pi/i));35 for(int j = 0;j < n;j+=(i<<1)){36 C w(1,0);37 for(int k = 0;k < i;k++,w*=wn){38 C x = a[j + k],y = w * a[j + k + i];39 a[j + k] = x + y,a[j + k + i] = x - y; 40 }41 }42 }43 if(f==-1) for(int i = 0;i < n;i++) a[i]/=n;44 }45 int main()46 { freopen("BZOJ3160.in","r",stdin); freopen("BZOJ3160.out","w",stdout);47 scanf("%s",s); len = strlen(s);48 m = 2 * (len - 1); for(n=1;n<=m;n<<=1)L++;49 for(int i = 1;i < n;i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(L-1));50 for(int i = bin[0] = 1;i < n;i++) bin[i] = (bin[i - 1] <<1) % mod;51 for(int i = 0;i < len;i++) a[i] = (s[i]=='a');52 FFT(a,1); 53 for(int i = 0;i < n;i++) b[i] = a[i] * a[i];54 //FFT(b,-1);55 for(int i = 0;i < len;i++) a[i] = (s[i]=='b');56 for(int i = len;i < n;i++) a[i] = 0;57 FFT(a,1); 58 for(int i = 0;i < n;i++) b[i] += a[i] * a[i];59 FFT(b,-1);60 for(int i = 2;i < 2 * len + 1;i++)f[i] = (ll)(b[i - 2].real() + 0.1);61 int ans = 0;62 for(int i = 2;i < 2 * len + 1;i++)ans = (ans + bin[(f[i] + 1)>>1] - 1)% mod;63 ans = (ans - manacher() + mod)%mod;64 printf("%d\n",ans);65 return 0;66 }//by tkys_Austin;
BZOJ 4827 [Hnoi2017]礼物
1 /*那么今天亲眼目睹了精度问题; 2 少加了0.1就交了七遍,下次小心点,要是考试时遇到这种情况就惨了; 3 平时要注意; 4 虽然还是没想出来,但是想了个大半; 5 首先c可以O(1)计算; 6 然后得到新的a,b;化一化式子,用卷积求a_{i + k} * b{i} 同上一个(还有一个权限的此类裸题)*/ 7 #include8 #include 9 #include 10 #include 11 #include 12 #define N 1<<2013 #define inf 0x3f3f3f3f14 using namespace std;15 typedef complex C;16 const double Pi = acos(-1);17 int n,m,L,rev[N],s1[N],s2[N];18 C a[N],b[N];19 char gc(){20 static char *p1,*p2,s[1000000];21 if(p1==p2) p2 = (p1=s) + fread(s,1,1000000,stdin);22 return(p1==p2)?EOF:*p1++;23 }24 int read(){25 int x = 0,f = 1; char c = gc();26 while(c<'0'||c>'9') { if(c=='-') f = -1; c = gc();}27 while(c>='0'&&c<='9') x = x * 10 + c - '0',c = gc();28 return x * f; 29 }30 int sqr(int x){ return x * x;}31 void FFT(C *a,int f){32 for(int i = 1;i < n;i++) if(i < rev[i]) swap(a[i],a[rev[i]]);33 for(int i = 1;i < n;i<<=1){34 C wn(cos(Pi/i),f * sin(Pi/i));35 for(int j = 0;j < n;j += (i<<1)){36 C w(1,0);37 for(int k = 0;k < i;k++,w*=wn){38 C x = a[j + k],y = w * a[j + k + i];39 a[j + k] = x + y,a[j + k + i] = x - y;40 }41 }42 }43 if(f==-1) for(int i = 0;i < n;i++) a[i]/=n;44 }45 int main()46 { freopen("BZOJ4827.in","r",stdin);47 freopen("BZOJ4827.out","w",stdout);48 n = read(); m = read(); 49 double sum = 0;50 for(int i = 0,t;i < n;i++) sum -= (t = read()),s1[i] = t;51 for(int i = 0,t;i < n;i++) sum += (t = read()),s2[i] = t;52 int c = round(sum / n);53 long long ans = 0;54 for(int i = 0;i < n;i++) s1[i] += c,ans += sqr(s1[i]) + sqr(s2[i]);55 for(int i = 0;i < n;i++) a[i] = s1[i],b[i] = s2[n - 1 - i];56 for(int i = n;i < 2 * n;i++) a[i] = s1[i - n];57 m = 3 * n; for(n = 1;n<=m;n<<=1)L++;58 for(int i = 1;i < n;i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(L-1));59 FFT(a,1); FFT(b,1);60 for(int i = 0;i < n;i++) a[i] *= b[i];61 FFT(a,-1);62 int mx = 0; m /= 3;63 for(int i = m - 1;i < 2 * m - 1;i++) mx = max((int)(a[i].real() + 0.1),mx);64 ans -= 2 * mx;65 printf("%lld\n",ans);66 return 0;67 }//by tkys_Austin;
大米飘香的总结:
FFT的话本身附带了很多数学知识,感觉只要耐心其实很简单,代码也很短,主要是题中对式子的处理会有一些技巧。但是我竟然写了一天半。。。。。。开学了,还不知道能不能停课,但是我还是要争取一下,一学期过去了,我已经,长大了。。。。。。。。。。
至少有十年不曾留泪 至少有十首歌给我安慰
可现在我会莫名地心碎 可现在我会莫名地哭泣
当我想你的时候
------汪峰《当我想你的时候》