hdu 5755 Gambler Bo 高斯消元 + 取余逆元

分析：
对每一个位置设一个未知变量$x$，每个位置都有一个结果变量$y$，表示要操作多少次可以把该位置变为$0$，这样对于每一个未知量可以对其周围的元素产生影响，列出一个现象方程组

$$MX=Y$$
$M$是系数矩阵，需要去根据每个元素的影响区构造。
高斯消元求解即可。
注意:这里的一切都是在模3,剩余系下的代数系统，所以每一步都要取余，出发要求逆元。
高斯消元复杂度是立方级别的，所以总的复杂度$O(N^3M^3)$

高斯消元模板

const int maxn = 100;
int a[maxn][maxn], X[maxn];//分别记录增广矩阵和解集
int free_x[maxn];//记录自由变量

int LCM(int x, int y){
    return x / __gcd(x, y) *y;
}

// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，但无整数解，-1表示无解，0表示唯一解，大于0表示无穷解，并返回自由变元的个数)
int Guass(int equ, int var) { //分别表示方程组的个数和变量的个数: n
    int i, j, k, col;
    memset(X, 0, sizeof(X));
    memset(free_x, 1, sizeof(free_x));
    for (k = 0,col = 0; k < equ && col < var; ++k, ++col){//枚举行列
        int max_r = k;//找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
        for (i = k + 1; i < equ; ++i) if (abs(a[i][col]) > abs(a[max_r][col])) max_r = i;
        if (max_r != k) for (i = k; i < var + 1; ++i) swap(a[k][i],a[max_r][i]);
        if (a[k][col] == 0) k--;//如果对应该列都为0，枚举该行的下一列  
        else { 
            for (i = k + 1; i < equ; ++i){//将k后边的col进行初等变换成行阶梯矩阵
                if (a[i][col] != 0){
                    int lcm = LCM(a[k][col], a[i][col]);
                    int ta = lcm / abs(a[i][col]), tb = lcm / abs(a[k][col]);
                    if (a[i][col] * a[k][col] < 0) tb = -tb;
                    for (j = col; j < var + 1; ++j){
                        //如果是在模剩余系的方程和你要取模
                        a[i][j] = ta*a[i][j] - tb*a[k][j];
                    }
                }
            }
        }
    }
    // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0). 即R(A) != R(A')无解
    for (i = k; i < equ; ++i){
        if (a[i][col] != 0) return -1;
    }
    // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.
    // 且出现的行数即为自由变元的个数.   即R(A) = R(A') < n
    if (k < var){
        // 首先，自由变元有var - k个，即不确定的变元至少有var - k个.
        int num = 0,freeidx;
        for (i = k - 1; i >= 0; --i){
            num = 0;// 用于判断该行中的不确定的变元的个数，如果超过1个，则无法求解，它们仍然为不确定的变元.
            int tmp = a[i][var];
            // 第i行一定不会是(0, 0, ..., 0)的情况，因为这样的行是在第k行到第equ行.
            // 同样，第i行一定不会是(0, 0, ..., a), a != 0的情况，这样的无解的.
            for (j = 0; j < var; ++j){
                if (a[i][j] != 0 && free_x[j]){
                    num++;
                    freeidx = j;
                }
            }
            if (num > 1) continue; // 无法求解出确定的变元.
            // 说明就只有一个不确定的变元free_index，那么可以求解出该变元，且该变元是确定的.
            tmp = a[i][var];
            for (j = 0; j < var; ++j){
                if (a[i][j] && j != freeidx) tmp -= a[i][j]*X[j];
            }
            X[freeidx] = tmp / a[i][freeidx];
            free_x[freeidx] = 0;
        }
        return var - k;
    }
    // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
    // 计算出Xn-1, Xn-2 ... X0.
    for (i = k - 1; i >= 0; --i){
        int tmp = a[i][var];
        for (j = i + 1; j < var; ++j){
            tmp -= a[i][j] * X[j];
        }
        //模剩余系求逆元
        X[i] = tmp / a[i][i];
    }
    return 0;
}

本题代码：

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define pr(x) cout << #x << ": " << x << "  " 
#define pl(x) cout << #x << ": " << x << endl;

struct jibancanyang
{
    int n, m;
    int A[909][909], X[909], equ, var;

    int lcm(int x, int y) {
        return x / __gcd(x, y) * y;
    }

    int extGcd(int a, int b, int &x, int &y) {
        int d = a;
        if (b) d = extGcd(b, a % b, y, x), y -= a / b * x;
        else x = 1, y = 0;
        return d;
    }

    int Guass() {
        int i, j, k, col;
        memset(X, 0, sizeof (X));
        for (k = 0, col = 0; k < equ && col < var; ++k, ++col) {
            int max_r = k;
            for (i = k + 1; i < equ && A[max_r][col] < 2; ++i) if (abs(A[i][col]) > abs(A[max_r][col])) max_r = i;
            if (max_r != k) for (i = k; i < var + 1; ++i) swap(A[k][i], A[max_r][i]);
            if (A[k][col] == 0)  k--;
            else {
                for (i = k + 1; i < equ; ++i) {
                    if (A[i][col]) {
                        int ta = A[i][col], tb = A[k][col];
                        for (j = col ; j < var + 1; ++j) {
                            A[i][j] = ((tb * A[i][j] - ta * A[k][j]) % 3 + 3) % 3;
                        }
                    }
                }
            }
        }
        for (i = k - 1; i >= 0; --i) {
            int tmp = A[i][var];
            for (j = i + 1; j < var; ++j){
                tmp =  ((tmp - A[i][j] * X[j]) % 3 + 3) % 3;
            }
            int x, y;
            extGcd(A[i][i], 3, x, y);
            X[i] = (x % 3 + 3) % 3 * tmp % 3;
        }
        return 0;
    }

    int getPla(int i, int j) {
        return i * m + j;
    }

    void fun() {
        int T;
        scanf("%d", &T);
        while (T--) {
            scanf("%d%d", &n, &m);
            memset(A, 0, sizeof(A));
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < m; j++) {
                    int x, pla = getPla(i, j);
                    scanf("%d", &x);
                    A[pla][pla] = 2;
                    if (i - 1 >= 0) A[getPla(i - 1, j)][pla] = 1;
                    if (i + 1 < n) A[getPla(i + 1, j)][pla] = 1;
                    if (j - 1 >= 0) A[getPla(i, j - 1)][pla] = 1;
                    if (j + 1 < m) A[getPla(i , j + 1)][pla] = 1;
                    A[pla][n * m] = (3 - x) % 3;
                }
            }
            equ = n * m, var = n * m;
            Guass();
            int ans = 0;
            for (int i = 0; i < n; ++i) {
                for (int j = 0; j < m; ++j) {
                    int pla = getPla(i, j);
                    ans += (X[pla] =  (X[pla] + 3) % 3);
                }
            }
            printf("%d\n", ans);
            for (int i = 0; i < n; ++i) {
                for (int j = 0; j < m; ++j) {
                    for (int k = X[getPla(i, j)]; k > 0; --k) {
                        printf("%d %d\n", i + 1, j + 1);
                    }
                }
            }
        }
    }

}ac;

int main()
{
#ifdef LOCAL
    freopen("in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
#endif
    ac.fun();
    return 0;
}