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 imes 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 imes 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=egin{pmatrix}s_n n^0 n^1 ... n^k\end{pmatrix}$$经过一次转移，我们希望得到$$vec f_{n+1}=egin{pmatrix}s_{n+1} n^0 n^1 ... n^k\end{pmatrix}$$我们需要设计一个矩阵(A)使得$$Avec f_n=vec f_{n+1}$$由之前推出的结论，我们可以得到$$A=egin{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^0end{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;
}