【矩阵快速幂】Matrix Power Series

2019年3月30日 0条评论 669次阅读 0人点赞于衡

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

#include<iostream>
#include<sstream>
#include<map>
#include<cmath>
#include<fstream>
#include<vector>
#include<sstream>
#include<cstring>
#include<cstdio>
#include<stack>
#include<queue>
#include<bitset>
#include<ctime>
#include<string>
#include<iomanip>
#include<algorithm>
using namespace std;
#define INT __int64
const int INF = 0x3f3f3f3f;
const double esp = 0.00000001;
const double PI = acos(-1.0);
const int MY = 8;
const int MX = 30 + 5;
INT n, k, mod;
struct M
{
  INT p[MX][MX];
  void zero()
  {
    memset(p, 0, sizeof(p));
  }
  void unit()
  {
    for (int i = 0; i < n; ++i)
      p[i][i] = 1;
  }
}ans;
M operator + (const M& a, const M& b)
{
  M c;
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
      c.p[i][j] = (a.p[i][j] + b.p[i][j]) % mod;
  return c;
}
M operator * (const M& a, const M& b)
{
  M c;
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
    {
      c.p[i][j] = 0;
      for (int q = 0; q < n; ++q)
        c.p[i][j] += (a.p[i][q] * b.p[q][j]) % mod;
    }
  return c;
}
M operator ^ (M a, INT k)
{
  M b;
  b.zero();
  b.unit();
  while (k)
  {
    if (k & 1)
      b = b * a;
    a = a * a;
    k >>= 1;
  }
  return b;
}
M sum(INT k)
{
  if (k == 1)
    return ans;
  else
  {
    M t = sum(k / 2);
    if (k & 1)
    {
      M b = ans ^ ((k >> 1) + 1);
      return t + b + t * b;
    }
    else
    {
      M b = ans ^ (k >> 1);
      return  t + t * b;
    }
  }
}
int main()
{
  while (~scanf("%I64d%I64d%I64d", &n, &k, &mod))
  {
    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
      {
        scanf("%I64d", &ans.p[i][j]);
        ans.p[i][j] %= mod;
      }
    M c = sum(k);
    for (int i = 0; i < n; ++i)
    {
      cout << c.p[i][0] % mod;
      for (int j = 1; j < n; ++j)
        cout << " " << c.p[i][j] % mod;
      cout << endl;
    }
  }
  return 0;
}