けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

TCO 2013 Round 2A Med TheMagicMatrix

mod. p での連立一次方程式の練習などに!!!

問題へのリンク

問題概要

 N ×  N のマス目があり、そのうちの  M マスについては既に  0 以上  9 以下の数値が入れられている。残りの  N^{2} - M マスに  0 以上  9 以下の数値を入れる方法のうち

  • どの行の和と、どの列の和をとっても、その一の位が等しい

を満たすものの個数を  1234567891 で割った余りを求めよ。

制約

  •  1 \le N \le 1000
  •  1 \le M \le 10

考えたこと

いかにも  {\rm mod}. 10 における連立方程式が立てられそうな問題ではある!!!

実際、各マス  (i, j) の値を  x_{i, j} として、各行や各列の和の一の位の値を  X とすると ( X = 0, 1, \dots, 9 についてすべて試す感じ)、

  • 各行の和が  {\rm mod}. 10 X であるという条件は

 \sum_{j = 1}^{N} x_{i, j} ≡ X

  • 各列の和が  {\rm mod}. 10 X であるという条件は

 \sum_{i = 1}^{N} x_{i, j} ≡ X

  • 指定された  M マスに関する条件は

 x_{i, j} ≡ a_{i, j}

といった具合に自然に連立一次方程式を立てることができる。問題となるのは以下の 2 点。

  1.  10合成数なので、掃き出し法のステップで詰まる部分がある
  2. 変数の個数が  N^{2} 個となり、 N = 1000 ではとても扱えるサイズではない

1. について

連立一次方程式  Ax ≡ b ({\rm mod}. 10) を解くためにまず

  •  Ax ≡ b ({\rm mod}. 2)
  •  Ax ≡ b ({\rm mod}. 5)

をそれぞれ解く。そして中国剰余定理により、


 Ax ≡ b ({\rm mod}. 10) の解  x と、

  •  Ax ≡ b ({\rm mod}. 2) の解  x_1
  •  Ax ≡ b ({\rm mod}. 5) の解  x_2

とから定まる組  (x_1, x_2) とは一対一対応する


ということがわかるので、結局

  •  Ax ≡ b ({\rm mod}. 2) の解の個数
  •  Ax ≡ b ({\rm mod}. 5) の解の個数

の積が答えになることがわかる。

2. について

少し考えてみると、「その行に数字が 1 個もない」という行と、「その列に数字が 1 個もない」という列の両方を含むような盤面が与えられたならば、

  • その行または列以外のすべてのマスを適当に数値を入れる

としたときに、各行や各列の和となる値  X を指定すると、残りの  2N-1 マスについてはただ一つに決まることがわかる。

よって、そのような場合については

  •  10^{(N-1)^{2} - M + 1} 通り

と求めることができる。

まとめ

以上より、

  • 数字が一つもない行と列があるときはすぐに求められる
  • そうでない場合は、 N \le M であることが保証されるので変数の個数は  O(M^{2}) 個であり、 {\rm mod}. 2 {\rm mod}. 5 でそれぞれ連立方程式の解の個数を求めればよく、計算量は  O(M^{6}) になる

という感じで求めることができる。

#include <iostream>
#include <vector>
#include <set>
#include <algorithm>
using namespace std;


const long long MOD = 1234567891LL;
long long modinv(long long a, long long mod) {
    long long b = mod, u = 1, v = 0;
    while (b) {
        long long t = a/b;
        a -= t*b; swap(a, b);
        u -= t*v; swap(u, v);
    }
    u %= mod;
    if (u < 0) u += mod;
    return u;
}

// matrix
template<int MOD> struct Matrix {
    vector<vector<long long> > val;
    Matrix(int n, int m, long long x = 0) : val(n, vector<long long>(m, x)) {}
    void init(int n, int m, long long x = 0) {val.assign(n, vector<long long>(m, x));}
    size_t size() const {return val.size();}
    inline vector<long long>& operator [] (int i) {return val[i];}
};

template<int MOD> ostream& operator << (ostream& s, Matrix<MOD> A) {
    s << endl; 
    for (int i = 0; i < A.size(); ++i) {
        for (int j = 0; j < A[i].size(); ++j) {
            s << A[i][j] << ", ";
        }
        s << endl;
    }
    return s;
}

template<int MOD> Matrix<MOD> operator * (Matrix<MOD> A, Matrix<MOD> B) {
    Matrix<MOD> R(A.size(), B[0].size());
    for (int i = 0; i < A.size(); ++i) 
        for (int j = 0; j < B[0].size(); ++j)
            for (int k = 0; k < B.size(); ++k) 
                R[i][j] = (R[i][j] + A[i][k] * B[k][j] % MOD) % MOD; 
    return R;
}

template<int MOD> Matrix<MOD> pow(Matrix<MOD> A, long long n) {
    Matrix<MOD> R(A.size(), A.size());
    for (int i = 0; i < A.size(); ++i) R[i][i] = 1;
    while (n > 0) {
        if (n & 1) R = R * A;
        A = A * A;
        n >>= 1;
    }
    return R;
}

template<int MOD> int GaussJordan(Matrix<MOD> &A, bool is_extended = false) {
    int m = A.size(), n = A[0].size();
    for (int row = 0; row < m; ++row)
        for (int col = 0; col < n; ++col)
            A[row][col] = (A[row][col] % MOD + MOD) % MOD;
    int rank = 0;
    for (int col = 0; col < n; ++col) {
        if (is_extended && col == n-1) break;
        int pivot = -1;
        for (int row = rank; row < m; ++row) {
            if (A[row][col] != 0) {
                pivot = row;
                break;
            }
        }
        if (pivot == -1) continue;
        swap(A[pivot], A[rank]);
        auto inv = modinv(A[rank][col], MOD);
        for (int col2 = 0; col2 < n; ++col2)
            A[rank][col2] = A[rank][col2] * inv % MOD;
        for (int row = 0; row < m; ++row) {
            if (row != rank && A[row][col]) {
                auto fac = A[row][col];
                for (int col2 = 0; col2 < n; ++col2) {
                    A[row][col2] -= A[rank][col2] * fac % MOD;
                    if (A[row][col2] < 0) A[row][col2] += MOD;
                }
            }
        }
        ++rank;
    }
    return rank;
}

template<int MOD> int linear_equation(Matrix<MOD> A, vector<long long> b, vector<long long> &res) {
    int m = A.size(), n = A[0].size();
    Matrix<MOD> M(m, n + 1);
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) M[i][j] = A[i][j];
        M[i][n] = b[i];
    }
    int rank = GaussJordan(M, true);

    // check if it has no solution
    for (int row = rank; row < m; ++row) if (M[row][n]) return -1;

    // answer
    res.assign(n, 0);
    for (int i = 0; i < rank; ++i) res[i] = M[i][n];
    return rank;
}



long long modpow(long long a, long long n, long long mod) {
    long long res = 1;
    while (n > 0) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

class TheMagicMatrix {
public:
    int find(int n, vector <int> rows, vector <int> cols, vector <int> vals) {
        int m = rows.size();

        // 数字のない行と列がある場合
        set<int> sr, sc;
        for (int i = 0; i < rows.size(); ++i) {
            sr.insert(rows[i]);
            sc.insert(cols[i]);
        }
        if (sr.size() < n && sc.size() < n) {
            long long ex = (n-1)*(n-1) - m + 1;
            return modpow(10LL, ex, MOD);
        }

        // 連立一次方程式を立てる
        Matrix<2> M2(n*2+m, n*n);
        Matrix<5> M5(n*2+m, n*n);
        vector<long long> b(n*2+m);

        // 行和
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                int id = i * n + j;
                M2[i][id] = M5[i][id] = 1;
            }
        }

        // 列和
        for (int j = 0; j < n; ++j) {
            for (int i = 0; i < n; ++i) {
                int id = i * n + j;
                M2[j + n][id] = M5[j + n][id] = 1;
            }
        }

        // 条件
        for (int k = 0; k < m; ++k) {
            int id = rows[k] * n + cols[k];
            M2[k + n*2][id] = M5[k + n*2][id] = 1;
            b[k + n*2] = vals[k];
        }

        // X = 0, 1, ..., 9
        long long res = 0;
        for (int X = 0; X < 10; ++X) {
            for (int i = 0; i < n*2; ++i) b[i] = X;

            vector<long long> v;
            int rank2 = linear_equation(M2, b, v);
            int rank5 = linear_equation(M5, b, v);
            if (rank2 == -1 || rank5 == -1) continue;
            long long tmp = modpow(2LL, n*n-rank2, MOD) * modpow(5LL, n*n-rank5, MOD) % MOD;
            res = (res + tmp) % MOD;
        }
        return res;
    }
};


// for test
int main() {
    int n, m;
    while (cin >> n >> m) {
        vector<int> r(m), c(m), v(m);
        for (int i = 0; i < m; ++i) cin >> r[i] >> c[i] >> v[i];
        TheMagicMatrix tmm;
        cout << tmm.find(n, r, c, v) << endl;
    }
}