AOJ 1327 One-Dimensional Cellular Automaton

行列を使う問題でした。nikeeshiくんががんばって考えました。
ありがとう。

#include <iostream>
#include <cstdio>
#include <iomanip>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <bitset>
#include <stack>
#include <utility>
#include <numeric>
#include <algorithm>
#include <functional>
#include <cctype>
#include <complex>
#include <string>
#include <sstream>

using namespace std;

#define all(c) c.begin(),c.end()
#define rall(c) c.rbegin(),c.rend()
#define mp(a,b) make_pair((a),(b))
#define eq ==

typedef long long ll;
typedef complex<double> point;
typedef pair<int,int> pii;

// →↑←↓
const int dx[] = {1,0,-1,0};
const int dy[] = {0,-1,0,1};


const double EPS = 1e-9;

namespace std{
    vector<vector<ll> > operator*(const vector<vector<ll> >& lhs,const vector<vector<ll> >& rhs){
        int N = lhs.size();
        vector<vector<ll> > ret(N,vector<ll>(N));
        for(int row=0;row<N;row++){
            for(int col=0;col<N;col++){
                int c = 0;
                for(int k=0;k<N;k++){
                    c+= rhs[row][k] * lhs[k][col];
                }
                ret[row][col] = c;
            }
        }
        return ret;
    }

    vector<vector<ll> > operator%(vector<vector<ll> > lhs,ll rhs){
        int N = lhs.size();
        vector<vector<ll> > ret = lhs;
        for(int i=0;i<N;i++){
            for(int j=0;j<N;j++){
                ret[i][j] = ret[i][j] % rhs;
            }
        }
        return ret;
    }
};

void display_matrix(vector<vector<ll> > mat){
    for(int i=0;i<mat.size();i++){
        for(int j=0;j<mat[0].size();j++){
            cerr << mat[i][j] << " ";
        }
        cerr << endl;
    }
}
void display_matrix(vector<ll> mat){
    for(int i=0;i<mat.size();i++){
        cerr << mat[i] << " ";
    }
    cerr << endl;
}

vector<vector<ll> >  mod_pow(vector<vector<ll> > x,ll n,ll mod){
    if(n==0){
        vector<vector<ll> > E(x.size(),vector<ll>(x.size()));
        for(int i=0;i<x.size();i++){
            E[i][i] = 1;
        }
        return E;
    }
    vector<vector<ll> > ret = mod_pow(x * x % mod,n/2,mod);
    if(n%2 == 1)  ret = ret * x % mod;
    return ret;
}

vector<ll> mat_multi(vector<vector<ll> > lhs,vector<ll> rhs,int mod){
    vector<ll> ret(rhs.size());
    int N = lhs.size();
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            ret[i] = (ret[i] + lhs[i][j] * rhs[j]) % mod;;
        }
    }
    return ret;
}

int main(){
    while(true){
        int N,M,A,B,C,T;
        cin >> N >> M >> A >> B >> C >> T;
        if(N == 0) break;
        vector<vector<ll> > matrix(N,vector<ll>(N));
        for(int i=0;i<N;i++){
            if(i != 0){
                matrix[i-1][i] = C;
            }
            matrix[i][i] = B;
            if(i != N-1){
                matrix[i+1][i] = A;
            }
        }
        vector<ll> S(N);
        for(int i=0;i<N;i++) cin >> S[i];
        vector<ll> ret = mat_multi(mod_pow(matrix,T,M),S,M);
        for(int i=0;i<N;i++){
            if(i != N-1){
                cout << ret[i] << " ";
            }else{
                cout << ret[i];
            }
        }
        cout << endl;
    }
    return 0;
}