BZOJ4451 [Cerc2015]Frightful Formula 多项式 FFT 递推 组合数学
原文链接http://www.cnblogs.com/zhouzhendong/p/8820963.html
题目传送门 - BZOJ4451
题意
给你一个$n\times n$矩阵的第一行和第一列，其余的数通过如下公式推出：
$$f_{i,j}=a\cdot f_{i,j-1}+b\cdot f_{i-1,j}+c$$
题解
利用$FFT$来解决此题
真是一道好题只是我太菜了。光$FFT$就调了好久。当然本题可以直接递推（将写在用$FFT$实现的算法的后面），用$FFT$只是用来给脑子偷懒的。
这题的主要思路很巧妙。
先假设$a=b=1$。
我们考虑对于每一个$f_{i,j}$计算它对$f_{n,n}$的贡献次数。
显然贡献次数就是在网格图中从$(i,j)$开始只能往右或者往下走的路径条数，即$\binom{n-i+n-j}{n-i}$。
注意第一行和第一列的对答案的贡献式有点不同。
因为第一行的格子不能对右边的格子做贡献，所以：
第一行的格子$(1,i)$的贡献次数就是$\binom{2n-i-2}{n-i}$。
同理，
第一列的格子$(i,1)$的贡献次数就是$\binom{2n-i-2}{n-i}$。
现在不考虑$a=b=1$。
我们还要乘一个权值。考虑从第$i$行走到第$n$行向下走了$n-i$次，所以乘了$n-i$次$b$，即乘了$b^{n-i}$。同理向右走要乘的权值就是$a^{n-j}$。所以综上所述，我们再重新写一下贡献次数的式子。
下面的三式满足$i,j>1$。注意$(1,1)$是没有贡献次数的。（一开始没注意到，被坑了很久）
$(1,i):\binom{2n-i-2}{n-i}a^{n-1}b^{n-i}$
$(i,1):\binom{2n-i-2}{n-i}a^{n-i}b^{n-1}$
$(i,j):\binom{2n-i-j}{n-i}a^{n-i}b^{n-j}$
考虑到$(i,j)$这种格子的值里面含有$(i,1)$和$(1,i)$这种的贡献，贡献不能重复计算，所以我们在计算格子$(i,j)$的时候就只考虑它在被贡献的时候产生的$c$对答案的贡献次数。
（就是上面的那个）
所以，我们可以列出答案的式子：
$$ans=\sum_{i=2}^{n}f_{i,1}\cdot\binom{2n-i-2}{n-i}a^{n-i}b^{n-1}\\+\sum_{i=2}^{n}f_{1,i}\cdot\binom{2n-i-2}{n-i}a^{n-i}b^{n-1}\\+\sum_{i=2}^{n}\sum_{j=2}^{n}c\cdot\binom{2n-i-j}{n-i}a^{n-i}b^{n-j}$$
各种逆元、阶乘以及幂的预处理就不说了。
于是前面的两个式子都可以$O(n)$搞定。
后面的那个式子，我们在给他稍稍变几个形：
$$\sum_{i=2}^{n}\sum_{j=2}^{n}c\cdot\binom{2n-i-j}{n-i}a^{n-i}b^{n-j}$$
把$c$提前，并修改求和指标，枚举$n-i$和$n-j$，得：
$$=c\sum_{i=0}^{n-2}\sum_{j=0}^{n-2}\frac{(i+j)!}{i!j!}a^ib^j$$
$$=c\sum_{i=0}^{n-2}i!\sum_{j=0}^{n-2}\frac{a^{(i-j)}}{(i-j)!}\cdot\frac{b^j}{j!}$$
设
$$f_i=\begin{cases}\large{\frac{a^i}{i!}}& \text{$(0\leq i\leq n-2)$}\\0& \text{$(n-2<i)$}\end{cases}$$
$$g_i=\begin{cases}\large{\frac{b^i}{i!}}& \text{$(0\leq i\leq n-2)$}\\0& \text{$(n-2<i)$}\end{cases}$$
则原式可以写为：
$$=c\sum_{i=0}^{n-2}(i+j)!\sum_{j=0}^{n-2}f_ig_j$$
发现这个符合多项式卷积形式。于是我们再换一个写法：
$$=c\sum_{i=0}^{2n-4}i!\sum_{j=0}^{i}f_jg_{i-j}$$
于是显然可以$FFT$快速计算了。注意，要取模，而且卷积得到的数字较大，可以用任意模数$FFT$或者拆系数$FFT$。我写的是拆系数。因为这样好写。
时间复杂度$O(n\log n)$。
但是这样不仅跑的慢，写起来也麻烦。
这题其实可以把$n$搞到$3000000$级别，把模数搞到$1e9+7$。
因为存在线性递推方式。
线性递推
根本没想到还有这样的操作。可能是本着练习$FFT$的目的然后就没仔细想了。
可以线性递推。
我们把之前的式子继续推下去。
$$c\sum_{i=0}^{n-2}\sum_{j=0}^{n-2}\frac{(i+j)!}{i!j!}a^ib^j$$
我们只需要考虑后面的式子，先不看$c$。
$$\sum_{i=0}^{n-2}\sum_{j=0}^{n-2}\binom{i+j}{i}a^ib^j$$
设
$$f_n=\sum_{i=0}^{n}\sum_{j=0}^{n}\binom{i+j}{i}a^ib^j$$
则我们要求的就是$f_{n-2}$。
推导：
$$\begin{eqnarray*}f_n&=&\sum_{i=0}^{n}\sum_{j=0}^{n}\binom{i+j}{i}a^ib^j\\&=&\sum_{i=0}^{n}a^i\sum_{j=0}^{n}\binom{i+j}{i}b^j\\&=&\sum_{i=0}^{n-1}a^i\sum_{j=0}^{n-1}\binom{i+j}{i}b^j+a^n\sum_{j=0}^{n}\binom{n+j}{n}b^j+\sum_{i=0}^{n}a^i\binom{i+n}{i}-\binom{2n}{n}a^nb^n\\&=&f_{n-1}+a^n\sum_{i=0}^{n}\binom{n+i}{n}b^i+b^n\sum_{i=0}^{n}\binom{n+i}{n}a^i-\binom{2n}{n}a^nb^n\end{eqnarray*}$$
我们只需要处理$\sum_{i=0}^{n}\binom{n+i}{n}b^i$和$\sum_{i=0}^{n}\binom{n+i}{n}a^i$，就可以搞定了。
由于这两个形式一样，我这里只对$a$进行推导。
设
$$ga_n=\sum_{i=0}^{n}\binom{n+i}{n}a^i$$
则
$$f_n=f_{n-1}+a^ngb_n+b^nga_n-\binom{2n}{n}a^nb^n$$
推一下 $ga_n$，
$$\begin{array} { c } { ga _ { n } = \sum _ { i = 0 } ^ { n - 1 } \left( \left( \begin{array} { c } { n - 1 + i } \\ { i } \end{array} \right) + \left( \begin{array} { c } { n - 1 + i } \\ { i - 1 } \end{array} \right) \right) a ^ { i } } \\ { = ga _ { n - 1 } + \left( \begin{array} { c } { 2 n - 2 } \\ { n - 1 } \end{array} \right) a ^ { n - 1 } + \sum _ { i = 0 } ^ { n - 2 } \left( \begin{array} { c } { n + i } \\ { i } \end{array} \right) a ^ { i + 1 } } \\ { = ga _ { n - 1 } + \left( \begin{array} { c } { 2 n - 2 } \\ { n - 1 } \end{array} \right) a ^ { n - 1 } + a ga _ { n } - \left( \begin{array} { c } { 2 n - 1 } \\ { n - 1 } \end{array} \right) a ^ { n } } \end{array}$$
移项，并两边同时除以$1-a$：
$$\begin{eqnarray*}&ga_n&=&ga_{n-1}+\binom{2n-1}{n}a^n+a\cdot ga_n-\binom{2n}{n}a^{n+1}\\ \Longrightarrow& (1-a)ga_n&=&ga_{n-1}+\binom{2n-1}{n}a^n-\binom{2n}{n}a^{n+1}\\ \Longrightarrow& ga_n&=&\frac{ga_{n-1}+\binom{2n-1}{n}a^n-\binom{2n}{n}a^{n+1}}{1-a}\end{eqnarray*}$$
然后就搞定了$ga$的递推，$gb$同理。
但是有特殊情况。
当$a=1$时，$1-a=0$，分母为$0$，这个公式挂掉了。
于是我们要解决这个问题。
方法是直接把$a=1$代入移项后的式子里面：
$$\begin{eqnarray*}&(1-a)ga_n&=&ga_{n-1}+\binom{2n-1}{n}a^n-\binom{2n}{n}a^{n+1}\\ \Longrightarrow& ga_{n-1}&=&\binom{2n}{n}-\binom{2n-1}{n}=\binom{2n-1}{n-1}\\ \Longrightarrow& ga_{n}&=&\binom{2n+1}{n}\end{eqnarray*}$$
就可以直接算$ga_n$了。
于是，
$$f_n=f_{n-1}+a^ngb_n+b^nga_n-\binom{2n}{n}a^nb^n$$
综合上面的两种递推式，写三支递推就可以搞定本题了。
而且一不留神卡常到第一了QAQ。
| No. | RunID | User | Memory | Time | Language | Code_Length | Submit_Time |
| 1 | 2705999 | zhouzhendong | 24728 KB | 580 MS | C++ | 1649 B | 2018-04-13 19:12:05 |
代码
代码1 - FFT实现
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1<<19;
double PI=acos(-1.0);
LL mod=1000003;
LL n,a,b,c,ans=0;
LL Fac[N],Inv[N],Powa[N],Powb[N];
LL fa[N],fb[N],tot[N];
int s,L,R[N];
struct C{double r,i;C(){}C(double a,double b){r=a,i=b;}C operator + (C x){return C(r+x.r,i+x.i);}C operator - (C x){return C(r-x.r,i-x.i);}C operator * (C x){return C(r*x.r-i*x.i,r*x.i+i*x.r);}
}w[N],A[N],B[N];
LL Pow(LL x,LL y){if (y==0)return 1LL;LL xx=Pow(x,y/2);xx=xx*xx%mod;if (y&1LL)xx=xx*x%mod;return xx;
}
LL inv(LL x){return Pow(x,mod-2);
}
LL func_C(int n,int m){return Fac[n]*Inv[m]%mod*Inv[n-m]%mod;
}
void FFT(C a[],int n){for (int i=0;i<n;i++)if (i<R[i])swap(a[i],a[R[i]]);for (int t=n>>1,d=1;d<n;d<<=1,t>>=1)for (int i=0;i<n;i+=(d<<1))for (int j=0;j<d;j++){C tmp=w[t*j]*a[i+j+d];a[i+j+d]=a[i+j]-tmp;a[i+j]=a[i+j]+tmp;}
}
void FFT_mul(C A[],C B[]){FFT(A,s),FFT(B,s);for (int i=0;i<s;i++)A[i]=A[i]*B[i],w[i].i*=-1.0;FFT(A,s);
}
int main(){scanf("%lld%lld%lld%lld",&n,&a,&b,&c);Fac[0]=Powa[0]=Powb[0]=1;for (int i=1;i<=n*2;i++){Fac[i]=Fac[i-1]*i%mod;Powa[i]=Powa[i-1]*a%mod;Powb[i]=Powb[i-1]*b%mod;}Inv[n*2]=inv(Fac[n*2]);for (int i=n*2-1;i>=0;i--)Inv[i]=Inv[i+1]*(i+1)%mod;for (int i=1,x;i<=n;i++){scanf("%d",&x);if (i>1)ans=(ans+func_C(n+n-i-2,n-i)*Powa[n-1]%mod*Powb[n-i]%mod*x)%mod;}for (int i=1,x;i<=n;i++){scanf("%d",&x);if (i>1)ans=(ans+func_C(n+n-i-2,n-i)*Powa[n-i]%mod*Powb[n-1]%mod*x)%mod;}for (int i=0;i<=n-2;i++)fa[i]=Powa[i]*Inv[i]%mod,fb[i]=Powb[i]*Inv[i]%mod;for (s=1,L=0;s<=2*n-4;s<<=1,L++);for (int i=0;i<s;i++){R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));w[i]=C(cos(2*i*PI/s),sin(2*i*PI/s));}LL ansdc=0;for (int i=0;i<s;i++)A[i]=C(fa[i]&1023,0),B[i]=C(fb[i]&1023,0);FFT_mul(A,B);for (int i=0;i<s;i++){LL v=((LL)(A[i].r/s+0.5))%mod;ansdc=(ansdc+Fac[i]*v)%mod;tot[i]=v,w[i].i*=-1.0;}for (int i=0;i<s;i++)A[i]=C(fa[i]>>10,0),B[i]=C(fb[i]>>10,0);FFT_mul(A,B);for (int i=0;i<s;i++){LL v=((LL)(A[i].r/s+0.5))%mod;ansdc=(ansdc+Fac[i]*48573%mod*v)%mod;tot[i]=(tot[i]+v)%mod,w[i].i*=-1.0;}for (int i=0;i<s;i++){A[i]=C((fa[i]>>10)+(fa[i]&1023),0);B[i]=C((fb[i]>>10)+(fb[i]&1023),0);}FFT_mul(A,B);for (int i=0;i<s;i++){LL v=((LL)(A[i].r/s+0.5))%mod;ansdc=(ansdc+Fac[i]*1024%mod*(v-tot[i]+mod))%mod;w[i].i*=-1.0;}ans=(ans+ansdc*c)%mod;printf("%lld",ans);return 0;
}
代码2 - 递推（没卡常）（代码3卡了一点常数的（跑的挺快））
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=200005,mod=1000003;
int read(){int x=0;char ch=getchar();while (!('0'<=ch&&ch<='9'))ch=getchar();while ('0'<=ch&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();return x;
}
LL n,a,b,c,ans=0;
LL Fac[mod],Inv[mod],Powa[N],Powb[N];
LL f[N],ga[N],gb[N];
LL Pow(LL x,LL y){if (y==0)return 1LL;LL xx=Pow(x,y/2);xx=xx*xx%mod;if (y&1LL)xx=xx*x%mod;return xx;
}
LL inv(int x){return Pow((x%mod+mod)%mod,mod-2);
}
LL C(int n,int m){return Fac[n]*Inv[m]%mod*Inv[n-m]%mod;
}
int main(){scanf("%lld%lld%lld%lld",&n,&a,&b,&c);Fac[0]=Powa[0]=Powb[0]=1;for (int i=1;i<=n;i++){Powa[i]=Powa[i-1]*a%mod;Powb[i]=Powb[i-1]*b%mod;}for (int i=1;i<=2*n;i++)Fac[i]=Fac[i-1]*i%mod;Inv[2*n]=Pow(Fac[2*n],mod-2);for (int i=2*n-1;i>=0;i--)Inv[i]=Inv[i+1]*(i+1)%mod;for (int x=read(),i=2;i<=n;i++)ans=(ans+C(n+n-i-2,n-i)*Powa[n-1]%mod*Powb[n-i]%mod*read())%mod;for (int x=read(),i=2;i<=n;i++)ans=(ans+C(n+n-i-2,n-i)*Powa[n-i]%mod*Powb[n-1]%mod*read())%mod;f[0]=ga[0]=gb[0]=1;LL inva=inv(1-a),invb=inv(1-b);for (int i=1;i<=n-2;i++){ga[i]=a==1?C(2*i+1,i):((ga[i-1]+C(2*i-1,i)*Powa[i]-C(2*i,i)*Powa[i+1])%mod*inva%mod);gb[i]=b==1?C(2*i+1,i):((gb[i-1]+C(2*i-1,i)*Powb[i]-C(2*i,i)*Powb[i+1])%mod*invb%mod);f[i]=(f[i-1]+Powa[i]*gb[i]+Powb[i]*ga[i]-C(2*i,i)*Powa[i]%mod*Powb[i])%mod;}printf("%lld",((ans+f[n-2]*c)%mod+mod)%mod);return 0;
}
代码3 - 卡常(580MS(截至2018-04-13代码运行速度Rk1))
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=200005,mod=1000003;
int read(){int x=0;char ch=getchar();while (!('0'<=ch&&ch<='9'))ch=getchar();while ('0'<=ch&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();return x;
}
LL n,a,b,c,ans=0;
LL Fac[mod],Inv[mod],Powa[N],Powb[N];
LL f[N],ga[N],gb[N];
LL Pow(LL x,LL y){if (y==0)return 1LL;LL xx=Pow(x,y/2);xx=xx*xx%mod;if (y&1LL)xx=xx*x%mod;return xx;
}
LL inv(int x){return Pow((x%mod+mod)%mod,mod-2);
}
LL C(int n,int m){return Fac[n]*Inv[m]%mod*Inv[n-m]%mod;
}
int main(){scanf("%lld%lld%lld%lld",&n,&a,&b,&c);Fac[0]=Powa[0]=Powb[0]=1;for (int i=1;i<=n;i++){Powa[i]=Powa[i-1]*a%mod;Powb[i]=Powb[i-1]*b%mod;}for (int i=1;i<=2*n;i++)Fac[i]=Fac[i-1]*i%mod;Inv[2*n]=Pow(Fac[2*n],mod-2);for (int i=2*n-1;i>=0;i--)Inv[i]=Inv[i+1]*(i+1)%mod;for (int x=read(),i=2;i<=n;i++)ans=(ans+C(n+n-i-2,n-i)*Powa[n-1]%mod*Powb[n-i]%mod*read())%mod;for (int x=read(),i=2;i<=n;i++)ans=(ans+C(n+n-i-2,n-i)*Powa[n-i]%mod*Powb[n-1]%mod*read())%mod;f[0]=ga[0]=gb[0]=1;LL inva=inv(1-a),invb=inv(1-b);for (int i=1;i<=n-2;i++){ga[i]=a==1?C(2*i+1,i):((ga[i-1]+C(2*i-1,i)*Powa[i]-C(2*i,i)*Powa[i+1])%mod*inva%mod);gb[i]=b==1?C(2*i+1,i):((gb[i-1]+C(2*i-1,i)*Powb[i]-C(2*i,i)*Powb[i+1])%mod*invb%mod);f[i]=(f[i-1]+Powa[i]*gb[i]+Powb[i]*ga[i]-C(2*i,i)*Powa[i]%mod*Powb[i])%mod;}printf("%lld",((ans+f[n-2]*c)%mod+mod)%mod);return 0;
}
转载于:https://www.cnblogs.com/zhouzhendong/p/BZOJ4451.html
BZOJ4451 [Cerc2015]Frightful Formula 多项式 FFT 递推 组合数学相关推荐
- BZOJ4451 : [Cerc2015]Frightful Formula
$(i,1)$对答案的贡献为$l_iC(2n-i-2,n-i)a^{n-1}b^{n-i}$. $(1,i)$对答案的贡献为$t_iC(2n-i-2,n-i)*a^{n-i}b^{n-1}$. $(i ...
- bzoj 4451 : [Cerc2015]Frightful Formula FFT
4451: [Cerc2015]Frightful Formula Time Limit: 10 Sec Memory Limit: 64 MB Submit: 177 Solved: 57 [S ...
- [FFT || 递推] BZOJ 4451 [Cerc2015]Frightful Formula
f1,i f_{1,i} 的贡献是 f1,i∗an−i∗bn−1∗Cn−i2n−i−2 f_{1,i}*a^{n-i}*b^{n-1}*C_{2n-i-2}^{n-i} fi,1 f_{i,1} 同理 ...
- 【bzoj 4451】[Cerc2015]Frightful Formula - 递推
才没有在做cerc2015呢 看到好像不少人这题写fft卡得死死的啊,不如 O(n) O(n) 递推(雾) 首先可以观察出 (i,1) (i,1)这个格子为 x x时对(n,n)(n,n) a,b a ...
- BZOJ 4451: [Cerc2015]Frightful Formula
Description 给你一个n*n矩阵的第一行和第一列,其余的数通过如下公式推出: \(f_{i,j}=a\times f_{i,j-1}+b\times f_{i-1,j}+c\) 求\(f_ ...
- bzoj 4451: [Cerc2015]Frightful Formula 数学+排列组合
题意 给你一个n*n矩阵的第一行和第一列,其余的数通过如下公式推出: F[i,j]=a*f[i,j-1]+b*f[i-1,j]+c 求f[n][n]%(10^6+3) 2<=n<=2000 ...
- CERC2015 Frightful Formula 神奇的模意义下分数
上网查了下这道题的正解是FFT......然而机智的xushu用待定系数法A了.... f[i,j]+k=a(f[i][j-1]+k)+b(f[i-1,j]+k) 解得 k=c/(a+b-1) 于是令 ...
- 【组合数学】递推方程 ( 常系数线性非齐次递推方程 的 非齐次部分是 多项式 与 指数 组合方式 | 通解的四种情况 )
文章目录 一.常系数线性非齐次递推方程 的 非齐次部分是 多项式 与 指数 组合方式 二.递推方程通解的四种情况 一.常系数线性非齐次递推方程 的 非齐次部分是 多项式 与 指数 组合方式 如果 &q ...
- matlab如何仿真递推型dft算法,递推dft算法
0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 伯格(Burg)递推算法 L-D算法缺点: 在计算相关函数估计时,对N个观测数据以 外的数据作零的假设,故谱估计误差较...... 第3 ...
最新文章
- javascript篇-console.log()打印object却显示为字符串[object object]
- Spring 详解（三）：AOP 面向切面的编程
- Ubuntu 16.04安装unrar解压RAR文件
- CocoaPods 错误解决 Attempt to read non existent folder
- mxnet基础到提高(21）－配置mxnet并运行第一个C＋＋程序
- Spring框架—SpringBean源码分析
- c 和java通讯大小端问题处理_记录一个如何解决java与C++socket通信的大小端问题...
- 计算机组成原理mbps,2016年湖北师范学院计算机组成原理(同等学力加试)复试笔试仿真模拟题...
- 【干货分享】流程DEMO-补打卡
- windows os x linux比较,对比测试：Ubuntu 11.04 vs Win7 vs OS X 10.7
- flash一个按钮控制动画_flutter闪屏过渡动画，闪光占位动画
- Volley源码学习1--volley结构图
- 使用JasperReport+iReport进行Web报表开发
- M1芯片CAD如何安装？M1 mac怎么安装AutoCAD？
- 机器学习及其MATLAB实现——BP神经网络
- Java基础知识总结（绝对经典）
- C++公历农历转换2020-2080年/除夕修正
- 为公司添加以网站作为邮箱后缀的企业邮箱
- GEE开发之Modis_LST地表温度数据分析
- OpenFOAM当中监测力和阻力系数
热门文章
- Windows Server 2003自带NAT功能,轻松实现不同网段互访
- matlab figure 怎么用,Matlab:figure的用法
- node.js安装模式 的区别_安装 若依 前后端 分离版
- 如何区分图像/物体“识别”和“分类”以及“检测”
- android 选择答题功能,Android APP编写简单答题器
- 【ROS】—— 机器人系统仿真 —Rviz中控制机器人模型运动与URDF集成Gazebo(十五)
- html编辑器自定义脚本,ckeditor自定义插件使用方法详解
- matlab的quiver函数的用法
- 设计模式 适配器模式 以手机充电器为例
- TTEFS-基于LayerFsd的文档透明加密SDK