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+b^n\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$，
$$ga_n = \sum_{i=0}^{n}\left(\binom{n-1+i}{i}+\binom{n-i+1}{i-1}\right)a^i\\=ga_{n-1} +\binom{2n-1}{n}a^n+\sum_{i=0}^{n-1} \binom{n+i}{i}a^{i+1}\\=ga_{n-1}+\binom{2n-1}{n}a^{n}+a\cdot ga_n-\binom{2n}{n}a^{n+1}$$
移项，并两边同时除以$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;
}
