bzoj5015: [Snoi2017]礼物

Description

热情好客的请森林中的朋友们吃饭，他的朋友被编号为 1～N，每个到来的朋友都会带给他一些礼物：。其中，第
一个朋友会带给他 1 个，之后，每一个朋友到来以后，都会带给他之前所有人带来的礼物个数再加他的编号的 K
次方那么多个。所以，假设 K=2，前几位朋友带来的礼物个数分别是：1,5,15,37,83假设 K=3，前几位朋友带来的
礼物个数分别是：1,9,37,111现在，好奇自己到底能收到第 N 个朋友多少礼物，因此拜托于你了。已知 N,K请输
出第 N 个朋友送的礼物个数 mod1000000007。

Input

第一行，两个整数 N,K
N≤10^18,K≤10

Output

一个整数，表示第 N 个朋友送的礼物个数 mod1000000007。

Sample Input

4 2

Sample Output

37

题解

记第\(n\)个朋友送的礼物数为\(a_n\)，前\(n\)个朋友送的礼物数之和为\(s_n\)，那么显然有\(s_n=2\times s_{n-1}+n^k\)，\(a_n=s_n-s_{n-1}\)
由于题目中\(N\)很大，直接递推显然会超时，考虑用矩阵快速幂加速。

注意到\(n^k\)这一项与\(n\)有关，因此我们在转移时要顺带考虑这一项如何转移。根据二项式定理$$(n+1)^k=\sum_{i=0}kC_k^i\times n^{k-i}$$我们只要知道了\(n^0\),\(n^1\),...,\(n^k\)，就可以算出\((n+1)^0\),\((n+1)^1\),...,\((n+1)^k\).那么我们可以让\(n^0\),\(n^1\),...,\(n^k\)随矩阵一同转移。
首先设计向量$$\vec f_n=\begin{pmatrix}s_n\ n^0\ n^1\ ...\ n^k\\end{pmatrix}$$经过一次转移，我们希望得到$$\vec f_{n+1}=\begin{pmatrix}s_{n+1}\ n^0\ n^1\ ...\ n^k\\end{pmatrix}$$我们需要设计一个矩阵\(A\)使得$$A\vec f_n=\vec f_{n+1}$$由之前推出的结论，我们可以得到$$A=\begin{bmatrix}2 & C_k^k & C_k^{k-1} & ... & C_k^0 \ 0 & C_0^0 & 0 & ... & 0\ 0 & C_1^1 & C_1^0 & ... & 0\ ... & ... & ... & ... & ...\ 0 & C_k^k & C_k^{k-1} & ... & C_k^0\end{bmatrix}$$

代码

#include<bits/stdc++.h>
#define MAXN 15
#define P 1000000007
#define LL long long
using namespace std;
struct Mat{
	LL a[MAXN][MAXN];
	int n;
	void Init(){for(int i=1;i<=n;i++)a[i][i]=1;}
	Mat(){memset(a,0,sizeof(a));}
}A,B,C;
struct Vec{
    LL a[MAXN];
    int n;
    Vec(){memset(a,0,sizeof(a));}
}f,g;
inline Mat operator * (Mat &x,Mat &y){
	Mat z;z.n=x.n;
	for(int i=1;i<=x.n;i++){
		for(int j=1;j<=x.n;j++){
			for(int k=1;k<=x.n;k++){
				z.a[i][j]=(z.a[i][j]+x.a[i][k]*y.a[k][j]%P)%P;
			}
		}
	}
	return z;
}
inline Vec operator * (Mat &x,Vec &y){
	Vec z;z.n=x.n;
	for(int i=1;i<=x.n;i++){
		for(int k=1;k<=x.n;k++){
			z.a[i]=(z.a[i]+x.a[i][k]*y.a[k]%P)%P;
		}
	}
	return z;
}
inline Mat qp(Mat x,LL y){
	Mat z;z.n=x.n;z.Init();
	while(y){
		if(y&1)z=x*z;
		x=x*x;y>>=1;
	}
	return z;
} 
LL N;
int K;
int main(){
	scanf("%lld%d",&N,&K);
	f.a[1]=1;f.n=K+2;
	for(int i=2;i<=K+2;i++)f.a[i]=1;
	A.a[2][2]=1;
	for(int i=3;i<=K+2;i++){
		A.a[i][2]=A.a[i][i]=1;
		for(int j=3;j<i;j++){
			A.a[i][j]=(A.a[i-1][j]+A.a[i-1][j-1])%P;
		}
	}
	A.a[1][1]=2;A.n=K+2;
	for(int i=2;i<=K+2;i++)A.a[1][i]=A.a[K+2][i];
	B=qp(A,N-2);
	f=B*f;g=A*f;
	printf("%lld",(g.a[1]-f.a[1]+P)%P);
	return 0; 
}