TopCoder SRM 620 Div1 Hard PerfectSquare

問題

nxn行列が与えられる。
この行列から以下の条件を満たすようにいくつか要素を選ぶ。
・それぞれの行の中に選んだ要素が奇数個ある
・それぞれの列の中に選んだ要素が奇数個ある
・選んだ要素の積が平方数である

要素の選び方は何通りあるか、mod 10^9 + 7で答えよ。

制約条件

n≦20
A[i][j]≦10^9

方針

mod2で以下の連立方程式を立てる。
要素(i, j)を選ぶことをxij = 1, 選ばないことをxij = 0と表現する変数xijを作る（0≦i, j＜n）

1番目の条件は、各行iについて、Σ{xij | j } ≡ 1 (mod2)
2番目の条件は、各列jについて、Σ{xij | i } ≡ 1 (mod2)
3番目の条件は、
A[i][j]を素因数分解したときに出てくる素因数全てを小さい順に並べた配列をp[k]とすると、
素因数p[k]について, Σ{xij | xij % p[k] == 0} ≡ 0 (mod2)

この連立方程式を吐き出し法で解く。
解がなかったら答えは0通り。
解がある場合、ランクrが変数n*nより少なかったら、その分だけ解の自由度が2通り増えるので、
2^(n*n - r)通りが答えになる。

ソースコード

typedef vector<vi> M;
int rank(vector<vi> A){
	
	int h = A.size(), w = A[0].size();
	int pos = 0;
	
	rep(i, h){
		int p = -1;
		
		while(p < 0 && pos < w - 1){
			for(int j = i; j < h; j++) if(A[j][pos]){
				p = j;
				goto OUT;
			}
			pos++;
		}
		OUT:
		
		if(p < 0){
			for(int j = i; j < h; j++) if(A[j][w - 1]) return -1;
			return i;
		}
		
		swap(A[p], A[i]);
		for(int j = i + 1; j < h; j++)if(A[j][pos])
			for(int k = pos; k < w; k++) A[j][k] ^= A[i][k];
		pos++;
	}
	return h;
}

class PerfectSquare {
	public:
	int ways(vector <int> x) {
		int n = x.size();
		n = sqrt(n);
		
		M A;
		map<int, int> row;
		
		rep(i, n){
			vi u, v;
			rep(j, n) rep(k, n){
				v.pb(i == j ? 1 : 0);
				u.pb(i == k ? 1 : 0);
			}
			u.pb(1); v.pb(1);
			A.pb(u); A.pb(v);
		}
		rep(i, n * n){
			int a = x[i];
			for(int j = 2; j * j <= a; j++) if(a % j == 0){
				
				int cnt = 0;
				while(a % j == 0) a /= j, cnt ^= 1;
				if(!cnt) continue;
				
				if(!row.count(j)){
					row[j] = A.size();
					A.pb(vi(n * n + 1));
				}
				A[row[j]][i] ^= 1;
			}
			if(a == 1) continue;
			
			if(!row.count(a)){
				row[a] = A.size();
				A.pb(vi(n * n + 1));
			}
			A[row[a]][i] ^= 1;
		}
		#if 0
		each(i, row) dbg(i->first);
		rep(i, n) rep(j, n) cerr<<x[i*n+j]<<(j==n-1?"\n":" ");
		rep(i, A.size()) rep(j, A[i].size()) cerr<<A[i][j]<<(j==A[i].size()-1?"\n":" ");
		#endif 
		
		int r = ::rank(A);
		int ans = 1;
		
		if(r < 0) return 0;
		rep(i, n * n - r) ans *= 2, ans %= (int)1e9 + 7;
		return ans;
	}
};