ShinriiTin's Blog - 博主是sb

codechef Random Number Generator

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

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

构造特征矩阵M的特征多项式f(x)=xk-∑i ci*xi

由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;
}

未分类 FFT 常系数齐次递推 Comments(0) 2015年6月12日 17:41

hdu4914 Linear recursive sequence

题意:给出序列f,fn=(n>0)?a*fn-p+b*fn-q:1,求第n项

叉姐论文: k阶矩阵M的任意幂次都能用M小于k次的幂次来线性表示

所以我们构造多项式fn(x),[xi]fn(x)表示用小于k次的幂次线性表示时i次幂前的系数

用多项式乘法+快速幂计算fn(x),每次乘法后溢出k-1次的部分暴力还原

这题 mod 119，中间结果不会太大，所以直接用复数域上的fft，做完一次乘法后再取模就好

Portal

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#define DFT(a) fft(a,len,0)
#define IDFT(a) fft(a,len,1)
#define rep(i,a,b) for(int i=a;i<(b);++i)
using namespace std;

const int MAXN=(1<<18)+5; 

typedef double db;
const db Pi=acos(-1.);

typedef long long ll;
const int Mod=119;

struct cp{
	db r,i;
	cp(db r=0,db i=0):r(r),i(i){}
	friend cp operator+(const cp&a,const cp&b){
		return cp(a.r+b.r,a.i+b.i);
	}
	friend cp operator-(const cp&a,const cp&b){
		return cp(a.r-b.r,a.i-b.i);
	}
	friend cp operator*(const cp&a,const cp&b){
		return cp(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);
	}
}w[2][MAXN];

inline cp Exp(db x){
	return cp(cos(x),sin(x));
}

int r[MAXN];

inline void init_fft(int n,int b){
	w[0][0]=w[0][n]=cp(1.,0);
	cp g=Exp(Pi*2./n);
	rep(i,1,n)w[0][i]=w[0][i-1]*g;
	rep(i,0,n+1)w[1][i]=w[0][n-i];
	rep(i,1,n){
		r[i]=(r[i>>1]>>1)|((i&1)<<(b-1));
	}
} 

inline void fft(cp*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){
				cp re=a[(i>>1)+j+k]*w[t][n/i*k];
				a[(i>>1)+j+k]=a[j+k]-re;
				a[j+k]=a[j+k]+re;
			} 
	if(t)rep(i,0,n)a[i].r/=n;
}

cp ta[MAXN],tb[MAXN];

int n,A,B,p,q,len,bit;

struct poly{
	int a[MAXN];
	inline void mul(const poly&o){
		rep(i,0,q)ta[i]=cp(a[i]),tb[i]=cp(o.a[i]);
		rep(i,q,len)ta[i]=tb[i]=cp();
		DFT(ta),DFT(tb);
		rep(i,0,len)ta[i]=ta[i]*tb[i];
		IDFT(ta);
		rep(i,0,len)a[i]=ll(ta[i].r+0.5)%Mod;
		for(int i=len;i>=q;--i)if(a[i]){
			a[i-p]=(a[i-p]+a[i]*A%Mod)%Mod;
			a[i-q]=(a[i-q]+a[i]*B%Mod)%Mod;
		}
	}
}a,b;

inline void qpow(poly&ans,poly&x,int k){
	rep(i,0,len)ans.a[i]=x.a[i]=0;
	ans.a[0]=x.a[1]=1;
	for(;k;x.mul(x),k>>=1)
		if(k&1)ans.mul(x);
}

int f[MAXN];

inline int F(int x){return x>0?f[x]:1;}

int main(){
	while(~scanf("%d%d%d%d%d",&n,&A,&B,&p,&q)){
		A%=Mod,B%=Mod;
		rep(i,1,q<<1|1){
			f[i]=(F(i-p)*A+F(i-q)*B)%Mod; 
		}
		if(n<=(q<<1)){
			printf("%d\n",f[n]);
			continue;
		}
		len=q<<1|1,bit=1;
		while(1<<bit<len)++bit;
		len=1<<bit;
		init_fft(len,bit);
		qpow(a,b,n-q);
		int ans=0;
		for(int i=q-1;~i;--i){
			ans=(ans+a.a[i]*f[i+q])%Mod;
		}
		printf("%d\n",ans);
	}
	return 0;
}

未分类 FFT 常系数齐次递推 Comments(0) 2015年6月12日 17:24

bzoj3456 城市规划

求n个点的简单无向连通图(无标号)的个数，对1004535809取模

令fn表示n个点的方案数，则有递推公式: fn=2C(n,2)-∑i(i>1&&i<n) C(n-1,i-1)*fi*2C(n-i,2) (用总方案数减去不连通的方案数)

两边都乘上(n-1)!-1然后再移项整理得:2C(n,2)*(n-1)-1=∑i(i>1&&i<=n)fi*(i-1)!*2C(n-i,2)*(n-i)!

令Ai=fi*(i-1)!,Bi=2C(i,2)*i!,Ci=2C(n,2)*(n-1)!-1,则有A×B=C,A=C*B-1

这样就只用算出B在mod xn+1意义下的逆元，就可以计算A了。

Portal

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define G 3
#define DFT(a,n) fft(a,n,0)
#define IDFT(a,n) fft(a,n,1)
#define INV(x,P) qpow(x,P-2,P)
#define rep(i,a,b) for(int i=a,s=b;i<s;++i)
using namespace std;

const int MAXN=(1<<18)+5;
typedef long long ll;
const int P=1004535809;

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

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

inline void init_fft(int n,int b){
	rn=INV(n,P);
	int g=qpow(G,(P-1)/n,P);
	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];
	rep(i,1,n+1){
		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 v=1ll*a[(i>>1)+j+k]*w[t][n/i*k]%P;
				a[(i>>1)+j+k]=(a[j+k]+P-v)%P;
				a[j+k]=(a[j+k]+v)%P;
			}
	if(t)rep(i,0,n)a[i]=1ll*a[i]*rn%P;
}

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

int n,l,b;

int fac[MAXN],inv_fac[MAXN];

inline void Init(){
	scanf("%d",&n);
	l=n+1,b=1;
	while(1<<b<l)++b;
	l=1<<b;
	fac[0]=1;
	rep(i,1,l)fac[i]=1ll*fac[i-1]*i%P;
	inv_fac[0]=inv_fac[1]=1;
	rep(i,2,l){
		int a=P/i,b=P%i;
		inv_fac[i]=1ll*a*(P-inv_fac[b])%P;
	}
	rep(i,2,l)
		inv_fac[i]=1ll*inv_fac[i]*inv_fac[i-1]%P;
}

inline ll C(ll n){return n*(n-1)>>1;}

int a[MAXN],c[MAXN];

inline void Solve(){
	rep(i,0,n+1)
		c[i]=qpow(2,C(i)%(P-1),P)*inv_fac[i]%P;
	inv(c,a,l,b);
	c[0]=0;
	rep(i,1,n+1)
		c[i]=qpow(2,C(i)%(P-1),P)*inv_fac[i-1]%P;
	rep(i,l,l<<1)a[i]=c[i]=0;
	init_fft(l<<1,b+1);
	DFT(a,l<<1),DFT(c,l<<1);
	rep(i,0,l<<1)a[i]=1ll*a[i]*c[i]%P;
	IDFT(a,l<<1);
	int ans=1ll*a[n]*fac[n-1]%P;
	printf("%d\n",ans);
}

int main(){
	Init(); Solve(); return 0;
}

未分类 FFT 多项式求逆 Comments(0) 2015年6月11日 15:27