BZOJ 5015: [Snoi2017]礼物

设 \(f_i\) 为答案，\(s_i\) 为前缀和
那么 \(s_i = 2s_{i-1}+i^k\)
这个根据二项式系数可以矩阵快速幂

#include <bits/stdc++.h>
#define pb push_back
#define fi first
#define se second
#define pii pair<int, int>
#define lp p << 1
#define rp p << 1 | 1
#define mid ((l + r) >> 1)
#define ll long long
#define db double
#define rep(i,a,b) for(int i=a;i<b;i++)
#define per(i,a,b) for(int i=b-1;i>=a;i--)
#define Edg int ccnt=1,head[N],to[N*2],ne[N*2];void addd(int u,int v){to[++ccnt]=v;ne[ccnt]=head[u];head[u]=ccnt;}void add(int u,int v){addd(u,v);addd(v,u);}
#define Edgc int ccnt=1,head[N],to[N*2],ne[N*2],c[N*2];void addd(int u,int v,int w){to[++ccnt]=v;ne[ccnt]=head[u];c[ccnt]=w;head[u]=ccnt;}void add(int u,int v,int w){addd(u,v,w);addd(v,u,w);}
#define es(u,i,v) for(int i=head[u],v=to[i];i;i=ne[i],v=to[i])
const int MOD = 1e9 + 7;
void M(int &x) {if (x >= MOD)x -= MOD; if (x < 0)x += MOD;}
int qp(int a, int b = MOD - 2) {int ans = 1; for (; b; a = 1LL * a * a % MOD, b >>= 1)if (b & 1)ans = 1LL * ans * a % MOD; return ans % MOD;}
int gcd(int a, int b) { while (b) { a %= b; std::swap(a, b); } return a; }

const int N = 15;
int n;
ll m;

struct Mat {
	int mat[N][N];
	Mat() { memset(mat, 0, sizeof(mat)); }
	Mat operator * (const Mat &p) const {
		Mat c;
		rep (i, 1, n +3) {
			rep (j, 1, n + 3) {
				ll temp = 0;
				rep (k, 1, n + 3) {
					temp += 1LL * mat[i][k] * p.mat[k][j];
					if (temp >= (1LL << 60)) temp %= MOD;
				}
				c.mat[i][j] = temp % MOD;
			}
		}
		return c;
	}
} base;

Mat qp(Mat a, ll b) {
	Mat c;
	rep (i, 1, n + 1) c.mat[i][i] = 1;
	while (b) {
		if (b & 1) c = a * c;
		a = a * a;
		b >>= 1;
	}
	return c;
}

int C[N][N];

int main() {
	C[0][0] = 1;
	rep (i, 1, N) {
		C[i][0] = C[i][i] = 1;
		rep (j, 1, i) M(C[i][j] = C[i - 1][j] + C[i - 1][j - 1]);
	}
	scanf("%lld%d", &m, &n);
	Mat base;
	base.mat[1][1] = 2;
	rep (j, 2, n + 3) base.mat[j][1] = C[n][j - 2];
	rep (j, 2, n + 3) {
		int cnt = 0;
		rep (i, j, n + 3) {
			base.mat[i][j]= C[n - j + 2][cnt++];
		}
	}
	Mat res = qp(base, m - 2);
	int ans = 0;
	rep (i, 1, n + 3)
		M(ans += res.mat[i][1]);
	res = base * res;
	int ans2 = 0;
	rep (i, 1, n + 3) 
		M(ans2 += res.mat[i][1]);
	M(ans2 = ans2 - ans);
	printf("%d\n", ans2);
	return 0;
}