AOJ 2624 : Graph Automata Player

問題リンク : Graph Automata Player | Aizu Online Judge

問題概要 :
有向グラフが与えられる ( 多重辺なし、自己ループあり )
有向グラフの各ノードには０か１が書かれている
１秒毎に各ノードの数値を更新する
更新の手順は以下のとおり

現在の時刻をtとする
各ノードについて、そのノードと直接有向辺で結ばれているノードに(t-1)秒の時点で書かれている数値の合計を求め、それを２で割った余りを書き込む

T秒後の有向グラフが与えられる
有向グラフの初期状態を求めよ
ただし、答えが一意に定まらない場合は "ambiguous"
答えがない場合は "none" と出力せよ

解法 :
入力で与えられる辺の状況を表す隣接行列をA、各ノードの数値の情報を行列bとする
t秒の各ノードの数値の情報を持った行列をx_tとすると
A * x_0 = x_1
A * x_1 = x_2
...
A * x_(T-1) = x_T

このx_Tが入力で与えられているbなので
初期状態x_0を求めるためには左からAの逆行列A_-1を掛けると良い

x_(T-1) = A_-1 * x_T
...
x_1 = A_-1 * x_2
x_0 = A_-1 * x_1

x_0 = A_-1 * ( A_-1 * ( A_-1 * .... A_-1 * ( A_-1 * x_T ) ... ) )
= (A_-1)^T * b

これより、 A_-1のT乗を求めれば良いことがわかる
A_-1を求める際に、A_-1が作れない場合であったり一意に定まらない場合は"none"や"ambiguous"を出力して終える

Aの逆行列が作れない => Aのランク!= A|b のランク
解が一意に定まらない=> Aのランク == A|b のランク == n ( n は A の行数 or 列数 )
Aのランクと A|b のランクが等しい時解をもつが、それが n と等しくない場合はなんとでもできる変数があるので解が一意に定まらない
詳しくは線形代数の教科書で ( ませまの薄いやつの後ろのほうにのってたはず

行列の掛け算等はライブラリ(たしかSpagettie source

コード:

#include<bits/stdc++.h>
 
#define REP(i,s,n) for(int i=s;i<n;i++)
#define rep(i,n) REP(i,0,n)
#define EPS (1e-8)
#define equals(a,b) (fabs((a)-(b))<EPS)
 
using namespace std;
 
const int IINF = INT_MAX;
 
int extgcd(int a,int b,int& x,int& y){
  int d = a;
  if(b != 0){
    d = extgcd(b,a%b,y,x);
    y -= (a/b)*x;
  } else x = 1,y = 0;
  return d;
}
 
int mod_inv(int a,int m) {
  int x,y;
  extgcd(a,m,x,y);
  return (m+x%m)%m;
}
 
vector<vector<int> > identity(int n) {
  vector<vector<int> > A(n,vector<int>(n,0));
  for(int i=0;i<n;i++) A[i][i] = 1;
  return A;
}

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

void mod_pow(vector<vector<int> > &A,int e, int mod) {
  vector<vector<int> > C = identity(A.size());
  while( e ) {
    if( e & 1 ) mod_multi(C, A, mod);
    mod_multi(A, A, mod);
    e >>= 1;
  }
  A = C;
}

bool mod_inv_matrix(vector<vector<int> > &A,int m,vector<vector<int> > &res){
  int n = A.size();
  vector<vector<int> > B(n,vector<int>(n+n));
  rep(i,n) rep(j,n) B[i][j] = A[i][j], B[i][n+j] = ( i == j );
  rep(x,n){
    int pivot = -1;
    REP(j,x,n) if( B[j][x] ) { pivot = j; break; }
    //if( pivot == -1 ) continue;
    if( pivot == -1 ) return false;
    swap(B[x],B[pivot]);
    rep(j,n) if( j != x ) {
      int t = B[j][x] * mod_inv(B[x][x],m) % m;
      if( t ) REP(k,x,n+n) B[j][k] = ( ( B[j][k] - B[x][k] * t ) % m + m ) % m;
    }
    int t = mod_inv(B[x][x],m);
    REP(j,x,n+n) B[x][j] = ( B[x][j] * t ) % m;
  }
 
  res.clear();
  res.resize(n,vector<int>(n));
  rep(i,n) rep(j,n) res[i][j] = B[i][n+j];
  return true;
}
 
int rank(vector<vector<int> > A,int mod) {
  int H = A.size(), W = A[0].size();
  int r = 0;
  rep(x,W) {  
    int pivot = IINF;
    REP(y,r,H) if(A[y][x]) {
      pivot = y;
      break;
    }
    if( pivot == IINF )continue;
    swap(A[r],A[pivot]);
    REP(y,r+1,H) {
      int value = A[y][x] * mod_inv(A[r][x],mod) % mod;
      if(value) REP(x2,x,W) A[y][x2] = ((A[y][x2] - A[r][x2] * value)%mod + mod ) % mod; 
    }
    r++;
  }
  return r;
}
 
int main(){
  int n,phase;
  scanf("%d",&n);
  vector<vector<int> > A(n,vector<int>(n));
  vector<int> b(n);
  rep(i,n) rep(j,n) scanf("%d",&A[i][j]);
  rep(i,n) scanf("%d",&b[i]);
  scanf("%d",&phase);
 
  mod_pow(A,phase,2);

  vector<vector<int> > Ab(n,vector<int>(n+1));
  rep(i,n) rep(j,n) Ab[i][j] = A[i][j];
  rep(i,n) Ab[i][n] = b[i];
  int rankA  = rank(A,2);
  int rankAb = rank(Ab,2);

  if( rankA != rankAb ) { puts("none"); return 0; }
  if( rankA != n      ) { puts("ambiguous"); return 0; }
 
  vector<vector<int> > invA;
  mod_inv_matrix(A,2,invA);
  vector<int> x = mod_multi(invA,b,2);
 
  rep(i,x.size()) {
    if( i ) printf(" ");
    cout << x[i];
  } puts("");
 
  return 0;
}