1308 - Awkward Lights

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1308

問題

ライツアウトで解法があるかどうか調べる。
普通のライツアウトはある点を押すとその点とマンハッタン距離１の点が反転するが、この問題ではその点とマンハッタン距離dの点が反転する。
解法があるなら1、無いなら0を出力せよ。

解法

連立合同方程式。詳しくは、 d:id:sune2:20120127:1327637624
下のソースは一般の連立合同方程式を解くためにいらんことしている。

ソース

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;
}

int main() {
  int m, n, d;
  while(cin>>m>>n>>d, m||n||d) {
    bool ba[n][m];
    mat A(n*m, vec(n*m));
    vec b(n*m);
    REP(i,n) {
      REP(j,m) {
        cin >> ba[i][j];
        A[i*m+j][i*m+j] = 1;
        REP(k,n) {
          REP(l,m) {
            int dis = abs(i-k)+abs(j-l);
            if (dis == d) {
              A[k*m+l][i*m+j] = 1;
            }
          }
        }
        b[i*m+j] = ba[i][j];
      }
    }
    vec x;
    modGaussElimination(A,b,2,x);
    cout << (x.size() != 0) << endl;
  }
}