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 点。

$10$ は合成数なので、掃き出し法のステップで詰まる部分がある
変数の個数が $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;
    }
}

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

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

TopCoder TCO 2013 Round 2A Medium - TheMagicMatrix

問題概要

制約

考えたこと

1. について

2. について

まとめ