题目大意：求长度为1~n的珠子串中恰有K种珠子的方案数（mod 1234567891）。
做法：用f[i][j]表示长度为i的珠子中恰有j种珠子的方案数，状态转移方程为：f[i][j]=f[i-1][j-1]*(k-j+1)+f[i-1][j]*j，结果是f[1][k]+...+f[n][k]。可是n可以达到10^9，传统的二维DP的时间复杂度为O(nk)，无法承受。于是我们考虑用矩阵优化。设(k+1)*(k+1)的矩阵F[i]和A，其中：
F[i]为：
f[i][0] f[i][1] ... f[i][k]
0 0 ... 0
.
.
.
0 0 ... 0
A为：
0 k 0 ... 0
0 1 k-1 ... 0
0 0 2 ... 0
.
.
.
0 0 0 ... k
因此状态转移方程为：F[i]=F[i-1]*A，因此结果为（F[1]+F[2]+...+F[n]）的第1行第k+1列的数。
进一步推算，F[1]+F[2]+...+F[n]=F[1]+F[1]*A+...+F[1]*A^(n-1)=F[1]*(E+A+...+A^(n-1))，因此只要二分求出A+...+A^(n-1)的结果后，代入公式计算即可。（F[1]的第1行第2列为k，其余全为0）
以下为本人代码：
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define mod 1234567891
using namespace std;
int T,n,k;
struct matrix
{
__int64 s[31][31];
}M[31],ans,E,F; //M[i]为A^(2^i)，E为单位矩阵，F为F[1]
void clear(matrix &A) //将A指向的矩阵清零
{
for(int i=0;i<=k;i++)
for(int j=0;j<=k;j++)
A.s[i][j]=0;
}
matrix pl(matrix A,matrix B) //矩阵加法
{
matrix S;
for(int i=0;i<=k;i++)
for(int j=0;j<=k;j++)
S.s[i][j]=(A.s[i][j]+B.s[i][j])%mod;
return S;
}
matrix mult(matrix A,matrix B) //矩阵乘法
{
matrix S;
clear(S);
for(int i=0;i<=k;i++)
for(int j=0;j<=k;j++)
for(int x=0;x<=k;x++)
S.s[i][j]=(S.s[i][j]+A.s[i][x]*B.s[x][j])%mod;
return S;
}
matrix power(int x) //矩阵快速幂求A^x
{
matrix S=E;
int i=0;
while(x!=0)
{
if (x&1) S=mult(S,M[i]);
i++;x>>=1;
}
return S;
}
matrix pow_sum(int x) //求(A+A^2+...+A^x)
{
if (x==0) {matrix S;clear(S);return S;}
if (x==1) return M[0];
if (x%2!=0) return pl(pow_sum(x-1),power(x));
else return mult(pow_sum(x/2),pl(power(x/2),E));
}
int main()
{
scanf("%d",&T);
while(T--)
{
scanf("%d%d",&n,&k);
clear(M[0]);clear(E);
for(int i=0;i<=k;i++)
{
M[0].s[i][i]=i;
if (i!=k) M[0].s[i][i+1]=k-i;
E.s[i][i]=1;
}
for(int i=1;i<=30;i++) M[i]=mult(M[i-1],M[i-1]);
clear(F);F.s[0][1]=k;
ans=mult(F,pl(pow_sum(n-1),E)); //计算答案
printf("%I64d
",ans.s[0][k]%mod);
}
return 0;
}