ぱーぽーの競プロ記

競技プログラミングに関することを書きます。

AOJ 1327 : One-Dimensional Cellular Automaton

概要


問題文はこちら
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1327

N個のセルからなるセルオートマトンS(i, t)があり、以下の更新則を持つ。

S(i, t + 1) = (A × S(i − 1, t) + B × S(i, t) + C × S(i + 1, t)) mod M

初期状態S(i, 0)が与えられたとき、S(i, t)の値を求めよ。(0 <= i < N)

解法


上記の更新則は以下のような行列で表すことができる。

(例) N=5のとき

{ 
 \begin{pmatrix} S[0][t + 1] \\ S[1][t + 1] \\ S[2][t + 1] \\ S[3][t + 1] \\ S[4][t + 1] \end{pmatrix} = \begin{pmatrix} b & c & 0 & 0 & 0 \\ a & b & c & 0 & 0 \\ 0 & a & b & c & 0 \\ 0 & 0 & a & b & c \\ 0 & 0 & 0 & a & b \end{pmatrix} \begin{pmatrix} S[0][t ] \\ S[1][t] \\ S[2][t] \\ S[3][t] \\ S[4][t] \end{pmatrix}
}


これをもう少し変形する。

{ 
 A = \begin{pmatrix} b & c & 0 & 0 & 0 \\ a & b & c & 0 & 0 \\ 0 & a & b & c & 0 \\ 0 & 0 & a & b & c \\ 0 & 0 & 0 & a & b \end{pmatrix} 
}

とすると、

{ 
 \begin{pmatrix} S[0][t] \\ S[1][t] \\ S[2][t] \\ S[3][t] \\ S[4][t] \end{pmatrix} = A^{t} \begin{pmatrix} S[0][0] \\ S[1][0] \\ S[2][0] \\ S[3][0] \\ S[4][0] \end{pmatrix}
}

となる。

t<=10^9なのでこれを高速に計算できれば解ける。

繰り返し二乗法を用いることでO(log t)で求めることができる。(蟻本参照)

ソースコード

#include <bits/stdc++.h>

#define REP(i, x, n) for(int i = x; i < (int)(n); i++)
#define rep(i, n) REP(i, 0, n)
#define all(x) (x).begin(), (x).end()
#define rall(x) (x).rbegin(), (x).rend()
#define F first
#define S second
#define mp make_pair
#define pb push_back

using namespace std;

typedef long long int lli;
typedef pair<int, int> Pii;

typedef vector<vector<int> > Matrix;

Matrix mul(const Matrix& A, const Matrix & B, const int mod) {
  Matrix C(A.size(), vector<int>(B[0].size()));
  
  for(int i = 0; i < A.size(); i++) {
    for(int k = 0; k < B.size(); k++) {
      for(int j = 0; j < B[0].size(); j++) {
        C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod;
      }
    }
  }

  return C;
}

// A^n
Matrix pow(Matrix A, int n, int mod) {
  Matrix B(A.size(), vector<int>(A.size()));
  
  for(int i = 0; i < A.size(); i++) {
    B[i][i] = 1;
  }

  while(n) {
    if(n & 1) B = mul(B, A, mod);
    A = mul(A, A, mod);
    n >>= 1;
  }

  return B;
}

int main() {
  int n, m, a, b, c, t;
  
  while(cin >> n >> m >> a >> b >> c >> t) {
    if((n | m | a | b | c | t) == 0) break;

    Matrix B(n, vector<int>(1));
    rep(i, n) cin >> B[i][0];

    Matrix A(n, vector<int>(n));
    rep(i, n) {
      if(i - 1 >= 0) A[i][i - 1] = a;
      A[i][i] = b;
      if(i + 1 < n) A[i][i + 1] = c;
    }

    A = pow(A, t, m);
    A = mul(A, B, m);
    
    rep(i, n) {
      cout << A[i][0] << (i + 1 == n ? '\n' : ' ');
    }
  }

  return 0;
}