码迷，mamicode.com
首页 > 其他好文 > 详细

学习高斯消元

时间：2015-03-03 15:15:51      阅读：220      评论：0      收藏：0      [点我收藏+]

标签：高斯消元

点击打开链接     习题  行列式

高斯消元问题类型:

用LCM 保持整型

1. 基本的高斯消元,裸模板 HDU3359

2. 开关问题,用^操作代替 -, 求x[i]时候一样用*  poj 1222 1830 1753

3. 枚举自由变元, return -1 是因为出现[0,0,0,0,a]这种情况,return 0 是唯一解,否则是有自由变元

4. 取模方程 (a1*x1+a2*x2...)%P=b在初等行变换后，每次进行消元的时候将所有值模P，最后求解回带的时候利用扩展欧几里得来对每一个ai求一个最小的可行解       

    3*a4 % P = t4，可以表示成3*a4 + K*P = t4  ,d=exgcd(a[i][i],mod,x,y); ans=ans/d*x;    poj 2947

5. 线性基问题,异或运算中求可组成的第K大数.每个数就是equ   二进制位数为var ,消元成线性基. 比如n个数字,m个线性基.组成2^(n-m)个数字.这个K的二进制形式每个位

    的权重刚好是线性基的顺序.如果有 0行,则0为第一小. 点击打开链接      SGU 275

6. 求电阻.依据KCL定理.每个节点电流流入流出相等.电路图是双向图,要从两个方面考虑.依据电流列等式,X[i]是电压.再设定一个0电压处.S流入为1,T流入为-1.

    电流流出为正,流入为负(相反也行),x[u][v]++,x[v][u]++,x[u][u]--,x[v][v]--; 表示u这个点流入流出和值为0电流是双向的,所以双向. HDU 5006 HDU 3976 点击打开链接 点击打开链接

    

求解线性方程的结果会出现三种情况：无解，多解和唯一解。看下图

        技术分享

当某一行出现（0,0，……0，a) 时，方程无解。因为x1*0+x2*0+……+xn*0=a（a不为0）等式无解

当从某一行开始至最后一行，出现（0,0，……0）说明方程存在自由元，则会出现多解

当矩阵的秩等于等式个数时，有唯一解

poj 1753
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define inf 0x3f3f3f3f
#define eps 1e-8
const int MAXN=17;
int N,M,T;
int a[MAXN][MAXN],x[MAXN];
int free_x[MAXN],x_num;
int equ,var;
int Gauss(int (*a)[MAXN]){
    int row,col;
    x_num=0;
    for(row=0,col=0;row<equ && col<var;row++,col++){
        int mxr = row;
        for(int i=row+1;i<equ;i++)
            if( fabs(a[i][col]) - fabs(a[mxr][col]) >eps )
                mxr=i;
        if(mxr != row)
            for(int i=row;i<var+1;i++) swap(a[row][i],a[mxr][i]);///
        if(fabs(a[row][col])< eps){
            free_x[x_num++]=col;row--;continue;
        }
        for(int i=row+1;i<equ;i++) if(fabs(a[i][col])>eps){
            for(int j=col;j<var+1;j++)
                a[i][j] ^= a[row][j];
        }
    }

    for(int i=row;i<equ;i++)
        if(a[i][var]>0) return -1;
    if(row<var) return var-row;
    for(int i=var-1;i>=0;i--) if(a[i][i]){
        x[i]=a[i][var];
        for(int j=i+1;j<var;j++) x[i]^=(a[i][j] * x[j]);
    }
    return 0;
}
char G[4][4];
int b[MAXN][MAXN];
int solve(int (*a)[MAXN]){
    int t=Gauss(a);
    int ans=0;
    if(t==-1) return inf;
    if(t==0) {
        for(int i=0;i<var;i++) if(x[i]) ans++;
        return ans;
    }
    ///枚举自由变元
    ans=inf;
    for(int i=0;i<(1<<t);i++){
        int cnt=0;
        for(int j=0;j<t;j++)
            if(i&(1<<j))
                x[free_x[j]]=1,cnt++;
            else x[free_x[j]]=0;
        for(int j=var-t-1;j>=0;j-- ){
            int idx;
            for(idx=j;idx<var;idx++) if(a[j][idx])
                break;
            x[idx]=a[j][var];
            for(int l=idx+1;l<var;l++)
                x[idx] ^=a[j][l] * x[l];
            cnt+=x[idx];
        }
        ans=min(ans,cnt);
    }
    return ans;
}
int main()
{
    equ=16,var=16;
    for(int i=0;i<4;i++)
       gets(G[i]);
    int ans=inf;
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    for(int i=0;i<16;i++){
        a[i][i]=b[i][i]=1;
        if(G[i/4][i%4]=='b') a[i][16]=0,b[i][16]=1;
        else a[i][16]=1,b[i][16]=0;


    }
    for(int i = 0;i < 4;i++)
        for(int j = 0;j < 4;j++){
            int t = i*4+j;
            if(i > 0)b[(i-1)*4+j][t]=a[(i-1)*4+j][t] = 1;
            if(i < 3)b[(i+1)*4+j][t]=a[(i+1)*4+j][t] = 1;
            if(j > 0)b[i*4+j-1][t]=a[i*4+j-1][t] = 1;
            if(j < 3)b[i*4+j+1][t]=a[i*4+j+1][t] = 1;
        }
//        for(int i=0;i<16;i++)printf("%d ",a[i][16]);cout<<endl;
//        for(int i=0;i<16;i++)printf("%d ",b[i][16]);cout<<endl;
    ans=solve(a);
    ans=min(ans,solve(b));
    if(ans==inf) printf("Impossible\n");
    else printf("%d\n",ans);

    return 0;
}


poj 2947 

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define inf 0x3f3f3f3f
#define eps 1e-8
const int MAXN=400;
int N,M,T;
int a[MAXN][MAXN],x[MAXN];
int free_x[MAXN],x_num;
int equ,var;
inline int gcd(int a,int b){
    return b==0?a:gcd(b,a%b);
}
inline int lcm(int a,int b){
    return a/gcd(a,b)*b;
}
int expgcd(int a,int b,int &x,int &y){
    int q,tmp;
    if(b==0){
        q=a;x=1;y=0;
        return q;
    }
    q=expgcd(b,a%b,x,y);
    tmp=x;x=y;
    y=tmp-(a/b)*y;
    return q;
}


int Gauss(int (*a)[MAXN]){
    int row,col;
    x_num=0;
    for(int i=0;i<=var;i++)
        free_x[i]=1;
    for(row=0,col=0;row<equ && col<var;row++,col++){
        int mxr = row;
        for(int i=row+1;i<equ;i++)
            if( fabs(a[i][col]) - fabs(a[mxr][col]) >eps )
                mxr=i;
        if(mxr != row)
            for(int i=row;i<var+1;i++) swap(a[row][i],a[mxr][i]);///
        if(fabs(a[row][col])< eps){
            free_x[x_num++]=col;row--;continue;
        }
        for(int i=row+1;i<equ;i++) if(fabs(a[i][col])>eps){
            int LCM=lcm(abs(a[i][col]),abs(a[row][col]));
            int ta=LCM/abs(a[i][col]);
            int tb=LCM/abs(a[row][col]);
            ///if(a[i][col]*a[row][col]<0) tb=-tb;
            for(int j=col;j<var+1;j++)
                a[i][j] = ((a[i][j]*ta-a[row][j]*tb)%7+7)%7;
        }///用LCM的方法是因为要保持整型
    }

    for(int i=row;i<equ;i++)
        if(a[i][var]>0) return -1;
    if(row<var){///这样标记自由变元的方法避免了列移动
        for(int k=row-1;k>=0;k--){
            x_num=0;
            int idx;
            for(int j=0;j<var;j++)
                if(a[k][j]>eps && free_x[j]) x_num++,idx=j;
            if(x_num>1) continue;
            int tmp=a[k][var];
            for(int j=0;j<var;j++) if(a[k][j]>eps && j!=idx)
                    tmp =((tmp - a[k][j]*x[j]%7)+7)%7;
            x[idx]=(tmp/a[k][idx])%7;
            free_x[idx]=0;
        }
        return var-row;
    }
    for(int i=var-1;i>=0;i--){
        int tmp=a[i][var];
        for(int j=i+1;j<var;j++)
            tmp =((tmp - a[i][j]*x[j]%7)+7)%7;
        while(tmp%a[i][i]) tmp+=7;
        x[i]=(tmp/a[i][i])%7;
        if(x[i]<3) x[i]+=7;
    }
    return 0;
}

int solve(int (*a)[MAXN]) {
    int t=Gauss(a);
    int ans=0;
    if(t==-1) return inf;
    if(t==0) {
        for(int i=0;i<var;i++) if(x[i]) ans++;
        return ans;
    }
    ///枚举自由变元
    ans=inf;
    for(int i=0;i<(1<<t);i++){
        int cnt=0;
        for(int j=0;j<t;j++)
            if(i&(1<<j))
                x[free_x[j]]=1,cnt++;
            else x[free_x[j]]=0;
        for(int j=var-t-1;j>=0;j-- ){
            int idx;
            for(idx=j;idx<var;idx++) if(a[j][idx])
                break;
            x[idx]=a[j][var];
            for(int l=idx+1;l<var;l++)
                x[idx] ^=a[j][l] * x[l];
            cnt+=x[idx];
        }
        ans=min(ans,cnt);
    }
    return ans;
}
int tans(char *s){
    if(s[0]=='M') return 1;
    else if(s[0]=='W') return 3;
    else if(s[0]=='F') return 5;
    else if(s[0]=='T' && s[1]=='U') return 2;
    else if(s[0]=='T' && s[1]=='H') return 4;
    else if(s[0]=='S' && s[1]=='A') return 6;
    else return 7;
}
int main() {
    while(~scanf("%d%d",&N,&M) && N+M){
        memset(a,0,sizeof(a));
        for(int i=0;i<M;i++){
            int n;
            char s1[10],s2[10];
            scanf("%d%s%s",&n,s1,s2);
            a[i][N]=(tans(s2)-tans(s1)+8)%7 ;
            for(int j=0;j<n;j++){
                int t;
                scanf("%d",&t);
                a[i][--t]++;
                a[i][t]%=7;
            }
        }
        equ=M,var=N;
        int ans=Gauss(a);
        if(ans==-1) puts("Inconsistent data.");
        else if(ans>0) puts("Multiple solutions.");
        else {
            for(int i=0;i<N;i++)
                if(i==N-1) printf("%d\n",x[i]);
                else printf("%d ",x[i]);
        }
    }

    return 0;
}

HDU3976

#include<stdio.h>
#include<string.h>
#include<math.h>
#include<algorithm>
const double eps = 1e-8 ;
using namespace std;

int T ,cas ,N,M;
double mat[55][55] ;
double ans[55] ;

void gauss(){
	int row, col ;
	row = col = 0 ;
	for( ;row<N && col<N;row++,col++){
		int max_r = row ;
		for(int i=row+1;i<N;i++){
			if(fabs(mat[max_r][col]) < fabs(mat[i][col]) )	max_r = i ;
		}
		if(max_r != row){
			for(int j=col;j<=N;j++){
				swap(mat[row][j] ,mat[max_r][j] );
			}
		}
		for(int i=row+1;i<N;i++){
			if(fabs(mat[i][col]) < eps)	continue ;
			double a = - mat[i][col] / mat[row][col] ;
			for(int j=col;j<=N;j++){
				mat[i][j] += mat[row][j]*a ;
			}
		}
	}
	memset(ans , 0 ,sizeof(ans)) ;
	for(int i=row-2;i>=0;i--){
		double res = mat[i][N] ;
		for(int j=i+1;j<N;j++){
			res -= mat[i][j]*ans[j];
		}
		ans[i] = res / mat[i][i] ;
	}
	printf("%.2f\n",ans[0]-ans[N-1]);
}
int main(){
	scanf("%d",&T);
	cas = 0 ;
	while(T--){
		++cas  ;
		scanf("%d %d",&N,&M);
		memset(mat , 0 ,sizeof(mat));
		for(int i=0;i<M;i++){
			int a,b, c; 
			scanf("%d %d %d",&a,&b,&c);
			a -- ; b-- ;
			mat[a][a] -= 1.0/(c*1.0) ;
			mat[b][b] -= 1.0/(c*1.0) ;
			mat[a][b] += 1.0/(c*1.0) ;
			mat[b][a] += 1.0/(c*1.0) ; 
		}
		mat[0][N] = -1 ;
		mat[N-1][N] = 1 ;
		printf("Case #%d: ",cas);
		gauss() ;
	}	
	return 0; 
}




学习高斯消元

标签：高斯消元

原文地址：http://blog.csdn.net/gg_gogoing/article/details/44037631
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！