【已整理】高斯消元法 总

使用高斯消元法的注意点，写在前面：

1.一定不要忘记初始化系数矩阵

2.注意零解是否符合要求

先放出以前的两篇博客链接

http://www.cnblogs.com/milesgo/articles/7219072.html

http://www.cnblogs.com/milesgo/articles/7219269.html

先讨论一下基础知识。

①齐次线性方程组

齐次线性方程组一定有零解。

有两种情况：（1）只有零解（2）有无穷多解，即除了零解，还有无穷多非零解

行数n表示的是方程个数，列数m表示的是变量个数。设矩阵的秩为r。高斯消元后，不全为0的行数即为矩阵的秩。如果r==m，则说明只有零解。如果r<m，则说明有m-r个自由元，对这m-r个自由元任意取值，可以得到无穷多个解。

如果系数矩阵是方阵，则只需要判断行列式是否等于0。若不等于0，说明满秩，则只有零解，若等于0，说明不满秩，则有无穷多解。

齐次线性方程组的解空间对于加法和数乘（线性运算）是封闭的。

②非齐次线性方程组

非齐次线性方程组的解的个数有三种情况：（1）无解（2）有唯一解（3）有无穷多解（如果题目对变量取值有限制，则不是有限多解）

将b向量加到系数矩阵的最后一列，得到增广矩阵。对增广矩阵进行高斯消元。结束后观察增广矩阵，若原系数矩阵的秩不等于增广矩阵的秩，则说明无解（即存在一行，原系数全部化为了零，但是b向量的对应元素不为0）。

如果两者秩相等，设为r，若r==m，则有唯一解，若r<m，则有无穷多解。

接下来总结一下模板，系数矩阵一般都是从0开始放，即0~n-1行，0~m-1列，非齐次方程组的b向量存在第m列。

高斯消元法的时间复杂度是n三次方的

（1）整数取模高斯消元模板 求出唯一解（只适用于有唯一解的情况，其它情况需要稍做更改）

const int MAXN=35;
const int mod=2;

int gcd(int a,int b){
    return b?gcd(b,a%b):a;
}

int lcm(int a,int b){
    return a/gcd(a,b)*b;
}

void Gauss(int a[][MAXN],int n,int m,int &r,int &c){//高斯消元
    r=c=0;//
    for(;r<n&&c<m;r++,c++){
        int maxi=r;
        for(int i=r+1;i<n;i++)
            if (abs(a[i][c])>abs(a[maxi][c]))
                maxi=i;
        if(maxi!=r){
            for(int i = r; i<m + 1; i++)
                swap(a[r][i], a[maxi][i]);
        }//每次都选当前列系数最大的行放到r行的位置
        if(a[r][c]==0){
            r--;//如果当前列系数最大的值也是0，则r在下一轮循环值不变
            continue;
        }
        for(int i=r+1;i<n;i++){//高斯消元的过程
            if(a[i][c]!=0){
                int x=abs(a[i][c]);
                int y=abs(a[r][c]);
                int LCM=lcm(x,y);
                int tx=LCM/x;
                int ty=LCM/y;
                if (a[i][c]*a[r][c]<0)
                    ty=-ty;
                for(int j=c;j<m+1;j++)
                    a[i][j]=((a[i][j]%mod*tx%mod-a[r][j]%mod*ty%mod)%mod+mod)%mod;
            }
        }
    }
}

int Rewind(int a[][MAXN],int x[],int r,int c){//回代求解
    for (int i=r-1;i>=0;i--){
        int t=a[i][c]%mod;
        for (int j=i+1;j<c;j++){
            if (a[i][j]!=0)
                t-=a[i][j]%mod*x[j]%mod;
        }
        x[i]=t/a[i][i]%mod;
        x[i]=(x[i]+mod)%mod;
    }
    return 0;
}

int a[MAXN][MAXN];
int x[MAXN];

使用方法

int main(){
    1.memset(a,0,sizeof(a)); 初始化系数矩阵
    2.系数矩阵赋值
    3.int r,c;
    4.Gauss(a,n,m,r,c);
    5.Rewind(a,x,r,c); x是解向量
    6.输出x
}

 能够直接套用该模板的题：

1.poj1222

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<queue>
#include<cmath>

using namespace std;

const int MAXN=35;
const int mod=2;

int gcd(int a,int b){
    return b?gcd(b,a%b):a;
}

int lcm(int a,int b){
    return a/gcd(a,b)*b;
}

void Gauss(int a[][MAXN],int n,int m,int &r,int &c){//高斯消元
    r=c=0;//
    for(;r<n&&c<m;r++,c++){
        int maxi=r;
        for(int i=r+1;i<n;i++)
            if (abs(a[i][c])>abs(a[maxi][c]))
                maxi=i;
        if(maxi!=r){
            for(int i = r; i<m + 1; i++)
                swap(a[r][i], a[maxi][i]);
        }//每次都选当前列系数最大的行放到r行的位置
        if(a[r][c]==0){
            r--;//如果当前列系数最大的值也是0，则r在下一轮循环值不变
            continue;
        }
        for(int i=r+1;i<n;i++){//高斯消元的过程
            if(a[i][c]!=0){
                int x=abs(a[i][c]);
                int y=abs(a[r][c]);
                int LCM=lcm(x,y);
                int tx=LCM/x;
                int ty=LCM/y;
                if (a[i][c]*a[r][c]<0)
                    ty=-ty;
                for(int j=c;j<m+1;j++)
                    a[i][j]=((a[i][j]%mod*tx%mod-a[r][j]%mod*ty%mod)%mod+mod)%mod;
            }
        }
    }
}

int Rewind(int a[][MAXN],int x[],int r,int c){//回代求解
    for (int i=r-1;i>=0;i--){
        int t=a[i][c]%mod;
        for (int j=i+1;j<c;j++){
            if (a[i][j]!=0)
                t-=a[i][j]%mod*x[j]%mod;
        }
        x[i]=t/a[i][i]%mod;
        x[i]=(x[i]+mod)%mod;
    }
    return 0;
}

int a[MAXN][MAXN];
int x[MAXN];

int main(){
    int t,p;
    scanf("%d",&t);
    for(int cas=1;cas<=t;cas++){
        memset(a,0,sizeof(a));
        for(int i=0;i<5;i++){
            for(int j=0;j<6;j++){
                scanf("%d",&p);
                if(i>0)    a[i*6+j][(i-1)*6+j]=1;
                if(i<4)    a[i*6+j][(i+1)*6+j]=1;
                if(j>0)    a[i*6+j][i*6+j-1]=1;
                if(j<5)    a[i*6+j][i*6+j+1]=1;
                a[i*6+j][i*6+j]=1;
                a[i*6+j][30]=p;
            }
        }
        int r,c;
        Gauss(a,30,30,r,c);
        Rewind(a,x,r,c);
        printf("PUZZLE #%d\n",cas);
        for(int i=0;i<5;i++){
            for(int j=0;j<6;j++){
                printf("%d%c",x[i*6+j],j==5?'\n':' ');
            }
        }
    }
    return 0;
}

 （2）整数取模高斯消元模板，判断是否有解，解是否唯一，求解的个数

只需对（1）中模板做略微更改

const int MAXN=35;
const int mod=2;

int gcd(int a,int b){
    return b?gcd(b,a%b):a;
}

int lcm(int a,int b){
    return a/gcd(a,b)*b;
}

bool Gauss(int a[][MAXN],int n,int m,int &r,int &c){//高斯消元
    r=c=0;
    for(;r<n&&c<m;r++,c++){
        int maxi=r;
        for(int i=r+1;i<n;i++)
            if (abs(a[i][c])>abs(a[maxi][c]))
                maxi=i;
        if(maxi!=r){
            for(int i = r; i<m + 1; i++)
                swap(a[r][i], a[maxi][i]);
        }//每次都选当前列系数最大的行放到r行的位置
        if(a[r][c]==0){
            r--;//如果当前列系数最大的值也是0，则r在下一轮循环值不变
            continue;
        }
        for(int i=r+1;i<n;i++){//高斯消元的过程
            if(a[i][c]!=0){
                int x=abs(a[i][c]);
                int y=abs(a[r][c]);
                int LCM=lcm(x,y);
                int tx=LCM/x;
                int ty=LCM/y;
                if (a[i][c]*a[r][c]<0)
                    ty=-ty;
                for(int j=c;j<m+1;j++)
                    a[i][j]=((a[i][j]%mod*tx%mod-a[r][j]%mod*ty%mod)%mod+mod)%mod;
            }
        }
    }
    for(int i=r;i<n;i++){
        if(a[i][m]!=0)
            return false;
    }//r表示矩阵的秩，则n-r为自由元个数，则quick（mod，n-r）即为解的个数
    return true;
}

int a[MAXN][MAXN];
int x[MAXN],y[MAXN];

int quick(int n,int k){
    int res=1;
    while(k){
        if(k&1)
            res*=n;
        n*=n;
        k>>=1;
    }
    return res;
}

将Gauss函数的返回值更改为bool类型，在Gauss函数的最后加上如下语句

for(int i=r;i<n;i++){
    if(a[i][m]!=0)
        return false;
}//r表示矩阵的秩，则n-r为自由元个数，则quick（mod，n-r）即为解的个数
return true;

在主函数中，加上如下语句

if(!Gauss(a,n,n,r,c))
        printf("Oh,it's impossible~!!\n");//无解
else
    printf("%d\n",quick(2,n-r))//解的个数

  能够直接套用该模板的题：

1.poj1830

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<queue>
#include<cmath>

using namespace std;

const int MAXN=35;
const int mod=2;

int gcd(int a,int b){
    return b?gcd(b,a%b):a;
}

int lcm(int a,int b){
    return a/gcd(a,b)*b;
}

bool Gauss(int a[][MAXN],int n,int m,int &r,int &c){//高斯消元
    r=c=0;
    for(;r<n&&c<m;r++,c++){
        int maxi=r;
        for(int i=r+1;i<n;i++)
            if (abs(a[i][c])>abs(a[maxi][c]))
                maxi=i;
        if(maxi!=r){
            for(int i = r; i<m + 1; i++)
                swap(a[r][i], a[maxi][i]);
        }//每次都选当前列系数最大的行放到r行的位置
        if(a[r][c]==0){
            r--;//如果当前列系数最大的值也是0，则r在下一轮循环值不变
            continue;
        }
        for(int i=r+1;i<n;i++){//高斯消元的过程
            if(a[i][c]!=0){
                int x=abs(a[i][c]);
                int y=abs(a[r][c]);
                int LCM=lcm(x,y);
                int tx=LCM/x;
                int ty=LCM/y;
                if (a[i][c]*a[r][c]<0)
                    ty=-ty;
                for(int j=c;j<m+1;j++)
                    a[i][j]=((a[i][j]%mod*tx%mod-a[r][j]%mod*ty%mod)%mod+mod)%mod;
            }
        }
    }
    for(int i=r;i<n;i++){
        if(a[i][m]!=0)
            return false;
    }//r表示矩阵的秩，则n-r为自由元个数，则quick（mod，n-r）即为解的个数
    return true;
}

int a[MAXN][MAXN];
int x[MAXN],y[MAXN];

int quick(int n,int k){
    int res=1;
    while(k){
        if(k&1)
            res*=n;
        n*=n;
        k>>=1;
    }
    return res;
}

int main(){
    int t,n,p,q;
    scanf("%d",&t);
    for(int cas=1;cas<=t;cas++){
        scanf("%d",&n);
        for(int i=0;i<n;i++)
            scanf("%d",x+i);
        for(int i=0;i<n;i++)
            scanf("%d",y+i);
        memset(a,0,sizeof(a));
        while(scanf("%d%d",&p,&q)){
            if(p==0&&q==0)
                break;
            a[q-1][p-1]=1;
        }
        for(int i=0;i<n;i++)
            a[i][i]=1;
        for(int i=0;i<n;i++)
            a[i][n]=x[i]^y[i];
        int r,c;
        if(!Gauss(a,n,n,r,c))
            printf("Oh,it's impossible~!!\n");//无解
        else
            printf("%d\n",quick(2,n-r));//解的个数
    }
    return 0;
}

 2.hdu5833 Zhu and 772002

题意：给你n（<=300）个数，每个数的最大素因子均不会超过2000，每个数的范围<=1e18。每次选出若干个数相乘，问你有多少种选法可以使得乘积为一个完全平方数

完全平方数，素因数分解后，所有素因数的个数均为偶数个。每个数选或者不选，相当于就是设一个变量xi，当xi为1时表示选这个数，当xi为0时表示不选这个数。而选了这个数的贡献ai就是系数，贡献ai取决于当前这个数有多少个prime[i]这个因子，如果是偶数个，那ai就是0，如果是奇数个，那ai就是1。

那本题就转化为了高斯消元法，唯一需要注意的一点是，要想获得最终的正确答案，需要减1，因为零解不符合题目要求。而本题构造出来的线性方程组是齐次线性方程组，一定有零解，所以这个1是一定要减的。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<queue>
#include<cmath>

using namespace std;
typedef long long int ll;
const int mod=2;
const int N=2005;

int prime[N];
bool vis[N];
int tot;

void get_prime(){
    tot=0;
    memset(vis,false,sizeof(vis));
    for(int i=2;i<N;i++){
        if(!vis[i]){
            prime[tot++]=i;
            for(int j=i+i;j<N;j+=i)
                vis[j]=true;
        }
    }
}

int gcd(int a,int b){
    return b?gcd(b,a%b):a;
}

int lcm(int a,int b){
    return a/gcd(a,b)*b;
}

bool Gauss(int a[][305],int n,int m,int &r,int &c){//高斯消元
    r=c=0;
    for(;r<n&&c<m;r++,c++){
        int maxi=r;
        for(int i=r+1;i<n;i++)
            if (abs(a[i][c])>abs(a[maxi][c]))
                maxi=i;
        if(maxi!=r){
            for(int i = r; i<m + 1; i++)
                swap(a[r][i], a[maxi][i]);
        }//每次都选当前列系数最大的行放到r行的位置
        if(a[r][c]==0){
            r--;//如果当前列系数最大的值也是0，则r在下一轮循环值不变
            continue;
        }
        for(int i=r+1;i<n;i++){//高斯消元的过程
            if(a[i][c]!=0){
                int x=abs(a[i][c]);
                int y=abs(a[r][c]);
                int LCM=lcm(x,y);
                int tx=LCM/x;
                int ty=LCM/y;
                if (a[i][c]*a[r][c]<0)
                    ty=-ty;
                for(int j=c;j<m+1;j++)
                    a[i][j]=((a[i][j]%mod*tx%mod-a[r][j]%mod*ty%mod)%mod+mod)%mod;
            }
        }
    }
    for(int i=r;i<n;i++){
        if(a[i][m]!=0)
            return false;
    }//r表示矩阵的秩，则n-r为自由元个数，则quick（mod，n-r）即为解的个数
    return true;
}

int a[2005][305];
const ll M=1e9+7;

ll quick(ll n,ll k){
    ll res=1;
    while(k){
        if(k&1)
            res=(res*n)%M;
        n=(n*n)%M;
        k>>=1;
    }
    return res;
}

int main(){
    int t,n;
    ll p;
    get_prime();
    scanf("%d",&t);
    for(int cas=1;cas<=t;cas++){
        memset(a,0,sizeof(a));
        scanf("%d",&n);
        for(int i=0;i<n;i++){
            scanf("%lld",&p);
            for(int j=0;j<tot&&p>1;j++){
                if(p%prime[j])  continue;
                int cur=0;
                while(p%prime[j]==0){
                    cur^=1;
                    p/=prime[j];
                }
                a[j][i]=cur;
            }
        }
        int r,c;
        if(!Gauss(a,tot,n,r,c))
            printf("Case #%d:\n0\n",cas);
        else
            printf("Case #%d:\n%lld\n",cas,((quick(2,n-r)-1)%M+M)%M);
    }
    return 0;
}