Matrix Power Series(POJ 3233)

• 原题如下：
  Matrix Power Series

  Time Limit: 3000MS		Memory Limit: 131072K
  Total Submissions: 28044		Accepted: 11440

  Description

  Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

  Input

  The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

  Output

  Output the elements of S modulo m in the same way as A is given.

  Sample Input

  2 2 4
  0 1
  1 1

  Sample Output

  1 2
  2 3

• 题解：构造矩阵：
  
  

  此时，令S_k=I+A+…+A^k-1，则有：
  
  通过计算这个矩阵的k次幂，就可求出A的累乘和，时间复杂度为O(n³logk)

• 代码：

  #include <cstdio>
  #include <cctype>
  #define number s-'0'
  #include <cstring>
  #include <vector>
  
  using namespace std;
  
  typedef vector<int> vec;
  typedef    vector<vec> mat;
  
  int n,k,m;
  mat A;
  
  void read(int &x)
  {
      char s;
      x=0;
      bool flag=0;
      while (!isdigit(s=getchar()))
          (s=='-')&&(flag=true);
      for (x=number; isdigit(s=getchar());x=x*10+number);
      (flag)&&(x=-x);
  }
  
  void write(int x)
  {
      if (x<0)
      {
          putchar('-');
          x=-x;
      }
      if (x>9) write(x/10);
      putchar(x%10+'0');
  }
  
  mat mul(mat &A, mat &B)
  {
      mat C(A.size(), vec(B[0].size()));
      for (int i=0; i<A.size(); i++)
      {
          for (int j=0; j<B[0].size(); j++)
          {
              for (int k=0; k<B.size(); k++)
              {
                  C[i][j]=(C[i][j]+A[i][k]*B[k][j])%m;
              }
          }
      }
      return C;
  }
  
  mat pow(mat A, int n)
  {
      mat B(A.size(), vec(A.size()));
      for (int i=0; i<A.size(); i++) B[i][i]=1;
      while (n>0)
      {
          if (n&1) B=mul(B, A);
          A=mul(A, A);
          n>>=1;
      }
      return B;
  }
  
  int main(int argc, char * argv[])
  {
      read(n);read(k);read(m);
      A=mat (n,vec(n));
      mat B(n*2, vec(n*2));
      for (int i=0; i<n; i++)
      {
          for (int j=0; j<n; j++)
          {
              read(A[i][j]);
              B[i][j]=A[i][j];
          }
          B[n+i][i]=B[n+i][n+i]=1;
      }
      B=pow(B,k+1);
      for (int i=0; i<n; i++)
      {
          for (int j=0; j<n; j++)
          {
              int a=B[n+i][j]%m;
              if (i==j) a=(a+m-1)%m;
              printf("%d%c", a, j+1==n?'
  ':' ');
          }
      }
  }