首页 > 其他 > 详细

学习高斯消元

时间：2015-03-03 15:13:22      阅读：295      评论：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年09月23日 (328)
2021年09月24日 (313)
2021年09月17日 (191)
2021年09月15日 (369)
2021年09月16日 (411)
2021年09月13日 (439)
2021年09月11日 (398)
2021年09月12日 (393)
2021年09月10日 (160)
2021年09月08日 (222)
最新文章
更多>
教程昨日排行
更多>
友情链接
汇智网   PHP教程   插件网  
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣，分享程序人生！