k次幂（构造矩阵的技巧）

第2题     k次幂 查看测评数据信息

输入格式

 

一行，两个整数，n和k。  1<=n<=1000000000,   1<=k<=50。

答案模1000000007。

 

输出格式

 

一个整数。

 

输入/输出例子1

输入：

5 1

 

输出：

15

 

样例解释

无

 

 

先考虑dp

f(i): 1^k+2^k+...+i^k
f(i)=f(i-1)+i^k

 

矩阵加速：
[f(i-1), i^k] * [ ] = [f(i), (i+1)^k]
(i+1)^k比较难求，我们展开（二项式定理）

(i+1)^k=C(k, 1)*i + C(k, 2)*i^2 + C(k, 3)*i^3 .....  +C(k, k-1)*i^(k-1) + C(k, k)*i^k
组合数可以预处理，但是i^(k-1)， i^(k-2)....i^0怎么算？

好像没法算。所以必须放在矩阵中维护。

[f(i-1), i^k， i^(k-1)，i^(k-2)，.....，i^0]
* [ ] =
[f(i), (i+1)^k， (i+1)^(k-1)，(i+1)^(k-2)，.....，(i+1)^0]

那么(i+1)^k就够算了嘛！

我们再看看(i+1)^（k-1）如何算

展开，发现条件全部已知了。

(i+1)^（k-1）=  C(k-1, 1)*i + C(k-1, 2)*i^2 + ...... + C(k-1, k-2)*i^(k-2) + C(k-1, k-1)*i^(k-1)

 

这个时候可以构造未知矩阵了：

第一列特殊一点，根据“f(i)=f(i-1)+i^k”来填

第二列往后开始有规律了

(i+1)^k，首先肯定不需要f(i-1)，然后按照二项式定理展开，

(i+1)^k=C(k, 1)*i + C(k, 2)*i^2 + C(k, 3)*i^3 .....  +C(k, k-1)*i^(k-1) + C(k, k)*i^k

根据上面的式子又可以继续填了。

然后到第三列，第四列同理。先展开，再填。

[1, 0,             0                          .. ]
[1, C(k, k),    0                           .. ]
[0, C(k, k-1), C(k-1, k-1)            .. ]
[0, C(k, k-2), C(k-1, k-1-1)           ]
[. .. .. .. ]
[. .. .. .. ]
[. .. .. .. ]
[0 C(k, k-k),   C(k-1, k-1-(k-1)) .. ]

 

当i=1时，
[0, 1, 1, 1......1] * [] = [结果]

 

 

注意一下，

本代码和构造矩阵有点不一样，是以

[1, 0,             0                           .. ]
[1, C(k, 0),    0                           .. ]
[0, C(k, 1), C(k-1, 0)                  .. ]
[0, C(k, 2), C(k-1, 1)                     ]
[. .. .. .. ]
[. .. .. .. ]
[. .. .. .. ]
[0 C(k, k),   C(k-1, k-1)              .. ]

构造的，但其实本质一样，因为C(n, r)=C(n, n-r)

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N=105, Mod=1000000007;

struct mp
{
	int n, m;
	int a[N][N];
	
	void init(int row, int col, bool isI)
	{
		n=row, m=col;
		memset(a, 0, sizeof a);
		
		if (isI)
			for (int i=1; i<=row; i++) a[i][i]=1;
	}	
}A, B;
mp operator * (const mp A, const mp B)
{
	mp C;
	C.init(A.n, B.m, 0);
	for (int i=1; i<=A.n; i++)
		for (int j=1; j<=B.m; j++)
			for (int k=1; k<=A.m; k++)
				C.a[i][j]=(C.a[i][j]+A.a[i][k]*B.a[k][j])%Mod;
		
	return C;
}
mp qpow_mp(mp A, int k)
{
	mp res;
	res.init(A.n, A.m, 1);
	while (k>0)
	{
		if (k&1) res=res*A;
		A=A*A;
		k>>=1;
	}
	
	return res;
}
int n, k, yh[N][N];
signed main()
{
	scanf("%lld%lld", &n, &k);
	A.init(1, 2+k, 0);
	A.a[1][1]=0;
	for (int i=2; i<=2+k; i++) A.a[1][i]=1;
	
	for (int i=1; i<=100; i++)
	{
		yh[i][i]=yh[i][1]=1;
		for (int j=2; j<i; j++) yh[i][j]=(yh[i-1][j]+yh[i-1][j-1])%Mod;
	}
	
	B.init(2+k, 2+k, 0);
	B.a[1][1]=B.a[2][1]=1;
	for (int i=2; i<=2+k; i++)
		for (int j=i; j<=2+k; j++)
			B.a[j][i]=yh[k-i+2+1][j-i+1];
			
	A=A*qpow_mp(B, n);
		
	printf("%lld", A.a[1][1]);
	return 0;
}