闲聊高斯消元入门

前言:前一段时间看了看高斯消元,一直想把学习过程写下来,但有点懒=-=现在补上


一、关于高斯消元的题目常见的类型有:

1 . 求唯一整数解

2 . 求唯一浮点数解

3 . 求解得个数,等同于求自由元个数,矩阵的秩。

然后其中一般都会让你判断无解的情况


二、高斯消元的步骤

第一步:将给定矩阵化成阶梯矩阵。

第二步:根据题目的要求 :

1 . 判断是否无解 

2 . 求出解的个数是否唯一

3 . 求出唯一解。

这是整体的思路,细节我就不详细讲了,难度都不大,而且网上都有好文章=-=


三、例题

注:以下代码可能差异很大,因为我也是一边学一边改进的。。


1、POJ 2947 Widget Factory


类型:求唯一整数解,同时判断无解和多解的情况。


代码:

#include<cstdio>
#include<cstring>
#include<algorithm>

using namespace std;

const int mod = 7;
const int MAXN = 355;
/****记得两个数组最好初始化为0****/
int mar[MAXN][MAXN];
int ans[MAXN];

inline int q_pow(int x,int y){
    int res = 1;
    while(y){
        if(y&1){
            res = res * x % mod;
        }y>>=1;
        x = x * x % mod;
    }return res ;
}

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

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

void debug(int n,int m){
    for(int i=0;i<n;i++){
        for(int j=0;j<=m;j++){
            printf("%d%c",mar[i][j],j!=m?' ':'\n');
        }
    }printf("\n");
}
/****** n个方程,m个变量******/
int gauss(int n,int m){
    //debug(n,m);
/*****1.转换为阶梯式*******/
    for(int col=0,row=0;col<=m&&row<n;col++,row++){
        int mx = row;
        for(int i=row+1;i<n;i++){
            if(mar[mx][col]<mar[i][col]){
                mx = i;
            }
        }
        if(mar[mx][col]==0){
            row--;
            continue;
        }
        for(int i=col;i<=m;i++){
            swap(mar[mx][i],mar[row][i]);
        }

        //debug(n,m);
        for(int i=row+1;i<n;i++){
            if(mar[i][col]){//debug(n,m);
                int d = lcm(mar[row][col],mar[i][col]);
                //printf("d=%d\n",d);
                int a = d/mar[row][col],b = d/mar[i][col];

                for(int j=col;j<=m;j++){//printf("j=%d\n",j);
                    mar[i][j] =( ( mar[row][j]*a - mar[i][j]*b ) % mod + mod ) % mod;
                }
            }
        }
    }
    //debug(n,m);
    //printf("in\n");
/******检查是否无解********/
int k = 0;/***k为无效的方程****/
for(int i=0;i<n;i++){
    for(int j=0;j<=m;j++){
        if(j!=m&&mar[i][j]){
            break;
        }else if(j==m){
            if(i<n&&mar[i][m]){
                return 0 ;
            }else if(!mar[i][m]){
                k++;
            }
        }
    }
}

if(n-k<m){
    return 2;
}


/******从最后一层阶梯向上依次求出未知量的值******/

for(int i=m-1;i>=0;i--){
    int val = mar[i][m];/***解为浮点数修改处***/
    for(int j=m-1;j>i;j--){
        val = ( ( val - ans[j] * mar[i][j] ) % mod +mod ) % mod;
    }ans[i] = (val * q_pow(mar[i][i],mod-2) % mod + mod)%mod;
}

return 1;

}

inline int get(char*s){
    switch(s[0]){
        case 'M':return 1;
        case 'T':return s[1]=='U'?2:4;
        case 'W':return 3;
        case 'F':return 5;
        case 'S':return s[1]=='A'?6:7;
    }
}

int main()
{
    int n,m;
    while(scanf("%d%d",&n,&m)==2&&n+m){
        memset(mar,0,sizeof(mar));
        memset(ans,0,sizeof(ans));

        char s1[5],s2[5];
        for(int i=0,t,a;i<m;i++){
            scanf("%d%s%s",&t,s1,s2);
            int t1 = get(s1),t2 = get(s2);
            mar[i][n]=(t2-t1+1+mod)%mod;
            while(t--){
                scanf("%d",&a);
                a--;
                mar[i][a]++;
                mar[i][a]%=mod;
            }
        }int k = gauss(m,n);
        if(k==1){
            for(int i=0;i<n;i++){
                printf("%d%c",ans[i]<3?ans[i]+7:ans[i],i!=n-1?' ':'\n');
            }
        }else if(k==0){
            printf("Inconsistent data.\n");
        }else {
            printf("Multiple solutions.\n");
        }
    }return 0;
}


2、HDU 2262 Where is the canteen


类型:求唯一浮点数解

代码:

#include<queue>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;

const int MAXN = 300;
double mat[MAXN][MAXN];
double ans[MAXN];

int gauss(int n,int m){
    for(int col=0,row=0;col<=m&&row<n;col++,row++){
        int mx = row;
        for(int i=row+1;i<n;i++){
            if(fabs(mat[mx][col])<fabs(mat[i][col])){
                mx = i;
            }
        }
        if(mat[mx][col]==0){
            row--;
            continue;
        }int sign = mat[mx][col]<0?-1:1;
        for(int i=col;i<=m;i++){
            swap(mat[mx][i],mat[row][i]);
            mat[row][i]*=sign;
        }

        for(int i=row+1;i<n;i++){
            if(mat[i][col] != 0){
                double t = mat[i][col]/mat[row][col];
                for(int j=col;j<=m;j++){
                    mat[i][j] -= mat[col][j]*t  ;
                }
            }
        }
    }

for(int i=m-1;i>=0;i--){
    double val = 1.0 * mat[i][m];
    for(int j=m-1;j>i;j--){
        val =  ( val - ans[j] * mat[i][j] );
    }ans[i] = val / mat[i][i] ;
}

return 1;

}


const int X[4]={-1,0,0,1};
const int Y[4]={0,1,-1,0};


int n,m;
char mp[20][20];
int id[20][20];;
bool is[20][20];
bool vis[20][20];
queue< pair<int,int> > Q;
#define m_p(x,y) make_pair(x,y)
#define fi first
#define se second


inline int init(){
    bool res = false ;
    while(!Q.empty()){
        int sz = Q.size();
        while(sz--){
            int x = Q.front().fi;
            int y = Q.front().se;
            if(mp[x][y]=='$'){
                is[x][y]=true;
                res = true;
            }Q.pop();
            for(int i=0,a,b;i<4;i++){
                a = x + X[i];
                b = y + Y[i];
                if(a>0&&b>0&&a<=n&&b<=m&&mp[a][b]!='#'){
                    if(!vis[a][b]){
                        vis[a][b] = true;
                        Q.push(m_p(a,b));
                    }
                }
            }
        }
    }int all = 0;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m;j++){
            if(vis[i][j]&&!is[i][j]){
                id[i][j]=all++;
            }
        }
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m;j++){
            if(vis[i][j]&&!is[i][j]){
                int e = id[i][j];
                double k = 0.;
                for(int l=0;l<4;l++){
                    int a = i+X[l];
                    int b = j+Y[l];
                    if(a>0&&b>0&&a<=n&&b<=m&&vis[a][b]){
                        if(is[a][b]){
                            k+=1.0;
                            continue;
                        }
                        int& t = id[a][b];
                        mat[e][t]=-1.0;
                        k++;
                    }
                }
                mat[e][e] = k;
                mat[e][all]= k;
            }
        }
    }
    if(!res){
        return false ;
    }
    gauss(all,all);
    return true;
}


int main()
{
    while(scanf("%d%d",&n,&m)==2){
        memset(vis,0,sizeof(vis));
        memset(is,0,sizeof(is));
        memset(mat,0,sizeof(mat));
        memset(ans,0,sizeof(ans));
        while(!Q.empty())Q.pop();
        int x0,y0;
        for(int i=1;i<=n;i++){
            scanf("%s",mp[i]+1);
            for(int j=1;j<=m;j++){
                if(mp[i][j]=='@'){
                    x0 = i; y0 = j;
                    vis[i][j]=true;
                    Q.push(m_p(i,j));
                }
            }
        }
        if(init()){
            printf("%.6f\n",ans[id[x0][y0]]);
        }else {
            printf("-1\n");
        }
    }return 0;
}


三、HDU 3364 Lanterns

类型:求解的个数

代码:

/***
注意:
    1.两个数组记得初始化
    2.n个方程,m个变量,mar[][0~m-1]存储对应的系数,mar[][m]存储方程的值
    3.关于返回值
        0:表示无解,阶梯式中出现0,0,0~x(x不为0)的情况。
        -1:表示唯一解。
        正数:阶梯式中出现的有效方程数目少于变量数目,表示有多组解,(变量数-有效方程数)个无法确定的变元。
    4.如果有mod,不会有负数的情况,取逆元进行运算。
    5.如果全程浮点数运算,看另一份代码。
***/



#include<cstdio>
#include<cstring>
#include<algorithm>

using namespace std;

const int mod = 2;
const int MAXN = 355;
int mar[MAXN][MAXN];
int ans[MAXN];

inline int q_pow(int x,int y){
    int res = 1;
    while(y){
        if(y&1){
            res = res * x % mod;
        }y>>=1;
        x = x * x % mod;
    }return res ;
}

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

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

void debug(int n,int m){
    printf("------------checking------------\n");
    for(int i=0;i<n;i++){
        for(int j=0;j<=m;j++){
            printf("%d%c",mar[i][j],j!=m?' ':'\n');
        }
    }printf("------------checking------------\n");
}
/****** n个方程,m个变量******/
int gauss(int n,int m){
/*****1.转换为阶梯式*******/
    for(int col=0,row=0;col<=m&&row<n;col++,row++){
        int mx = row;
        for(int i=row+1;i<n;i++){
            if(mar[mx][col]<mar[i][col]){
                mx = i;
            }
        }
        if(mar[mx][col]==0){
            row--;
            continue;
        }
        for(int i=col;i<=m;i++){
            swap(mar[mx][i],mar[row][i]);
        }

        for(int i=row+1;i<n;i++){
            if(mar[i][col]){
                int d = lcm(mar[row][col],mar[i][col]);
                int a = d/mar[row][col],b = d/mar[i][col];
                for(int j=col;j<=m;j++){
                    mar[i][j] =( ( mar[row][j]*a - mar[i][j]*b ) % mod + mod ) % mod;/***mod修改处***/
                }
            }
        }
    }
/******检查是否无解或者多解********/
int k = 0;/***k为无效的方程数目,无效方程为全为0的方程****/
for(int i=0;i<n;i++){
    for(int j=0;j<=m;j++){
        if(j!=m&&mar[i][j]){
            break;
        }else if(j==m){
            if(i<n&&mar[i][m]){
                return 0 ;
            }else if(!mar[i][m]){
                k++;
            }
        }
    }
}
/******若有效的方程数目小于变量数,则有(变量数-有效方程数)个无法确定的变元****/
if(n-k<m){
    return m-n+k;
}


/******从最后一层阶梯向上依次求出未知量的值******/

for(int i=m-1;i>=0;i--){
    int val = mar[i][m];/***解为浮点数修改处***/
    for(int j=m-1;j>i;j--){
        val = ( ( val - ans[j] * mar[i][j] ) % mod +mod ) % mod;/***mod修改处***/
    }ans[i] = (val * q_pow(mar[i][i],mod-2) % mod + mod)%mod;/***mod修改处***/
}

return -1;

}


int A[MAXN][MAXN];

int main()
{
    int T,t=0;scanf("%d",&T);
    while(T--){
        int n,m;
        memset(mar,0,sizeof(mar));
        memset(ans,0,sizeof(ans));
        memset(A,0,sizeof(A));
        scanf("%d%d",&n,&m);
        for(int i=0,t,x;i<m;i++){
            scanf("%d",&t);
            while(t--){
                scanf("%d",&x);
                x--;
                A[x][i]++;
            }
        }
        printf("Case %d:\n",++t);
        int q;scanf("%d",&q);
        while(q--){
            for(int i=0;i<n;i++){
                for(int j=0;j<m;j++){
                    mar[i][j]=A[i][j];
                }
            }memset(ans,0,sizeof(ans));
            for(int i=0;i<n;i++){
                scanf("%d",&mar[i][m]);
            }//printf("begin------\n");
            //debug(n,m);
            int k = gauss(n,m);
            if(!k){
                printf("0\n");
            }else if(k==-1){
                printf("1\n");
            }else {
                printf("%I64d\n",1LL<<k);
            }
        }
    }return 0;
}

四、 HDU5833 Zhu and 772002


类型 : 求解的个数。

代码:

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<map>
#include<set>
using namespace std;
/******
!!!!!!!!
    Remember !!!;
        use init() function
/***** o(n) ******/

namespace prime_table{

const int MAX_N = 2000;
const long long mod = 2;
const long long p = 1e9+7;
int all=0;
int pr[MAX_N/10+1000];
bool isp[MAX_N+10];

inline void init(){
    all = 0;
    memset(isp,0,sizeof(isp));
    for(int i=2;i<=MAX_N;i++){
        if(!isp[i]){
            pr[all++] = i;
        }for(int j=0;j<all;j++){
            long long t = 1LL*pr[j]*i ;
            if(t<=MAX_N){
                isp[t] = true;
                if(i%pr[j]==0)break;
            }else {
                break;
            }
        }
    }

//printf("all=%d\n",all);


return ;
}

}
using namespace prime_table;

/**** ×éºÏ*****/
long long qpow(long long x,long long y){
    long long res =1;
    while(y){
        if(y&1){
            res = res * x % p;
        }x = x * x % p;
        y >>= 1;
    }return res ;
}

long long ni[305];
long long pre[305];

void iniit(){
    ni[0]=0;
    pre[1] = pre[0] = ni[1]=1;
    for(int i=2;i<=300;i++){
        pre[i] = pre[i-1] * i % p;
        ni[i] = qpow(i,p-2);
    }
}




/***
注意:
    1.两个数组记得初始化
    2.n个方程,m个变量,mar[][0~m-1]存储对应的系数,mar[][m]存储方程的值
    3.关于返回值
        0:表示无解,阶梯式中出现0,0,0~x(x不为0)的情况。
        -1:表示唯一解。
        正数:阶梯式中出现的有效方程数目少于变量数目,表示有多组解,(变量数-有效方程数)个无法确定的变元。
    4.如果有mod,一般不会有负数的情况,然后取逆元。
    5.如果全程浮点数运算,看另一份代码。
***/




const int MAXN = 355;
bool vis[MAXN];
long long e[MAXN][MAXN];
long long mar[MAXN][MAXN];
long long ans[MAXN];

inline void oper(int id,long long x){
    id--;
    for(int i=0;i<all&&pr[i]<=x;i++){
        int k = 0;
        while(x%pr[i]==0){
            vis[i] = true;
            e[i][id]=!e[i][id];
            x/=pr[i];
        }
    }
    return ;
}

inline int work(int n){
    int cnt = 0;
    for(int i=0;i<303;i++){
        if(vis[i]){
            for(int j=0;j<n;j++){
                mar[cnt][j] = e[i][j];
            }mar[cnt++][n]=0;
        }
    }return cnt;
}

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

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

void debug(int n,int m){
    printf("------------checking------------\n");
    for(int i=0;i<n;i++,printf("\n")){
        for(int j=0;j<=m;j++){
            printf("%d%c",mar[i][j],j!=m?' ':'\n');
        }
    }printf("------------checking------------\n");
}
/****** n个方程,m个变量******/
int gauss(int n,int m){//debug(n,m);
/*****1.转换为阶梯式*******/
    for(int col=0,row=0;col<=m&&row<n;col++,row++){
        int mx = row;
        for(int i=row+1;i<n;i++){
            if(mar[mx][col]<mar[i][col]){
                mx = i;
            }
        }
        if(mar[mx][col]==0){
            row--;
            continue;
        }
        for(int i=col;i<=m;i++){
            swap(mar[mx][i],mar[row][i]);
        }

        for(int i=row+1;i<n;i++){
            if(mar[i][col]){
                for(int j=col;j<=m;j++){
                    mar[i][j] =mar[row][j] ^ mar[i][j] ;/***mod修改处***/
                }
            }
        }
    }//debug(n,m);
/******检查是否无解或者多解********/
int k = 0;/***k为无效的方程数目,无效方程为全为0的方程****/
for(int i=0;i<n;i++){
    for(int j=0;j<=m;j++){
        if(j!=m&&mar[i][j]){
            break;
        }else if(j==m){
            if(i<n&&mar[i][m]){
                return -1 ;
            }else if(!mar[i][m]){
                k++;
            }
        }
    }
}
/******若有效的方程数目小于变量数,则有(变量数-有效方程数)个无法确定的变元****/

    return m-n+k;


}


int main(){
    //freopen("data.in","r",stdin);
    //freopen("my_ans.out","w",stdout);
    init();iniit();
    int T,t=0;
    scanf("%d",&T);
    while(T--){
        int n;
        scanf("%d",&n);
        long long x ;
        memset(vis,0,sizeof(vis));
        memset(e,0,sizeof(e));
        for(int i=1;i<=n;i++){
            scanf("%I64d",&x);
            oper(i,x);
        }
        int res = gauss(work(n),n);
        //printf("ans = %d\n",res);
        if(res<=0){
            res = 0;
        }else {
            res = (int)(( qpow(2,1LL*res)-1 + p)%p);
        }
        printf("Case #%d:\n%d\n",++t,res);
    }return 0;
}

/***

15
1345 947 1979 1032 1286 1262 879 510 1671 1291 1552 802 1470 1828 1524








0
***/



评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值