Codeforces 392C. Yet Another Number Sequence 题解

题目链接：C. Yet Another Number Sequence

题目大意：令 (F_1=1,F_2=2,F_n=F_{n-1}+F_{n-2}(nge 3),A_i(k)=F_i imes i^k) ，求 (sum_{i=1}^{n} A_i(k)) 。

(nleq 10^{17},kleq 40)

Bonus: (kleq 200)

---

题解：我们令 (f(n,k)=sum_{i=1}^{n} F_i imes i^k) ，接下来划一波式子：

为了方便起见，令 (F_0=1,F_{-1}=0)

[f(n,j)=sum_{i=1}^{n} (F_{i-1}+F_{i-2}) imes i^k \ =sum_{i=1}^{n} F_{i-1} imes i^k +sum_{i=1}^{n} F_{i-2} imes i^k \ =sum_{i=1}^{n} F_{i-1}sum_{j=0}^{k} inom{k}{j} (i-1)^j+sum_{i=1}^{n} F_{i-2} sum_{j=0}^{k} inom{k}{j} (i-2)^j 2^{k-j} \ =sum_{j=0}^{k} inom{k}{j}sum_{i=0}^{n-1} F_{i} imes i^j + 1 +sum_{j=0}^{k}inom{k}{j} 2^{k-j}sum_{i=1}^{n-2}F_i imes i^j +2^k \ =sum_{j=0}^{k} inom{k}{j} f(n-1,j)+sum_{j=0}^{k} inom{k}{j}2^{k-j}f(n-2,j) +2^k+1]

划到这一步我们就可以直接矩阵快速幂了，矩阵大小为 ((2k+3) imes (2k+3)) ，时间复杂度 (O(k^3log n)) 。

接下来我们考虑 Bonus 中的数据范围怎么做。

我们令 (G) 为斐波那契数列的转移矩阵， (f(n,k)=sum_{i=1}^{n} G^{i} imes i^k)。

那么我们考虑从 (f(n,k)) 拓展到 (f(2n,k)) ：

[f(2n,k)=sum_{i=0}^{2n}G^{i} imes i^k \ =sum_{i=0}^{n}G^i imes i^k + G^nsum_{i=0}^{n}G^i imes (i+n)^k \ =f(n,k)+G^nsum_{i=0}^{n}G^isum_{j=0}^{k}inom{k}{j}i^jn^{k-j} \ =f(n,k)+G^nsum_{j=0}^{k}inom{k}{j}n^{k-j}sum_{i=0}^{n}G^ii^j \ =f(n,k)+G^nsum_{j=0}^{k}inom{k}{j}n^{k-j}f(n,j)]

直接递归下去处理即可，时间复杂度 (O(k^2log n + 8klog n)) 。

接下来是代码：

(O(k^3log n))

#include <cstdio>
typedef long long ll;
const int Maxk=40;
const int Mod=1000000007;
int C[Maxk+5][Maxk+5];
int pow_2[Maxk+5];
ll n;
int k;
int len;
void init(){
	pow_2[0]=1;
	C[0][0]=1;
	for(int i=1;i<=k;i++){
		C[i][0]=C[i][i]=1;
		for(int j=1;j<i;j++){
			C[i][j]=(C[i-1][j]+C[i-1][j-1])%Mod;
		}
	}
	for(int i=1;i<=k+1;i++){
		pow_2[i]=(pow_2[i-1]<<1)%Mod;
	}
}
struct Matrix{
	int a[2*Maxk+5][2*Maxk+5];
	void init(){
		for(int i=0;i<len;i++){
			a[i][i]=1;
		}
	}
	friend Matrix operator *(Matrix a,Matrix b){
		Matrix ans;
		for(int i=0;i<len;i++){
			for(int j=0;j<len;j++){
				ans.a[i][j]=0;
				for(int k=0;k<len;k++){
					ans.a[i][j]=(ans.a[i][j]+1ll*a.a[i][k]*b.a[k][j])%Mod;
				}
			}
		}
		return ans;
	}
}trans,ans;
Matrix quick_power(Matrix a,ll b){
	Matrix ans;
	ans.init();
	while(b){
		if(b&1){
			ans=ans*a;
		}
		a=a*a;
		b>>=1;
	}
	return ans;
}
int main(){
	scanf("%lld%d",&n,&k);
	init();
	len=(k<<1)+3;
	ans.a[0][0]=1;
	trans.a[0][0]=1;
	for(int i=0;i<=k;i++){
		ans.a[0][i+1]=1;
	}
	for(int i=k+2;i<len;i++){
		trans.a[i-(k+1)][i]=1;
	}
	for(int i=0;i<=k;i++){
		for(int j=0;j<=i;j++){
			trans.a[j+1][i+1]=C[i][j];
		}
		for(int j=0;j<=i;j++){
			trans.a[j+(k+2)][i+1]=1ll*C[i][j]*pow_2[i-j]%Mod;
		}
		trans.a[0][i+1]=(pow_2[i]+1)%Mod;
	}
	if(n==1){
		printf("%d
",ans.a[0][k+1]);
		return 0;
	}
	ans=ans*quick_power(trans,n-1);
	printf("%d
",ans.a[0][k+1]);
	return 0;
}

(O(k^2log n + 8klog n))

#include <cstdio>
typedef long long ll;
const int Maxk=40;
const int Mod=1000000007;
int C[Maxk+5][Maxk+5];
ll n;
int k;
void init(){
	C[0][0]=1;
	for(int i=1;i<=Maxk;i++){
		C[i][0]=C[i][i]=1;
		for(int j=1;j<i;j++){
			C[i][j]=(C[i-1][j]+C[i-1][j-1])%Mod;
		}
	}
}
struct Matrix{
	int a[2][2];
	void init(){
		a[0][0]=a[1][1]=1;
		a[0][1]=a[1][0]=0;
	}
	void clear(){
		a[0][0]=a[0][1]=a[1][0]=a[1][1]=0;
	}
	friend Matrix operator *(Matrix a,Matrix b){
		Matrix ans;
		for(int i=0;i<2;i++){
			for(int j=0;j<2;j++){
				ans.a[i][j]=0;
				for(int k=0;k<2;k++){
					ans.a[i][j]=(ans.a[i][j]+1ll*a.a[i][k]*b.a[k][j])%Mod;
				}
			}
		}
		return ans;
	}
	friend Matrix operator +(Matrix a,Matrix b){
		Matrix ans;
		for(int i=0;i<2;i++){
			for(int j=0;j<2;j++){
				ans.a[i][j]=(a.a[i][j]+b.a[i][j])%Mod;
			}
		}
		return ans;
	}
	friend Matrix operator *(Matrix a,int b){
		Matrix ans;
		for(int i=0;i<2;i++){
			for(int j=0;j<2;j++){
				ans.a[i][j]=1ll*a.a[i][j]*b%Mod;
			}
		}
		return ans;
	}
}fib,f[Maxk+5],t[Maxk+5],tmp;
int pow_m[Maxk+5];
void solve(ll n){
	if(n==1){
		for(int i=0;i<=k;i++){
			f[i]=fib;
		}
		tmp=fib;
		return;
	}
	ll m=(n>>1);
	solve(m);
	for(int i=0;i<=k;i++){
		t[i]=f[i];
		f[i].clear();
	}
	pow_m[0]=1;
	for(int i=1,tmp_m=m%Mod;i<=k;i++){
		pow_m[i]=1ll*pow_m[i-1]*tmp_m%Mod;
	}
	for(int i=0;i<=k;i++){
		for(int j=0;j<=i;j++){
			f[i]=f[i]+t[j]*(1ll*C[i][j]*pow_m[i-j]%Mod);
		}
		f[i]=f[i]*tmp;
	}
	for(int i=0;i<=k;i++){
		f[i]=f[i]+t[i];
	}
	tmp=tmp*tmp;
	if(n&1){
		tmp=tmp*fib;
		pow_m[0]=1;
		for(int i=1,tmp_m=n%Mod;i<=k;i++){
			pow_m[i]=1ll*pow_m[i-1]*tmp_m%Mod;
		}
		for(int i=0;i<=k;i++){
			f[i]=f[i]+tmp*pow_m[i];
		}
	}
}
int main(){
	fib.a[0][0]=fib.a[0][1]=fib.a[1][0]=1;
	fib.a[1][1]=0;
	scanf("%lld%d",&n,&k);
	init();
	solve(n);
	printf("%d
",f[k].a[0][0]);
	return 0;
}