[BZOJ 2854] civilization 高斯消元

题意

　　给定 $n$, $A_{n\times n}, B = \left\{ b_i \right\}$ .

　　解方程 $Ax = B^T$ .

　　$n \le 200, A_{ij} \in [- {10} ^ 9, {10} ^ 9]$ .

　　保证 $x_i \in [- {10} ^ {18}, {10} ^ {18}]$ .

分析

　　大素数取模下高斯消元, 通过 CRT 进行合并.

　　http://blog.csdn.net/owen_hzt/article/details/41493637

实现

　　首先, 我们测试一下 long long 大致可以存到多少.

　　long long 的最大值大约是 $9 \times {10} ^ {18}$ .

　　我们要选取两个合适的大素数取模, 然后进行 CRT .

　　保证这两个素数大于 ${10} ^ 9$ , 相乘大于 ${10} ^ {18}$ , 且不会超过 long long 的最大值.

　　而且, 这两个素数要保证系数矩阵满秩.

　　系数矩阵是否满秩, 必须要通过高斯消元来发现, 所以我们先准备多几个素数, 取能成功的前两个.

　　我们可以写一个判断是否是素数的程序, 进行检验.

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>

inline bool Prime(int x) {
    if (x <= 1) return false;
    for (int i = 2; i * i <= x; i++)
        if (x % i == 0) return false;
    return true;
}

int main(void) {
    for (int x; ~scanf("%d", &x); )
        puts(Prime(x) ? "YES" : "NO");
    return 0;
}

　　我们找到了这样几个数: 1000000007, 1000000009, 1000000021, 1000000087.

　　最后合并答案的时候注意要写快速乘法.

　　因为对于 mod 一个 long long 大小的数的操作, 直接乘会爆 long long !

　　最终代码:

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
using namespace std;

#define F(i, a, b) for (register int i = (a); i <= (b); i++)

#define LL long long

const int N = 205;
const LL P[4] = {(int)1e9+7, (int)1e9+9, (int)1e9+21, (int)1e9+87};

int n; LL A[N][N], B[N][4];
LL X[N][4]; int z1, z2;

inline LL rd(void) {
    LL f = 1; char c = getchar(); for (; !isdigit(c); c = getchar()) if (c == ‘-‘) f = -1;
    LL x = 0; for (; isdigit(c); c = getchar()) x = x*10+c-‘0‘; return x*f;
}

LL Inv(LL x, LL M) { LL res = 1; for (LL y = M-2; y > 0; y >>= 1, x = x * x % M) if (y & 1) res = res * x % M; return res; }
bool Gauss(int K) {
    LL M = P[K];
    static LL a[N][N];
    F(i, 1, n) {
        F(j, 1, n) a[i][j] = A[i][j];
        a[i][n+1] = B[i][K];
    }

    F(i, 1, n) {
        if (!a[i][i]) {
            bool G = false;
            for (int j = i+1; j <= n && !G; j++) if (a[j][i] > 0) {
                G = true;
                F(k, 1, n+1) swap(a[i][k], a[j][k]);
            }
            if (!G) return false;
        }
        LL R = Inv(a[i][i], M);
        F(j, 1, n) if (i != j && a[j][i] > 0) {
            LL tmp = R * a[j][i] % M;
            F(k, 1, n+1)
                a[j][k] = ((a[j][k] - a[i][k] * tmp) % M + M) % M;
        }
    }

    F(i, 1, n) {
        LL R = Inv(a[i][i], M);
        X[i][K] = a[i][n+1] * R % M;
    }
    return true;
}

inline LL Mul(LL x, LL y, LL P) {
    LL res = 0;
    for (; y > 0; y >>= 1, x = (x + x) % P)
        if (y & 1) res = (res + x) % P;
    return res;
}
inline void Print(LL X1, LL P1, LL X2, LL P2) {
    printf("%lld\n", (Mul(X1 * P2, Inv(P2 % P1, P1), P1 * P2) + Mul(X2 * P1, Inv(P1 % P2, P2), P1 * P2)) % (P1 * P2));
}

int main(void) {
    #ifndef ONLINE_JUDGE
        freopen("bzoj2854.in", "r", stdin);
        freopen("bzoj2854.out", "w", stdout);
    #endif

    n = rd();
    F(i, 1, n) {
        F(j, 1, n) A[i][j] = rd();
        static char s[50]; scanf("%s", s+1);
        for (int j = 0, L = strlen(s+1); j < 4; j++)
            F(k, 1, L) B[i][j] = (B[i][j] * 10 + s[k] - ‘0‘) % P[j];
    }

    z1 = z2 = -1;
    for (int k = 0; k < 4; k++) {
        if (!Gauss(k)) continue;
        if (z1 == -1) z1 = k;
        else if (z2 == -1) z2 = k;
    }

    F(i, 1, n)
        Print(X[i][z1], P[z1], X[i][z2], P[z2]);

    return 0;
}