PKU 3233 Matrix Power Series

总的思想就是二分，求矩阵的高次幂需要用到二分的思想，然后，求矩阵幂的和，再用一次二分。

 

S = A^1 + A^2 + ... + A^k

 

如果k是奇数，则

S = (A^1 + A^2 + ... + A^(k / 2)) + A^((k + 1) / 2) + A^((k + 1) / 2) * ((A^1 + A^2 + ... + A^(k / 2)) )

   = A^((k + 1) / 2) + ( A^((k + 1) / 2) + I)(A^1 + A^2 + ... + A^(k / 2))

 

如果k是偶数，则

S = (A^1 + A^2 + ... + A^(k / 2)) + A^(k / 2) * ((A^1 + A^2 + ... + A^(k / 2)) )

   = ( A^(k/ 2) + I)(A^1 + A^2 + ... + A^(k / 2))

 

#include <iostream>

 

using namespace std;

 

const int N = 31;

int n, k, m;

int A[N][N];

 

//用于计算矩阵的幂

int pResult[N][N];

int pTransform[N][N];

int pTemp[N][N];

void MatPow(const int&P)

{

      int i, j, k;

      int p = P;

      memset(pResult, 0, sizeof(pResult));

      for(i = 1; i <= n; ++i)

           pResult[i][i] = !pResult[i][i];

      memcpy(pTransform, A, sizeof(pTransform));

      while(p > 0)

      {

           if(p & 1)

           {

                 memset(pTemp, 0, sizeof(pTemp));

                 for(i = 1; i <= n; ++i)

                 {

                      for(j = 1; j <= n; ++j)

                      {

                            for(k = 1; k <= n; ++k)

                            {

                                  pTemp[i][j] += pTransform[i][k] * pResult[k][j];

                                  if(pTemp[i][j] >= m)

                                  {

                                       pTemp[i][j] = pTemp[i][j] %m;

                                  }

                            }

                      }

                 }

                 memcpy(pResult, pTemp, sizeof(pResult));

           }

 

           p = p >> 1;

           memset(pTemp, 0, sizeof(pTemp));

           for(i = 1; i <= n; ++i)

           {

                 for(j = 1; j <= n; ++j)

                 {

                      for(k = 1; k <= n; ++k)

                      {

                            pTemp[i][j] += pTransform[i][k] * pTransform[k][j];

                            if(pTemp[i][j] >= m)

                            {

                                  pTemp[i][j] = pTemp[i][j] %m;

                            }

                      }

                 }

           }

           memcpy(pTransform, pTemp, sizeof(pTransform));

      }

      return;

}

 

//用于计算矩阵幂的和

int sResult[N][N];

int sTransform[N][N];

int sTemp[N][N];

void MatPowSum(const int& P)

{

      int i, j, k;

      int p = P;

      memset(sResult, 0, sizeof(sResult));

      memset(sTransform, 0, sizeof(sTransform));

      for(i = 1; i <= n; ++i)

      {

           sTransform[i][i] = !sTransform[i][i];

      }

 

      while(p > 0)

      {

           if(p & 1)

           {

                 MatPow((p + 1) / 2);

                 memset(sTemp, 0, sizeof(sTemp));

                 for(i = 1; i <= n; ++i)

                 {

                      for(j = 1; j <= n; ++j)

                      {

                            for(k = 1; k <= n; ++k)

                            {

                                  sTemp[i][j] += sTransform[i][k] * pResult[k][j];

                                  if(sTemp[i][j] >= m)

                                  {

                                       sTemp[i][j] = sTemp[i][j] %m;

                                  }

                            }

                      }

                 }

 

                 for(i = 1; i <= n; ++i)

                 {

                      for(j = 1; j <= n; ++j)

                      {

                            sResult[i][j] += sTemp[i][j];

                            if(sResult[i][j] >= m)

                            {

                                  sResult[i][j] %= m;

                            }

                      }

                 }

                 //A的(k+1)/2次方与单位阵相加

                 for(i = 1; i <= n; ++i)

                 {

                      pResult[i][i] += 1;

                 }

                 //原来的转移矩阵* ( A的(k+1)/2次方+ I )

                 memset(sTemp, 0, sizeof(sTemp));

                 for(i = 1; i <= n; ++i)

                 {

                      for(j = 1; j <= n; ++j)

                      {

                            for(k = 1; k <= n; ++k)

                            {

                                  sTemp[i][j] += sTransform[i][k] * pResult[k][j];

                                  if(sTemp[i][j] >= m)

                                  {

                                       sTemp[i][j] = sTemp[i][j] %m;

                                  }

                            }

                       }

                 }

                 memcpy(sTransform, sTemp, sizeof(sTransform));

           }

           else

           {

                 MatPow(p / 2);

                 //A的k/2次方与单位阵相加

                 for(i = 1; i <= n; ++i)

                 {

                      pResult[i][i] += 1;

                 }

                 //原来的转移矩阵* ( A的(k+1)/2次方+ I )

                 memset(sTemp, 0, sizeof(sTemp));

                 for(i = 1; i <= n; ++i)

                 {

                      for(j = 1; j <= n; ++j)

                      {

                            for(k = 1; k <= n; ++k)

                            {

                                  sTemp[i][j] += sTransform[i][k] * pResult[k][j];

                                  if(sTemp[i][j] >= m)

                                  {

                                       sTemp[i][j] = sTemp[i][j] %m;

                                  }

                            }

                      }

                 }

                 memcpy(sTransform, sTemp, sizeof(sTransform));

           }

 

           p = p >> 1;

          

      }

      return;

}

 

int main()

{

      int i, j;

      scanf("%d %d %d", &n, &k, &m);

      for(i = 1; i <= n; ++i)

      {

           for(j = 1; j <= n; ++j)

           {

                 scanf("%d", &A[i][j]);

           }

      }

     

      MatPowSum(k);

 

      //output

      for(i = 1; i <= n; ++i)

      {

           for(j = 1; j <= n; ++j)

           {

                 printf("%d ", sResult[i][j]);

           }

           printf("/n");

      }

      return 0;

}