亜種含めライツアウトの解法を求める（AOJ 1308） - 競技プログラミング+αなブログ

ICPC2010アジア予選の問題で、
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1308
というものがあった。
普通のライツアウトはある点を押すとその点とマンハッタン距離１の点が反転するが、この問題ではその点とマンハッタン距離dの点が反転する。
全ての点に対して押したかどうかを表す変数を盤面の数(w*h)だけ用意すると、その点を押すことでどの点が反転するかを考えることによって、mod2の合同式がw*h個できる。
行列で表せば、
Ax==b(mod 2)
である。ここで、bは与えられた盤面と対応する。
つまりライツアウトの問題を連立合同方程式の問題に帰着できるのである。
この問題では、与えられたライツアウトの盤面が解くことができるかできないかを求めればいいので、連立合同方程式が解を持つかどうかを判定すれば良い。
普通Ax=bを解くときはガウスの消去法で解くのだが、蟻本やネットを見てもどうも解が一意に求まる場合しか書いていない。
試行錯誤してもよく分からないので、線形代数の教科書を引っ張り出してきて読むと、
Ax=bが解を持つ必要十分条件は、
rank(A) = rank(A|b)
であり、
Ax=bが一意の解を持つ必要十分条件は、
rank(A) = rank(A|b) = n
と書いてあった。これが連立合同方程式のときも使うことができる。
rank(A)とrank(A|b)が求まればいいのだが、自分は次のように考えた。

手元にあるガウスの消去法ではピボットが見つからなかったとき一意の解無しと判定している
これを単純にcontinueするだけでは、上三角行列にはなるのだが、綺麗な階段状にならない
綺麗な階段状（階段行列）になるように工夫する（詳しくは割愛）
このとき、階段の高さがrank(A)である。
もしAに対応する部分では0だけの行（階段の下）で、bに対応する部分で0以外になっていれば、rank(A) != rank(A|b)
それ以外ではrank(A) = rank(A|b)

このアジア予選の問題を解くだけならこれだけで終わるのだが、どうせなら解xを求めてどこを押せばいいのかを知りたい。
解xの求め方を示す。
rank(A) = rank(A|b) = n なら普通の後退代入でいいのだが、nでないなら少し変わった後退代入をする。

まずx=vec(n,INF)とする（INFは任意を表す）。
上三角行列の下から見て行って、その行の一番左の０以外を見つけ、そこの変数の値を決定する。
このとき、計算のため必要となる変数で任意の値であるものは、0としてしまう。
最後にxの要素でINFがあれば0とする。

この方法はライツアウトのいかなる亜種にも対応できるはずで、例えば、斜め45度が全て反転する
http://www.gamedesign.jp/flash/tacoyaki/tacoyaki.html
でも解くことができた。

ソース

int modNorm(int a, int m) {
  return (a%m+m)%m;
}

// ax+by=gcd(a,b)
int extgcd(int a, int b, int &x, int &y) {
  int g = a;
  x = 1; y = 0;
  if (b) {
    g = extgcd(b, a%b, y, x);
    y -= (a/b) * x;
  }
  return g;
}

int invMod(int a, int p) {
  int x, y;
  if (extgcd(a,p,x,y) == 1) return (x+p)%p;
  return 0;
}

typedef vector<int> vec;
typedef vector<vec> mat;

bool modGaussElimination(const mat &A, const vec &b, int m, vec &res) {
  int n = A.size();
  mat B(n, vec(n+1));
  REP(i,n) REP(j,n)
    B[i][j] = A[i][j];
  REP(i, n) B[i][n] = b[i];
  
  int nowy = 0;
  REP(x, n) {
    int pivot = -1;
    for (int j=nowy; j<n; ++j)
      if (B[j][x]) {
        pivot = j;break;
      }
    if (pivot == -1) continue;
    swap(B[nowy], B[pivot]);

    for (int j=nowy+1; j<n; ++j) {
      int t = B[j][x] * invMod(B[nowy][x], m) % m;
      if (t)
        for (int k=x; k<=n; ++k)
          B[j][k] = modNorm(B[j][k] - B[nowy][k] * t, m);
    }
    nowy++;
    
  }

  res.clear();
  for (int y=nowy; y<n; ++y)
    if (B[y][n])                // rank(A) != rank(A|b)
      return 0;
  if (nowy != n) {              // rank(A) == rank(A|b) != n
    // 解が一意でない。ひとつ求める
    res = vec(n,INF);           // INFは任意を表す
    for (int y=n-1; y>=0; --y) {
      int x;
      for (x=y; x<n; ++x) {
        if (B[y][x])
          break;
      }
      if (x==n) continue;
      int sum = B[y][n];
      
      for (int i=x+1; i<n; ++i) {
        if (res[i] == INF) res[i] = 0; // この時点でres[i]が決まってなかったら0とする
        sum = modNorm(sum - res[i] * B[y][i], m);
      }
      res[x] = sum * invMod(B[y][x], m) % m;
    }
    REP(i, n) if (res[i]==INF) res[i] = 0;
    return 0;
  }
  // 解が一意に決まる。普通に後退代入
  res.resize(n);
  for (int x=n-1; x>=0; --x) {
    int sum = B[x][n];
    for (int i=n-1; i>x; --i) {
      sum = modNorm(sum - res[i] * B[x][i], m); 
    }
    res[x] = sum * invMod(B[x][x], m) % m;
  }
  
  return 1;
}

関数の解説

bool modGaussElimination(const mat &A, const vec &b, int m, vec &res)

Ax==b(mod m)を解く
resに解が格納される
返り値は、解が一意ならtrue, 一意でないならfalse

ライツアウトの問題を解くだけならmod 2限定なのだが、せっかくなので他のmodにも対応できるようにした。そのために逆元を計算するinvModなどが必要になった。
(A|b)をまとめてBとしている。
実際の使い方は、別記事(d:id:sune2:20120127:1327638017)参照。