【HDU2294】Pendant-DP矩阵优化

题目大意：求长度为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\n",ans.s[0][k]%mod);
  }
  
  return 0;
}