hdu 5755（高斯消元——模线性方程组模板）

PS. 看了大神的题解，发现确实可以用m个未知数的高斯消元做。因为确定了第一行的情况，之后所有行的情况都可以根据第一行推。 这样复杂度直接变成O(m*m*m)

 

知道了是高斯消元后，其实只要稍加处理，就可以解决带模的情况。

1 是在进行矩阵行变化的时候，取模。

2 最后的除法用逆元。（因为a[i][i]必定非0 且小于模数） 

然后对于无穷多解的情况，只需要将那些列全为0的未知数定义一个固定值。（这里设的是0）其余操作不变。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <string>

#define CL(a,num) memset((a),(num),sizeof(a))
#define iabs(x)  ((x) > 0 ? (x) : -(x))
#define Min(a,b) (a) > (b)? (b):(a)
#define Max(a,b) (a) > (b)? (a):(b)

#define ll long long
#define inf 0x7f7f7f7f
#define MOD 100000007
#define lc l,m,rt<<1
#define rc m + 1,r,rt<<1|1
#define pi acos(-1.0)
#define test puts("<------------------->")
#define maxn 100007
#define M 100007
#define N 1010
using namespace std;
//freopen("din.txt","r",stdin);

int a[N][N],X[N];//分别记录增广矩阵和解集
int free_x[N];//记录自由变量
int equ,var;//分别表示方程组的个数和变量的个数

int GCD(int x,int y){
    if (y == 0) return x;
    return GCD(y,x%y);
}
int LCM(int x,int y){
    return x/GCD(x,y)*y;
}
void Debug(void)
{
    int i, j;
    for (i = 0; i < equ; i++)
    {
        for (j = 0; j < var + 1; j++)
        {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}
// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，
//但无整数解,-1表无整数解,-1表示无解，0表示唯一解，大于0表示无穷解，
//并返回自由变元的个数)

//ax + by = gcd(a,b)
//传入固定值a,b.放回 d=gcd(a,b), x , y
void extendgcd(ll a,ll b,ll &d,ll &x,ll &y)
{
    if(b==0){d=a;x=1;y=0;return;}
    extendgcd(b,a%b,d,y,x);
    y-=x*(a/b);
}

//Ax=1(mod M)，gcd(A,M)==1
//输入：10^18>=A,M>=1
//输出：返回x的范围是[1,M-1]
ll GetNi(ll A,ll MM)
{
    ll rex=0,rey=0;
    ll td=0;
    extendgcd(A,MM,td,rex,rey);
    return (rex%MM+MM)%MM;
}

int Guass(){
    int i,j,k,col;
    CL(X,0); CL(free_x,1);
    
    for (k = 0,col = 0; k < equ && col < var; k++, col++){
        int max_r = k;
        for (i = k + 1; i < equ; ++i){
            if (iabs(a[i][col]) > iabs(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){
            //可以随意定义的变量
            X[col] = 0;//强制赋值为0
            free_x[col] = 0;
            k--;
            //cout<<k<<endl;
            continue;
        }
        for (i = k + 1; i < equ; ++i){
            if (a[i][col] != 0){
                int lcm = LCM(a[k][col],a[i][col]);
                int ta = lcm/iabs(a[i][col]); int tb = lcm/iabs(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])%3+3)%3;
                }
            }
        }
    }
    //Debug();
    // 1. 无解的情况:
    for (i = k; i < equ; ++i){
        if (a[i][col] != 0) return -1;
    }
    // 2. 无穷解的情况:
    if (k < var){
        
        int num = 0,freeidx=0;
        for (i = k - 1; i >= 0; --i){
            num = 0;
            int tmp = a[i][var];
            for (j = 0; j < var; ++j){
                if (a[i][j] != 0 && free_x[j]){
                    num++;
                    freeidx = j;
                }
            }
            if (num > 1) continue;
            tmp = a[i][var];
            for (j = 0; j < var; ++j){
                if (a[i][j] && j != freeidx) tmp -= a[i][j]*X[j];
            }
            //这里也要
            
            int k2 = (tmp%3+3)%3;
            int k1 = (a[i][freeidx]%3+3)%3;
            if(k1!=0)
            {
                X[freeidx] = k2*(int)GetNi(k1, 3);
            }
            else
            {
                X[freeidx] = 0;
                //printf("X[%d]为任意？\n",i);
            }
            
            //X[freeidx] = tmp/a[i][freeidx];
            free_x[freeidx] = 0;
        }
        return var - k;
    }
    // 3. 唯一解的情况:
    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];
        }
        //强行搞一发
        int k2 = (tmp%3+3)%3;
        int k1 = (a[i][i]%3+3)%3;
        if(k1!=0)
        {
            X[i] = k2*(int)GetNi(k1, 3);
        }
        else
        {
            X[i] = 0;
            //printf("X[%d]为任意？\n",i);
        }
        //X[i] = tmp/a[i][i];//不整除？
    }
    return 0;
}


int g[33][33];

int getid(int i,int j,int n,int m)
{
    return (i-1)*m + (j-1);
}
int up[4]={1,-1,0,0};
int rl[4]={0,0,1,-1};

int main(){
    //freopen("din.txt","r",stdin);
    int i,j;
    
    int T;
    cin>>T;
    while(T--)
    {
        int n,m;
        scanf("%d%d",&n,&m);
        for(i=1;i<=n;i++)
            for(j=1;j<=m;j++)
                scanf("%d",&g[i][j]);
        equ = n*m;
        var = n*m;
        
        memset(a,0,sizeof(a));
        
        int id = 0;
        for(i=1;i<=n;i++)
            for(j=1;j<=m;j++)
            {
                for(int k=0;k<4;k++)
                {
                    int ti = i+up[k];
                    int tj = j+rl[k];
                    if( ti>=1&&ti<=n && tj>=1&&tj<=m )
                    {
                        int tid = getid(ti, tj, n, m);
                        a[id][tid] = 1;
                    }
                }
                
                a[id][id] = 2;
                a[id][n*m] = (3-g[i][j])%3;
                id++;
            }
        CL(X,0);
        CL(free_x,0);
        //for (i = 0; i < equ; ++i)
            //for (j = 0; j < var + 1; ++j) scanf("%d",&a[i][j]);
        //Debug();
        
        int free_num = Guass();
        
        if (free_num == -1) printf("无解！\n");
        else if (free_num == -2) printf("无整数解\n");
        else {
//        else if (free_num > 0){
        
            //printf("无穷多解! 自由变元个数为%d\n", free_num);
            //我懂了，我需要确定free_num个数
//            for (i = 0; i < var; ++i){
//                if (free_x[i]) printf("X%d 是不确定的\n",i + 1);
//                else printf("X%d %d\n",i + 1,X[i]);
//            }
            
            
//        }
//        else{
            //我觉得可以！
            int ans = 0;
            for (i = 0 ; i < var; ++i){
                if(X[i]%3 != 0) ans += (X[i]%3+3)%3;
                //printf("X%d %d\n",i + 1,X[i]);
            }
            printf("%d\n",ans);
            for(i=0;i<var;i++)
            {
                int tmp =  (X[i]%3+3)%3;
                for(j=0;j<tmp;j++)
                {
                    int x,y;
                    x = i/m;
                    y = i%m;
                    x++; y++;
                    printf("%d %d\n",x,y);
                }
            }
        }
        //printf("\n");
    }
    return 0;
}
/*
 10
 3 3
 0 0 0
 0 0 0
 0 0 1
 1 2
 0 0
 2 30
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 */

 

2. m个未知数的情况

//
//  main.cpp
//  hdu5755.1
//
//  Created by New_Life on 16/8/4.
//  Copyright © 2016年 chenhuan001. All rights reserved.
//

#include <iostream>
#include <string.h>
#include <stdio.h>
#include <algorithm>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <string>

#define CL(a,num) memset((a),(num),sizeof(a))
#define iabs(x)  ((x) > 0 ? (x) : -(x))
#define Min(a,b) (a) > (b)? (b):(a)
#define Max(a,b) (a) > (b)? (a):(b)

#define ll long long
#define inf 0x7f7f7f7f
#define MOD 100000007
#define lc l,m,rt<<1
#define rc m + 1,r,rt<<1|1
#define pi acos(-1.0)
#define test puts("<------------------->")
#define maxn 100007
#define M 100007
#define N 33
using namespace std;
//freopen("din.txt","r",stdin);

int a[N][N],X[N];//分别记录增广矩阵和解集
int free_x[N];//记录自由变量
int equ,var;//分别表示方程组的个数和变量的个数

int GCD(int x,int y){
    if (y == 0) return x;
    return GCD(y,x%y);
}
int LCM(int x,int y){
    return x/GCD(x,y)*y;
}
void Debug(void)
{
    int i, j;
    for (i = 0; i < equ; i++)
    {
        for (j = 0; j < var + 1; j++)
        {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}
// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，
//但无整数解,-1表无整数解,-1表示无解，0表示唯一解，大于0表示无穷解，
//并返回自由变元的个数)

//ax + by = gcd(a,b)
//传入固定值a,b.放回 d=gcd(a,b), x , y
void extendgcd(ll a,ll b,ll &d,ll &x,ll &y)
{
    if(b==0){d=a;x=1;y=0;return;}
    extendgcd(b,a%b,d,y,x);
    y-=x*(a/b);
}

//Ax=1(mod M)，gcd(A,M)==1
//输入：10^18>=A,M>=1
//输出：返回x的范围是[1,M-1]
ll GetNi(ll A,ll MM)
{
    ll rex=0,rey=0;
    ll td=0;
    extendgcd(A,MM,td,rex,rey);
    return (rex%MM+MM)%MM;
}

int Guass(){
    int i,j,k,col;
    CL(X,0); CL(free_x,1);
    
    for (k = 0,col = 0; k < equ && col < var; k++, col++){
        int max_r = k;
        for (i = k + 1; i < equ; ++i){
            if (iabs(a[i][col]) > iabs(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){
            //可以随意定义的变量
            X[col] = 0;//强制赋值为0
            free_x[col] = 0;
            k--;
            //cout<<k<<endl;
            continue;
        }
        for (i = k + 1; i < equ; ++i){
            if (a[i][col] != 0){
                int lcm = LCM(a[k][col],a[i][col]);
                int ta = lcm/iabs(a[i][col]); int tb = lcm/iabs(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])%3+3)%3;
                }
            }
        }
    }
    //Debug();
    // 1. 无解的情况:
    for (i = k; i < equ; ++i){
        if (a[i][col] != 0) return -1;
    }
    // 2. 无穷解的情况:
    if (k < var){
        
        int num = 0,freeidx=0;
        for (i = k - 1; i >= 0; --i){
            num = 0;
            int tmp = a[i][var];
            for (j = 0; j < var; ++j){
                if (a[i][j] != 0 && free_x[j]){
                    num++;
                    freeidx = j;
                }
            }
            if (num > 1) continue;
            tmp = a[i][var];
            for (j = 0; j < var; ++j){
                if (a[i][j] && j != freeidx) tmp -= a[i][j]*X[j];
            }
            //这里也要
            
            int k2 = (tmp%3+3)%3;
            int k1 = (a[i][freeidx]%3+3)%3;
            if(k1!=0)
            {
                X[freeidx] = k2*(int)GetNi(k1, 3);
            }
            else
            {
                X[freeidx] = 0;
                //printf("X[%d]为任意？\n",i);
            }
            
            //X[freeidx] = tmp/a[i][freeidx];
            free_x[freeidx] = 0;
        }
        return var - k;
    }
    // 3. 唯一解的情况:
    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];
        }
        //强行搞一发
        int k2 = (tmp%3+3)%3;
        int k1 = (a[i][i]%3+3)%3;
        if(k1!=0)
        {
            X[i] = k2*(int)GetNi(k1, 3);
        }
        else
        {
            X[i] = 0;
            //printf("X[%d]为任意？\n",i);
        }
        //X[i] = tmp/a[i][i];//不整除？
    }
    return 0;
}



int g[33][33];
int save[33][33][33];

int getni(int x)
{
    if(x==1) return 2;
    if(x==2) return 1;
    return 0;
}

void func(int x,int y)
{
    g[x][y] = (g[x][y]+2)%3;
    g[x+1][y] = (g[x+1][y]+1)%3;
    g[x][y+1] = (g[x][y+1]+1)%3;
    g[x-1][y] = (g[x-1][y]+1)%3;
    g[x][y-1] = (g[x][y-1]+1)%3;
}

int main() {
    int T;
    cin>>T;
    while(T--)
    {
        int n,m;
        scanf("%d%d",&n,&m);
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
                scanf("%d",&g[i][j]);
        memset(save,0,sizeof(save));
        
        for(int i=1;i<=m;i++)
        {
            save[1][i][i] = 1;
        }
        for(int i=2;i<=n;i++)
        {
            for(int j=1;j<=m;j++)
            {
                for(int k=0;k<=m;k++)
                {
                    int tmp = (save[i-1][j][k]*2+save[i-1][j+1][k]+save[i-1][j-1][k]) + save[i-2][j][k];
                    if(k == 0) tmp += g[i-1][j];
                    tmp %= 3;
                    save[i][j][k] = getni(tmp);
                }
            }
        }
        
        //然后构建矩阵
        memset(a,0,sizeof(a));
        equ = m;
        var = m;
        for(int j=1;j<=m;j++)
        {
            for(int k=0;k<=m;k++)
            {
                int tmp = save[n][j][k]*2+save[n][j+1][k]+save[n][j-1][k] + save[n-1][j][k];
                if(k==0)
                {
                    tmp += g[n][j];
                    a[j-1][m] = getni(tmp%3);
                }
                else a[j-1][k-1] = tmp%3;
            }
        }
        CL(X,0);
        CL(free_x,0);
        int free_num = Guass();
        if(free_num == -1 || free_num == -2){ cout<<free_num<<" 错误"<<endl; continue;}
        int ans = 0;
        int saveans[33][33];
        memset(saveans,0,sizeof(saveans));
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
            {
                int tmp = 0;
                for(int k=0;k<=m;k++)
                {
                    if(k==0) tmp += save[i][j][0];
                    else tmp += save[i][j][k]*X[k-1];
                    tmp %= 3;
                }
                ans += tmp;
                saveans[i][j] = tmp;
            }
        
        printf("%d\n",ans);
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
            {
                for(int k=0;k<saveans[i][j];k++)
                {
                    printf("%d %d\n",i,j);
                    //func(i,j);
                }
            }
        
//        for(int i=1;i<=n;i++)
//        {
//            for(int j=1;j<=m;j++) printf("%d ",g[i][j]);
//            printf("\n");
//        }
    }
    return 0;
}
/*
 10
 1 3
 0 0 1
 3 3
 0 0 0
 0 0 0
 0 0 1
 1 2
 0 0
 2 30
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 */