hdu 2294 Pendant DP+矩阵优化

题意：有K种珍珠，每种N颗，求这些珍珠组成的长度在1~N之间的项链中包含K种珍珠的有多少种（答案模1234567891）？

其中，1<=K<=30,1<=N<=1000,000,000

共有十组case。

首先，这题是DP题，得出状态转移方程 dp[i][j] = j*dp[i-1][j]+(k-j+1)*dp[i-1][j-1].其中 dp[i][j]表示包含j种珍珠且长度为i的项链的种数。若用滚动数组求解，则答案则为dp[1][k]+dp[2][k]+...+dp[N][k],但由于N太大，这样求明显会超时。于是对于数组dp[i]转移到dp[i+1]过程，可以构造矩阵完成。这样对于后面求和，便可以用二分+矩阵快速幂取模完成。注意，二分求和的过程要写成非递归形式，否则栈溢出，至少我的程序是这样的。

源代码：

#include <stdio.h>
#include <string.h>
#define MOD 1234567891
typedef struct
{
  __int64 matrix[31][31];
}Matrix;

int K;

Matrix ProE()
{
  Matrix res;
  int i,j;

  for(i = 1;i<=K;i++)
   for(j = 1;j<=K;j++)
   {
    if(i==j)
     res.matrix[i][j] = 1;
    else
     res.matrix[i][j] = 0;
   }
  return res;
}

Matrix multi(Matrix x,Matrix y)
{
  int i,j,k;
  Matrix res;

  memset(res.matrix,0,sizeof(res.matrix));
  for(i = 1;i<=K;i++)
   for(j = 1;j<=K;j++)
   {
    if(x.matrix[i][j])
      for(k = 1;k<=K;k++)
     res.matrix[i][k]=(res.matrix[i][k]+x.matrix[i][j]*y.matrix[j][k])%MOD;
   }
 return res;
}

Matrix powermod(Matrix x,int n)
{
 Matrix res;

 res = ProE();
 for(;n;n>>=1)
 {
  if(n&1)
    res = multi(res,x);
  x = multi(x,x);
 }
 return res;
}

Matrix add(Matrix x,Matrix y)
{
 int i,j;
 Matrix res;

 for(i = 1;i<=K;i++)
  for(j = 1;j<=K;j++)
   res.matrix[i][j] = (x.matrix[i][j]+y.matrix[i][j])%MOD;
 return res;
}

Matrix powersum(Matrix x,int n)
{
  Matrix keep,res;
  memset(res.matrix,0,sizeof(res.matrix));
  keep = ProE();
  for(;n;n>>=1)
  {
   if(n&1)
   res = add(res,multi(powermod(x,n),keep));
   keep = multi(keep,add(powermod(x,n/2),ProE()));
  }
  return res;
}

int main()
{
   int T;
   __int64 N,i;
   Matrix res;

   scanf("%d",&T);
   while(T--)
   {
    scanf("%I64d%I64d",&N,&K);
    if(N==1)
    {
     printf("%I64d\n",(N*K)%MOD);
     continue;
    }
    memset(res.matrix,0,sizeof(res.matrix));
       for(i = 1;i<K;i++)
       {
         res.matrix[i][i] = K+1-i;
         res.matrix[i][i+1] = K-(K+1-i)+1;
       }
    res.matrix[K][K] = 1;
     res = add(ProE(),powersum(res,N-1));
    printf("%I64d\n",(res.matrix[1][K]*K)%MOD);
   }
   return 0;
}