@KirinBill
2018-07-13T21:15:00.000000Z
字数 6414
阅读 1569
线性代数
目录
对于常系数齐次递推式,虽然我们可以用矩阵快速幂求了,但是如果递推式的阶数比较大,那就很麻烦了,我们这里介绍一种神奇的优化。
模板题:[BZOJ 4161] Shlw loves matrixI
朴素版:
#include <cstdio>
#include <cstring>
#define mset(a,x,n) memset(a,x,(n)*sizeof(*(a)))
#define mcpy(a,b,n) memcpy(a,b,(n)*sizeof(*(a)))
const int MAXN=2005,P=1e9+7;
int a[MAXN],b[MAXN];
int f[MAXN],g[MAXN],h[MAXN];
inline void inc(int &x,int y){x+=y,x= x>=P ? x-P:x;}
inline int pow(int x,int y){
int ret=1;
for(;y;x=(long long)x*x%P,y>>=1){
if(y&1) ret=(long long)ret*x%P;
}
return ret;
}
inline int inv(int x){return pow(x,P-2);}
inline void poly_mul(int a[],int b[],int p[],int m){
static int c[MAXN<<1];
mset(c,0,(m<<1)-1);
for(int i=0;i<m;++i){
for(int j=0;j<m;++j)
inc(c[i+j],(long long)a[i]*b[j]%P);
}
for(int i=(m<<1)-2;i>=m;--i){
int k=(long long)c[i]*inv(p[m])%P;
for(int j=0;j<=m;++j)
inc(c[i-j],P-(long long)k*p[m-j]%P);
}
mcpy(a,c,m);
}
inline void poly_pow(int x[],int y,int p[],int m,int ret[]){
ret[0]=1;
for(;y;poly_mul(x,x,p,m),y>>=1){
if(y&1) poly_mul(ret,x,p,m);
}
}
int main(){
int n,k;
scanf("%d%d",&n,&k);
for(int i=1;i<=k;++i){
scanf("%d",&a[i]);
a[i]=(a[i]+P)%P;
}
for(int i=0;i<k;++i){
scanf("%d",&b[i]);
b[i]=(b[i]+P)%P;
}
if(n<k){
printf("%d",b[n]);
return 0;
}
for(int i=1;i<=k;++i) f[k-i]=P-a[i];
f[k]=1;
g[1]=1;
poly_pow(g,n,f,k,h);
int ans=0;
for(int i=0;i<k;++i)
inc(ans,(long long)h[i]*b[i]%P);
printf("%d",ans);
return 0;
}
代码:
#include <cstdio>
#include <algorithm>
#include <cstring>
#define mset(a,x,n) memset(a,x,(n)*sizeof(*(a)))
#define mcpy(a,b,n) memcpy(a,b,(n)*sizeof(*(a)))
using std::min;
const int MAXN=1005,P=998244353;
int n,q,K;
int powq[MAXN];
int f[MAXN][MAXN],g[MAXN],h[MAXN<<1];
int a[MAXN],b[MAXN];
inline void inc(int &x,int y){x+=y,x= x>=P ? x-P:x;}
inline int pow(int x,int y){
int ret=1;
for(;y;x=(long long)x*x%P,y>>=1){
if(y&1) ret=(long long)ret*x%P;
}
return ret;
}
inline int inv(int x){return pow(x,P-2);}
inline void DP(int K){
mset(g,0,K+1);
g[0]=1;
for(int i=K;i;--i){
int m=min(n,K/i);
mset(f[i],0,m+1);
for(int j=1;j<=m;++j){
for(int k=1;k<=j;++k){
inc(f[i][j],(long long)g[k-1]*g[j-k]%P);
if(k!=K) inc(f[i][j],(long long)g[k-1]*f[i][j-k]%P);
}
f[i][j]=(long long)f[i][j]*(1-q+P)%P;
}
for(int j=1;j<=m;++j)
g[j]=(long long)(g[j]+f[i][j])%P*powq[j]%P;
}
}
inline void poly_mul(int a[],int b[],int ret[],int p[],int m){
static int c[MAXN<<1];
mset(c,0,(m<<1)-1);
for(int i=0;i<m;++i){
for(int j=0;j<m;++j)
inc(c[i+j],(long long)a[i]*b[j]%P);
}
for(int i=(m<<1)-1;i>=m;--i){
int k=(long long)c[i]*p[m]%P;
for(int j=0;j<=m;++j)
inc(c[i-j],P-(long long)k*p[m-j]%P);
}
mcpy(ret,c,m+1);
}
inline void poly_pow(int ret[],int y,int p[],int m){
static int x[MAXN];
mset(x,0,m+1),mset(ret,0,m);
x[1]=1,ret[0]=1;
for(;y;y>>=1){
if(y&1) poly_mul(ret,x,ret,p,m);
poly_mul(x,x,x,p,m);
}
}
inline int cal(int K){
if(K==0) return pow((1-q+P)%P,n);
DP(K);
mcpy(h,g,K+1);
for(int i=K+1;i;--i)
g[i]=(long long)(1-q+P)*g[i-1]%P;
for(int i=1;i<=K;++i){
for(int j=1;j<=i;++j)
inc(h[i],(long long)h[i-j]*g[j]%P);
}
if(n<=K) return h[n];
a[K+1]=1;
for(int i=K+1;i;--i)
a[K+1-i]=P-g[i];
poly_pow(b,n,a,K+1);
int ret=0;
for(int i=0;i<=K;++i)
inc(ret,(long long)b[i]*h[i]%P);
return ret;
}
inline void pre_cal(){
powq[0]=1;
for(int i=1;i<=K;++i)
powq[i]=(long long)powq[i-1]*q%P;
}
int main(){
int x,y;
scanf("%d%d%d%d",&n,&K,&x,&y);
q=(long long)x*inv(y)%P;
pre_cal();
printf("%d",(cal(K)-cal(K-1)+P)%P);
return 0;
}