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$$

  求$f_{n,n}\mod (10^6+3)$。

题解

 利用$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;
}
上一篇:pom.xml文件中,添加自定义参数Properties


下一篇:UVa 557 (概率 递推) Burger