codechef Random Number Generator

题意:求常系数齐次递推数列的第n项

当递推系数非零项很多的时候暴力复位一次是O(k²)的

构造特征矩阵M的特征多项式f(x)=x^k-∑_i c_i*xⁱ

由f(M)=0,那么我们复位时只需要mod f(x)就好

Portal

#include <stdio.h>
#include <string.h>
#include <algorithm>

#define g() getchar()
template<class Q>void Scan(Q&x){
	char c; while(c=g(),c<48||c>57);
	for(x=0;c>47&&c<58;c=g())x=10*x+c-48;
}

#define G 3
#define maxn 65536
#define DFT(a,n) fft(a,n,0)
#define IDFT(a,n) fft(a,n,1)
#define INV(x) qpow(x,P-2)
#define rep(i,a,b) for(int i=a;i<(b);++i)

using namespace std;

const int MAXN=(1<<16)+5;
const int P=104857601;
typedef long long ll;

inline ll qpow(ll x,int k){
	ll ans=1;
	for(;k;x=x*x%P,k>>=1)
		if(k&1)ans=ans*x%P;
	return ans;
}

int r[MAXN],w[2][MAXN],rn;

inline void init_fft(int n,int b){	
	rn=INV(n);
	rep(i,1,n){
		r[i]=(r[i>>1]>>1)|((i&1)<<(b-1));
	}
}

inline void fft(int*a,int n,int t){
	rep(i,1,n)if(i<r[i])swap(a[i],a[r[i]]);
	for(int i=2;i<=n;i<<=1)
		for(int j=0;j<n;j+=i)
			for(int k=0;k<i>>1;++k){
				int re=1ll*a[(i>>1)+j+k]*w[t][maxn/i*k]%P;
				a[(i>>1)+j+k]=(a[j+k]+P-re)%P;
				a[j+k]=(a[j+k]+re)%P;
			}
	if(t)rep(i,0,n)a[i]=1ll*a[i]*rn%P;
}

int k,len,bit;

inline void inv(int*a,int*b){
	static int c[MAXN],l;
	b[0]=1,b[1]=0,l=1;
	rep(j,1,bit+1){
		l<<=1;
		rep(i,0,l)c[i]=a[i];
		rep(i,l,l<<1)b[i]=c[i]=0;
		init_fft(l<<1,j+1);
		
		DFT(b,l<<1),DFT(c,l<<1);
		rep(i,0,l<<1){
			b[i]=1ll*b[i]*(2+P-1ll*b[i]*c[i]%P)%P;
		}
		IDFT(b,l<<1);
		rep(i,l,l<<1)b[i]=0;
	} 
}

int m[MAXN],dm[MAXN],im[MAXN];

int ta[MAXN],ra[MAXN];

inline void mod(int*a){
	DFT(ta,len),DFT(a,len);
	rep(i,0,len)a[i]=1ll*a[i]*ta[i]%P;
	IDFT(a,len); 
	
	rep(i,0,k-1)ra[i]=a[(k-1<<1)-i];
	rep(i,k-1,len)ra[i]=0;
	DFT(ra,len);
	rep(i,0,len)ra[i]=1ll*ra[i]*im[i]%P;
	IDFT(ra,len);

	rep(i,0,k-1)ta[i]=ra[k-2-i];
	rep(i,k-1,len)ta[i]=0;
	DFT(ta,len);
	rep(i,0,len)ta[i]=1ll*ta[i]*dm[i]%P;
	IDFT(ta,len);

	rep(i,0,k)a[i]=(a[i]+P-ta[i])%P;
	rep(i,k,len)a[i]=0;
}

inline void Init(){
	int n=maxn;
	int g=qpow(G,(P-1)/n);
	w[0][0]=w[0][n]=1;
	rep(i,1,n)w[0][i]=1ll*w[0][i-1]*g%P;
	rep(i,0,n+1)w[1][i]=w[0][n-i];
}

ll n;

int a[MAXN],ans[MAXN],x[MAXN];

int main(){
	Init();
	
	Scan(k),Scan(n),--n;
	
	len=k+1,bit=1;
	while(1<<bit<len)++bit;
	len=1<<bit;
	
	rep(i,0,k)Scan(a[i]);
	for(int i=k-1;~i;--i){
		Scan(m[i]),m[i]=(P-m[i])%P;
	}
	
	if(k==1){
		int res=qpow((P-m[0])%P,n%(P-1))*a[0]%P;
		printf("%d\n",res);
		return 0;
	} 
	
	m[k]=1;
	rep(i,0,k+1)dm[i]=m[k-i];
	inv(dm,im);
	
	len<<=1,++bit;
	rep(i,0,k+1)dm[i]=m[i];
	rep(i,k+1,len)dm[i]=im[i]=0;
	DFT(dm,len),DFT(im,len);
	
	for(ans[0]=x[1]=1;n;n>>=1){
		if(n&1){
			rep(i,0,len)ta[i]=x[i]; 
			mod(ans);
		}
		rep(i,0,len)ta[i]=x[i];
		mod(x);
	}
	
	int res=0;
	rep(i,0,k){
		res=(res+1ll*a[i]*ans[i])%P;
	}
	printf("%d\n",res);
	
	return 0;
}