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}(n\ge 3),A_i(k)=F_i\times i^k\) ，求 \(\sum_{i=1}^{n} A_i(k)\) 。

\(n\leq 10^{17},k\leq 40\)

Bonus: \(k\leq 200\)

---

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

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

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

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

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

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

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

\[f(2n,k)=\sum_{i=0}^{2n}G^{i}\times i^k \\ =\sum_{i=0}^{n}G^i\times i^k + G^n\sum_{i=0}^{n}G^i\times (i+n)^k \\ =f(n,k)+G^n\sum_{i=0}^{n}G^i\sum_{j=0}^{k}\binom{k}{j}i^jn^{k-j} \\ =f(n,k)+G^n\sum_{j=0}^{k}\binom{k}{j}n^{k-j}\sum_{i=0}^{n}G^ii^j \\ =f(n,k)+G^n\sum_{j=0}^{k}\binom{k}{j}n^{k-j}f(n,j)\]

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

接下来是代码：

\(O(k^3\log 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\n",ans.a[0][k+1]);
		return 0;
	}
	ans=ans*quick_power(trans,n-1);
	printf("%d\n",ans.a[0][k+1]);
	return 0;
}

\(O(k^2\log n + 8k\log 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\n",f[k].a[0][0]);
	return 0;
}