hdu4914 Linear recursive sequence

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

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

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

用多项式乘法+快速幂计算f_n(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;
}