zoukankan      html  css  js  c++  java
  • 高斯消元

    0. 前言

    写这篇博客的缘由是感觉高斯消元有很多细节需要注意。所以不会讲原理,权当作一篇可供复习的资料吧。

    1. 浮点方程组

    例题

    1. 需要注意输出时对于 (0) 的特判。因为当 (x) 为一个接近 (0) 却为负值的时候,直接输出会变成 -0.00
    2. 但有些数据会专门卡你,所以在找主元时进行随机化也是不错的处理方案。
    #include <cstdio>
    
    #define print(x,y) write(x),putchar(y) 
    
    template <class T> inline T read(const T sample) {
        T x=0; int f=1; char s;
        while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
        while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
        return x*f;
    }
    template <class T> inline void write(const T x) {
        if(x<0) return (void) (putchar('-'),write(-x));
        if(x>9) write(x/10);
        putchar(x%10^48);
    }
    
    #include <cmath>
    #include <iostream>
    using namespace std;
    
    const int maxn=55;
    const double eps=1e-9;
    
    double a[maxn][maxn];
    
    void Gauss(int equ,int var) {
    	int pos,row,col;
    	bool che=0;
    	for(row=col=1;row<=equ and col<var;++row,++col) {
    		pos=row;
    		for(int i=row+1;i<=equ;++i)
    			if(fabs(a[i][col])-fabs(a[pos][col])>eps)
    				pos=i;
    		if(pos^row)
    			for(int i=col;i<=var;++i)
    				swap(a[pos][i],a[row][i]);
    		if(fabs(a[row][col])<eps) {
    			che=1,--row;
    			continue;
    		}
    		for(int i=row+1;i<=equ;++i)
    			if(fabs(a[i][col])>eps) {
    				double tmp=a[i][col]/a[row][col];
    				for(int j=col;j<=var;++j)
    					a[i][j]-=tmp*a[row][j];
    				a[i][col]=0;
    				// 直接赋值为 0,减小精度误差的影响 
    			}
    	}
    	for(int i=row;i<=equ;++i)
    		if(a[i][var]>eps) return (void)(puts("-1"));
    	if(che) return (void)(puts("0"));
    	/*
    		equ-row 就是此时的自由变元个数 
    	*/
    	for(int i=equ;i>=1;--i) {
    		for(int j=i+1;j<=equ;++j)
    			a[i][var]-=a[i][j]*a[j][var];
    		a[i][var]/=a[i][i];
    	}
    	for(int i=1;i<=equ;++i) {
    		if(fabs(a[i][var])<eps)
    			a[i][var]=0;
    		printf("x%d=%.2f
    ",i,a[i][var]);
    	}
    }
    
    int main() {
    	int n=read(9);
    	for(int i=1;i<=n;++i)
    		for(int j=1;j<=n+1;++j)
    			a[i][j]=read(9);
    	Gauss(n,n+1);
    	return 0;
    }
    

    2. 异或线性方程组

    例题

    1. 可以用 ( m bitset) 优化。
    2. 题目可能询问解为 (1) 的最小个数。你可以从后往前枚举自由变元的值,从而依次解出 (x)。但是在不能确定自由变元个数的情况时,最好直接使用 ( ext{meet-in-middle})
    3. 如何理解自由变元?若没有自由变元,高斯消元回代时只有一个未知数。当出现自由变元,你会发现在这一个方程中的 (2) 个未知数都可以被视为自由变元,我们可以任取一个来赋任意值使得未知数个数又变成 (1)
    const int maxn=1005;
    
    int n,m,x[maxn];
    bitset <maxn> a[maxn<<1];
    
    void Gauss() {
    	int tmp,row=1,col=1,ans=0;
    	for(;row<=m && col<n;++row,++col) {
    		tmp=row;
    		while(tmp<m && !a[tmp][col]) ++tmp;
    		if(!a[tmp][col]) return (void)puts("Cannot Determine");
    		ans=Max(ans,tmp);
    		if(tmp^row) swap(a[tmp],a[row]);
    		rep(i,row+1,m) {
    			if(!a[i][col]) continue;
    			a[i]^=a[row];
    		} 
    	}
    	fep(i,n-1,1) {
    		x[i]=a[i][n];
    		rep(j,i+1,n-1) {
    			if(a[i][j]) x[i]^=x[j];
    		}
    	}
    	print(ans,'
    ');
    	rep(i,1,n-1) puts((x[i]&1)?"?y7M#":"Earth");
    }
    
    int main() {
    	n=read(9)+1,m=read(9);
    	rep(i,1,m) {
    		int x;
    		rep(j,1,n-1) scanf("%1d",&x),a[i][j]=x;
    		a[i][n]=read(9);
    	}
    	Gauss();
    	return 0;
    }
    

    3. 整数线性方程组

    1. 防止精度误差,用 ( m lcm) 来消元。
    for(int i=row+1;i<equ;i++)
    	if(a[i][col]!=0){
    		int lcm=LCM(fabs(a[i][col]),fabs(a[row][col]));
    		int ta=lcm/abs(a[i][col]),tb=lcm/abs(a[row][col]);
    		if(a[i][col]*a[row][col]<0) tb=-tb; // 注意正负
    		for(int j=col;j<=var;j++)
    			a[i][j]=a[i][j]*ta-a[row][j]*tb;
    	}
    

    4. 模线性方程组

    1. 对于模数为质数,在除以 (a[i][i]) 时可以用费马小定理求解逆元。
    2. 对于模数非质数,题目就要保证 (a[i][i])( m mod) 互质,然后用 ( m exgcd) 求解。求 (a)(p) 意义下的逆元可以这样:用 ( m exgcd) 解出 (ax+bp=1),移个项就能发现 (x) 即为所求。
    void gauss() {
    	int tmp;
    	for(int i = 0; i <= n; ++ i) {
    		/*
    			用 lcm 来消元,所以主元只需要不为 0 即可
    		*/
    		if(! a[i][i]) {
    			tmp = 0;
    			for(int j = i + 1; j <= n; ++ j) if(a[j][i]) {tmp = j; break;}
    			if(! tmp) continue;
    			for(int j = i; j <= n + 1; ++ j) swap(a[tmp][j], a[i][j]);
    		}
    		for(int j = i + 1; j <= n; ++ j) {
    			tmp = a[j][i];
    			if(! a[j][i]) continue;
    			for(int k = i; k <= n + 1; ++ k) a[j][k] = (0ll + 1ll * a[j][k] * a[i][i] % mod - 1ll * tmp * a[i][k] % mod + mod) % mod;
    		}
    	}
    	for(int i = n; ~i; -- i) {
    		for(int j = i + 1; j <= n; ++ j) a[i][n + 1] = (0ll + a[i][n + 1] - 1ll * ans[j] * a[i][j] % mod + mod) % mod;
    		ans[i] = 1ll * a[i][n + 1] * qkpow(a[i][i], mod - 2) % mod;
    	}
    }
    

    5. 矩阵行列式

    戳这

    int Gauss() {
    	int j,tmp,r=1; bool f=0;
    	for(int i=1;i<=n;++i) {
    		for(j=i;j<=n;++j)
    			if(a[j][i]) break;
    		if(j>n) return 0;
    		if(i^j) swap(a[i],a[j]),f^=1;
    		r=1ll*r*a[i][i]%mod;
    		tmp=inv(a[i][i]);
    		for(j=i;j<=n;++j)
    			a[i][j]=1ll*a[i][j]*tmp%mod;
    		for(j=i+1;j<=n;++j) {
    			a[j][i]=mod-a[j][i];
    			for(int k=i+1;k<=n;++k)
    				a[j][k]=inc(a[j][k],1ll*a[j][i]*a[i][k]%mod);
    			a[j][i]=0;
    		}
    	}
    	return f?mod-r:r;
    }
    

    如果模数不为质数呢?可以模拟辗转相除法:将 (a_{i,i})(a_{j,i}) 对消,直到其中之一为 (0) 即可。所以实际上这是 (mathcal O(n^3log v)),可能稍微会慢一点。

    最后 if(!a[i][i]) return 0; 其实等价于判断是否第 (i) 列全为 (0)

    int Det() {
    	int pos,ret=1,tmp; bool f=0;
    	for(int i=1;i<=n;++i) {
    		for(int j=i+1;j<=n;++j)
    			while(a[j][i]) {
    				int t=a[i][i]/a[j][i];
    				for(int k=i;k<=n;++k) {
    					a[i][k]=(a[i][k]-1ll*a[j][k]*t%mod+mod)%mod;
    					swap(a[i][k],a[j][k]);
    				}
    				f^=1;
    			}
    		if(!a[i][i]) return 0;
    		ret=1ll*ret*a[i][i]%mod;
    	}
    	return f?mod-ret:ret;
    }
    
  • 相关阅读:
    二元查找树转化成排序的双向链表——要求不创建新的节点
    MySQL 通配符学习小结
    HDU 1596 find the safest road (最短路)
    webapp开发调试环境--weinre配置
    全局钩子具体解释
    英尺到米的换算
    apache2.2 虚拟主机配置
    HTTP Digest authentication
    前端project师的修真秘籍(css、javascript和其他)
    CODE:BLOCK中的CreateProcess: No such file or directory
  • 原文地址:https://www.cnblogs.com/AWhiteWall/p/12642332.html
Copyright © 2011-2022 走看看