这道题是一道十分标准的矩阵快速幂题目。

那么我们先来看一下两个矩阵是怎么相乘的。

首先，两个矩阵 $A$ 与 $B$ 相乘的必要条件为 $A$ 的列等于 $B$ 的行数，相乘的结果矩阵 $C$ 的行数等于 $A$ 的行数，$C$ 的列数等于 $B$ 的列数。

在做这道题时我们可以先构造一个关于矩阵的结构体，代码如下：

struct mat
{
    long long m,n;
    long long a[30][30]; 
    mat operator * (mat& b)  //重构乘法运算符
    {
        mat c;
        memset(c.a,0,sizeof(c.a));
        c.m=m;
        c.n=b.n;
        for(int i=1;i<=m;i++)
            for(int j=1;j<=b.n;j++)
                for(int p=1;p<=n;p++)
                {
                    c.a[i][j]+=(a[i][p]*b.a[p][j])%1000000007;
                    c.a[i][j]%=1000000007;
        }
        return c;
    }
};

但是我们会有一个问题：矩阵里面应该放什么变量呢？

首先我们会想到放入第 $i$ 个人的礼物数量，

其次我们会用到之前所有人的礼物数量之和，因为计算第 $i$ 个人礼物数量会用到。

我们在计算第 $i$ 个人的礼物数量需要用到 $i^{k}$,所以在矩阵中还要有 $i^{k}$ 。

同时需要计算 $i^{k}$ 还需要 $i^{k-1}$ ，而 $i^{k-1}$ 还需要 $i^{k-2}$。

所以我们找到了矩阵中需要放入的变量：当前朋友的礼物、$i^{k}$、$i^{k-1} \cdots i$、$1$ 。

我构造的矩阵是这样的：

$\begin{vmatrix} num_{i} &i^{k} &i^{k-1} &\cdots &i &1 &S_{i} \end{vmatrix}$

那么接下来的一步就变成矩阵快速幂了，

那么这种矩阵应该乘一个怎样的矩阵才能变化到下一步呢？

这个矩阵的下一步是这样的：

$\begin{vmatrix} num_{i+1} &\left(i+1\right)^{k} &\left(i+1\right)^{k-1} &\cdots &\left(i+1\right) &1 &S_{i+1} \end{vmatrix}$

我们来一步一步地看：

$num_{i+1}=S_{i}+\left(i+1\right)^{k}$,

$S_{i+1}=S_{i}+num_{i}$

那么 $\left(i+1\right)^{k}$、$\left(i+1\right)^{k-1}$ 等又应该如何构造呢？

$\because$ 二项式定理，所以这个地方应该是一个杨辉三角，我们可以写一个函数来求。

代码如下：

long long yh[20][20];
long long yhsj(long long x,long long y)
{
    if(!yh[x][y])
    {
        if(y==1||y==x)
            yh[x][y]=1;
        else
            yh[x][y]=yh[x-1][y-1]+yh[x-1][y];
    }
    return yh[x][y];
}
void nb(mat& a,long long k)
{
    a.m=a.n=k+3;
    memset(a.a,0,sizeof(a.a));
    a.a[k+3][1]=1;
    a.a[k+3][k+3]=2;
    for(int i=1;i<=k+1;i++)
        for(int j=1;j<=i;j++)
             a.a[k+3-j][k+3-i]=yhsj(i,j);
    for(int i=2;i<=k+2;i++)
        a.a[i][1]=a.a[i][k+3]=a.a[i][2];
}

其中 $\operatorname{yhsj}$ 函数用于计算杨辉三角，数组 $yh$ 用于记忆化搜索杨辉三角的某一项。

接下来就是快速幂的部分了。

实际上这个部分很简单，这要把整数的快速幂稍微修修改可以了。

递归形式：

mat mut(mat a,long long t)
{
    if(t==1)
        return a;
    if(t==0)
    {
        mat b;
        b.m=b.n=k+3;
        memset(b.a,0,sizeof(b.a));
        for(int i=1;i<=k+3;i++)
            b.a[i][i]=1;
        return b;
    }
    mat ans;
    ans.m=ans.n=3;
    memset(ans.a,0,sizeof(ans.a));
    ans=mut(a,t/2);
    ans=ans*ans;
    if(t%2==1)
        ans=ans*a;
    return ans;
}

$while$ 循环形式：

mat mut(mat a,long long t)
{
    mat ans;
    memset(ans.a,0,sizeof(ans.a));
    ans.m=ans.n=k+3;
    for(int i=1;i<=k+3;i++)
    	ans.a[i][i]=1;
    while(t)
    {
        if(t&1)
        	ans=ans*a;
        a=a*a;
        t=t/2;
    }
    return ans;
}

下面给出最终代码：

#include<cstring>
#include<cstdio>
using namespace std;
long long yh[20][20];
long long n,k;
struct mat
{
	long long m,n;
	long long a[101][101];
	mat operator * (mat& b)
	{
		mat c;
		memset(c.a,0,sizeof(c.a));
		c.m=m;
		c.n=b.n;
		for(int i=1;i<=m;i++)
			for(int j=1;j<=n;j++)
				for(int p=1;p<=n;p++)
				{
					c.a[i][j]+=(a[i][p]*b.a[p][j])%1000000007;
					c.a[i][j]%=1000000007;
				}
		return c;
	}
	mat operator + (mat& b)
	{
		mat c;
		c.n=b.n;
		c.m=b.m;
		for(int i=1;i<=m;i++)
			for(int j=1;j<=n;j++)
				c.a[i][j]=a[i][j]+b.a[i][j];
		return c;
	}
};
long long yhsj(long long x,long long y)
{
	if(!yh[x][y])
	{
		if(y==1||y==x)
			yh[x][y]=1;
		else
			yh[x][y]=yh[x-1][y-1]+yh[x-1][y];
	}
	return yh[x][y];
}
void nb(mat& a,long long k)
{
	a.m=a.n=k+3;
	memset(a.a,0,sizeof(a.a));
	a.a[k+3][1]=1;
	a.a[k+3][k+3]=2;
	for(int i=1;i<=k+1;i++)
		for(int j=1;j<=i;j++)
			a.a[k+3-j][k+3-i]=yhsj(i,j);
	for(int i=2;i<=k+2;i++)
		a.a[i][1]=a.a[i][k+3]=a.a[i][2];
}
void nbnb(mat& a,long long k)
{
	memset(a.a,0,sizeof(a.a));
	a.m=1;
	a.n=k+3;
	for(int i=1;i<=k+3;i++)
		a.a[1][i]=1;
}
mat mut(mat a,long long t)
{
	if(t==1)
		return a;
	if(t==0)
	{
		mat b;
		b.m=b.n=k+3;
		memset(b.a,0,sizeof(b.a));
		for(int i=1;i<=k+3;i++)
			b.a[i][i]=1;
		return b;
	}
	mat ans;
	ans.m=ans.n=3;
	memset(ans.a,0,sizeof(ans.a));
	ans=mut(a,t/2);
	ans=ans*ans;
	if(t%2==1)
		ans=ans*a;
	return ans;
}
int main()
{
	scanf("%lld%lld",&n,&k);
	mat ans;
	nb(ans,k);
	mat op;
	nbnb(op,k);
	mat num=mut(ans,n-1);
	num=op*num;
	printf("%lld",num.a[1][1]);
	return 0;
}

不愧是我，最终代码正好打了100行。

本人第一次写题解，欢迎各位大佬指出我的错误。

建议先学二项式系数。

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define TP template<typename T>
#define TP_ template<typename T,typename ... T_>
TP void read(T &x)
{
	x=0;int f=1;char ch=getchar();
	for(;ch<'0'||ch>'9';ch=getchar())if(ch=='-')f=-1;
	for(;ch>='0'&&ch<='9';ch=getchar())x=(x<<1)+(x<<3)+(ch^48);
	x*=f;
}
TP_ void read(T &x,T_&...y){read(x);read(y...);}
TP void write(T x){if(x<0){putchar('-'),x=-x;}if(x>9)write(x/10);putchar(48+x%10);}
TP void writeln(T x){if(x<0){putchar('-'),x=-x;}if(x>9)write(x/10);putchar(48+x%10);puts("");}
TP_ void writeln(const T x,T_ ...y){write(x);putchar(' ');writeln(y...);}
typedef long long LL;
constexpr int P=1e9+7;
LL c[15][15],n,k;
struct matrix 
{
	LL a[21][21];
	matrix()
	{
		memset(a,0,sizeof(a));
	}
	friend matrix operator *(const matrix &A,const matrix &B)
	{
		matrix C;
		for(int i=1;i<=k+3;i++)
			for(int j=1;j<=k+3;j++)
				for(int l=1;l<=k+3;l++)
					C.a[i][j]=(C.a[i][j]+A.a[i][l]*B.a[l][j]%P)%P;
		return C;
	}
}F,f;
matrix qpow(matrix a,LL b)
{
	matrix ans;for(int i=1;i<=k+3;i++)ans.a[i][i]=1;
	for(;b;b>>=1)
	{
		if(b&1)ans=ans*a;
		a=a*a;
	}
	return ans;
}
void init()
{
	c[0][0]=1;
	for(int i=1;i<=k+1;i++)
		c[i][0]=1;
	for(int i=1;i<=k+1;i++)
		for(int j=1;j<=i;j++)
			c[i][j]=(c[i-1][j]+c[i-1][j-1])%P;
}
int main()
{
	read(n,k);
	init();LL x=1;
	for(int i=1;i<=k+1;i++)//f.a[1][1]~f.a[1][k+1] 存 2^0,2^1……2^k,f.a[1][k+2]存答案,f.a[1][k+3]存前面所有人的总数。 
		F.a[1][i]=x,x*=2;
	F.a[1][k+2]=1;
	F.a[1][k+3]=0;
	for(int i=1;i<=k+1;i++)
		for(int j=1;j<=i;j++)
			f.a[j][i]=c[i-1][j-1];
	
	f.a[k+3][k+2]=1;
	f.a[k+1][k+2]=1;
	f.a[k+2][k+2]=1;
	
	f.a[k+2][k+3]=1;
	f.a[k+3][k+3]=1;
	writeln((F*qpow(f,n-1)).a[1][k+2]);
	return 0;
}