作用：
1.解n元一次线性方程组
2.可以用来求矩阵的秩
3.求可逆方阵的逆矩阵
其本质就是通过初等行变换，把线性方程组的增广矩阵化为行阶梯型矩阵。
未优化的高斯消元：（顺序消去法）
基本思想：目标就是把系数矩阵的增广矩阵通过初等行变换一步一步的转化为行阶梯型矩阵，然后从后到前把已求出的解代入前面的方程依次求解出各个解。同时注意无解的情况（见代码）。
洛谷P3389【模板】
#include <bits/stdc++.h>
using namespace std;
double G[105][105];
int n;
inline void read(double &x)
{
x=0;
char ch=getchar();//cout<<ch<<endl;
int f=1;
while(ch<'0'||ch>'9')
{
if(ch=='-')
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9')
{
x=x*10+ch-'0';
ch=getchar();
}
x*=f;//cout<<x<<endl;
}
bool Guass()
{//化成上三角矩阵
for(int i=1;i<=n;i++)//遍历系数所在的n列
{
int tmp=i;
for(int j=i+1;j<=n;j++)//找出当前列（i~n行）的最大系数所在的行
{
if(fabs(G[j][i])>fabs(G[tmp][i]))
tmp=j;
}
for(int j=i;j<=n+1;j++)//把当前列最大的系数所在的行放在第i行
swap(G[i][j],G[tmp][j]);
if(G[i][i]==0)//无解
return false;
for(int j=n+1;j>=i;j--)
G[i][j]/=G[i][i];//把当前行的第一个元素化为1
for(int j=i+1;j<=n;j++)
{
for(int k=i+1;k<=n+1;k++)//开始转化上三角
G[j][k]-=G[j][i]*G[i][k];
G[j][i]=0;
}
}
//后向替代算法，逆向求解
for(int i=n;i>=1;i--)//把解放在第n+1列中
{
for(int j=i+1;j<=n;j++)
G[i][n+1]-=G[i][j]*G[j][n+1];
G[i][n+1]/=G[i][i];
G[i][i]=1;
}
return true;
}
int main()
{
while(scanf("%d",&n)!=EOF)
{
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)
read(G[i][j]);
if(Guass())
for(int i=1;i<=n;i++)
printf("%.2lf\n",G[i][n+1]);
else
printf("No Solution\n");
}
return 0;
}
优化的高斯消元：（列主元消去法）
顺序消去法的编程难度较小，但存在弊端。
1.每次运算时，必须保证主对角线上的元素必须非0，否则无法计算（因为有除的操作）。
2.即使主对交线上的元素不为0，如果很小，当它作为除数时，可能导致的误差也是很大的，从而影响算法的稳定性。
列主元消元法
基本思想是：尽可能地选取绝对值大的主元作为除数。
如何判断线性方程组的解的情况：（唯一解，无穷解无解）
1.线性方程组有解的充要条件：
系数阵的秩=增广矩阵的秩
2.线性方程组的解的结构：（唯一解/无穷解）
对于齐次方程组（即常数项=0）：
那么零向量肯定为一组解，因此讨论是否有非零解。
如果r(A)<n,那么它有非零解，而且线性无关的非零解的个数为n-r(A)。
否则只有零解。
对于非齐次方程组：
如果r(A)=r(A,b)=n,那么方程组有唯一解。
如果r(A)=r(A,b)<n,那么方程组有无穷解。
解异或方程组：
通常用在开关问题；
和普通的加减方程组一样，只要把加减操作换成异或操作，高斯消元化简后求出自由变量的个数。（自由变量就是对其取值对结果没有影响）
如果矩阵只有0和1，可以用bitset加速。
相关操作：
poj1830:
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
int n,fnum;
int s[35],e[35];
int equ[35][35];
int fre[35];
void init()
{
fnum=0;
memset(equ,0,sizeof(equ));
memset(fre,0,sizeof(fre));
for(int i=1;i<=n;i++)
equ[i][i]=1;
}
int Guass()
{
int i,v;
for(i=1,v=1;i<=n&&v<=n;i++,v++)//n为方程个数，未知数的个数
{
int pos=i;
for(int j=i+1;j<=n;j++)//找到当前列中最大的系数
if(abs(equ[j][v])>abs(equ[pos][v]))
pos=j;
if(pos!=i)
for(int j=v;j<=n+1;j++)//交换
swap(equ[pos][j],equ[i][j]);
if(equ[i][v]==0)/判断是否为自由变量
{//如果是自由变量，那么处理当前行的下一个变量
i--;
fre[++fnum]=v;//记下自由变量
continue;
}
for(int j=i+1;j<=n;j++)
if(equ[j][v])
for(int k=v;k<=n+1;k++)
equ[j][k]^=equ[i][k];//类似加减法的处理，根据异或的性质
}
for(int j=i;j<=n;j++)//判断无解
if(equ[j][v])//此时
return -1;
if(i<=n)//自由变量的个数
return (n+1-i);
for(int j=n;j>=1;j--)
for(int k=j+1;k<=n;k++)
equ[j][n+1]^=(equ[k][n+1]&&equ[j][k]);
return 0;
}
int main()
{
int K;
scanf("%d",&K);
while(K--)
{//方程：x1^x2^...^xn=s[1]^e[1]
scanf("%d",&n);
init();
for(int i=1;i<=n;i++)
scanf("%d",&s[i]);
for(int i=1;i<=n;i++)
scanf("%d",&e[i]);
for(int i=1;i<=n;i++)
equ[i][n+1]=(s[i]^e[i]);
int i,j;
while(scanf("%d%d",&i,&j),i||j)
equ[j][i]=1;
int t=Guass();
if(t==-1)//无解
printf("Oh,it's impossible~!!\n");
else
printf("%d\n",(1<<t));
}
return 0;
}
(?)hdu4827
bitset加速高斯消元
详解
#include <cstdio>//如果有多组解，求其中一组，对于自由变量可以自己赋值
#include <cstring>
#include <bitset>
#include <cmath>
using namespace std;
const int N=1010;
bitset<N> A[N];
int n,m,fnum;
int fre[N],ans[N];
void Set()
{
for(int i=1;i<=n;i++)
A[i].reset();
int a,b;
for(int i=1;i<=m;i++)
{
scanf("%d%d",&a,&b);
A[a].set(b);
A[b].set(a);
}
for(int i=1;i<=n;i++)
{
if(A[i].count()&1)
{
A[i].set(i);
A[i].set(n+1);
}
}
}
int Guass()
{
int a,b;
fnum=0;
for(a=1,b=1;a<=n&&b<=n;a++,b++)
{
int pos=a;
for(int i=a+1;i<=n;i++)
{
if(A[i][b]>A[pos][b])
pos=i;
}
if(pos!=a)
swap(A[pos],A[a]);
if(A[a][b]==0)
{
a--;
fre[++fnum]=b;
continue;
}
for(int i=a+1;i<=n;i++)
{
if(A[i][b])
A[i]^=A[a];
}
}
int Rank=a-1;//printf("Rank=%d\n",Rank);for(int i=1;i<=fnum;i++)printf("*%d\n",fre[i]);
for(int i=Rank;i>=1;i--)
{
if(A[i][n+1])
{
int j=1;
while(!A[i][j])
j++;
ans[j]=1;
for (int k=i-1;k>=1;k--)
{
if (A[k][j])
{
A[k][j]=0;
A[k].flip(n+1);
}
}
}
}
}
int main()
{
int t,cas=0;
scanf("%d",&t);
while(t--)
{
scanf("%d%d",&n,&m);
for(int i=1;i<=n;i++)
ans[i]=0;
printf("Case #%d:\n",++cas);
Set();
Guass();
for(int i=1;i<=n;i++)
printf("%d",ans[i]);
puts("");
}
return 0;
}
poj2965
关键是如何转化为异或的线性方程组
因为每个开关的状态取决于自己和影响他的开关的操作，因此可以对每一个开关的状态建立有个方程，就可以得到一个16个方程，16个未知数的方程组。（状态：1为开，0,为关；操作：1改变，0不改变）。
其中每一个方程的系数，能影响他的开关系数为1，否则为0，自身为1。常数项取值为最后要求的开关状态异或当前的状态。
然后进行高斯消元，把自由变量全部认为不操作即0，解出非自由变量即可。
用bitset加速后用时32ms，不用用时125ms。
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <bitset>
using namespace std;
char pic[5][5];
bitset<20> A[20];
int ans[20];
int cnt=0;
void Set(int x,int y)
{//建立方程组
int d=(x-1)*4+y;
for(int i=1;i<=4;i++)
{
int t=(x-1)*4+i;
A[d].set(t);
}
for(int i=1;i<=4;i++)
{
int t=(i-1)*4+y;
if(i!=x)
A[d].set(t);
}
}
void init()
{
for(int i=1;i<=16;i++)
A[i].reset();
memset(ans,0,sizeof(ans));
for(int i=1;i<=4;i++)
{
for(int j=1;j<=4;j++)
{
int dx=4*(i-1)+j;
if(pic[i][j]=='+')
A[dx].set(17);
Set(i,j);
}
}
}
void Guass()
{//高斯消元
int a,b;
cnt=0;
for(a=1,b=1;a<=16&&b<=16;a++,b++)
{
int pos=a;
for(int i=a+1;i<=16;i++)
{
if(A[i][b]>A[pos][b])
pos=i;
}
if(pos!=a)
swap(A[pos],A[a]);
if(A[a][b]==0)
{
a--;
continue;
}
for(int i=a+1;i<=16;i++)
if(A[i][b])
A[i]^=A[a];
}
for(int i=a-1;i>=1;i--)
{//求解
if(!A[i][17])
continue;
int j=1;
while(!A[i][j])
j++;
ans[j]=1;
for(int k=i-1;k>=1;k--)
{
if(A[k][j])
{
A[k][j]=0;
A[k].flip(17);
}
}
}
cnt=0;
for(int i=1;i<=16;i++)
if(ans[i])
cnt++;
printf("%d\n",cnt);
for(int i=1;i<=4;i++)
{
for(int j=1;j<=4;j++)
{
if(ans[(i-1)*4+j])
printf("%d %d\n",i,j);
}
}
}
int main()
{
while(scanf("%s",pic[1]+1)!=EOF)
{
for(int i=2;i<=4;i++)
scanf("%s",pic[i]+1);
init();//建立方程组
Guass();
}
return 0;
}
poj3185
开关问题，因为是一个序列，所有要么从左边开始，要么从右边开始，只要进行两次模拟，然后比较大小即可。
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
int blow[22];
int main()
{
while(scanf("%d",&blow[1])!=EOF)
{
for(int i=2;i<=20;i++)
scanf("%d",&blow[i]);
int tmp[22];
for(int i=1;i<=20;i++)
tmp[i]=blow[i];
int a=0,b=0;
for(int i=1;i<=20;i++)
{
if(tmp[i])
{
if(i==20)
a=1000;
a++;
tmp[i+1]=!tmp[i+1];
tmp[i+2]=!tmp[i+2];
}
}
for(int i=20;i>=1;i--)
{
if(blow[i])
{
if(i==1)
b=1000;
b++;
blow[i-1]=!blow[i-1];
blow[i-2]=!blow[i-2];
}
}
printf("%d\n",min(a,b));
}
return 0;
}
hdu4592(变态好题)
看网上的题解都写的很麻烦，没有找到好的思路，记下来，以后再看。
hdu3359
浮点数高斯消元模板题
注意输出要求，不要输出空格！
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
int w,h,d,num;
double grid[12][12],A[110][110];
void Set()
{
num=w*h;
memset(A,0,sizeof(A));
for(int i=1;i<=h;i++)
{
for(int j=1;j<=w;j++)
{
int t=(i-1)*w+j;
int cnt=0;
for(int a=max(i-d,1);a<=i+d&&a<=h;a++)
{
for(int b=max(j-d,1);b<=j+d&&b<=w;b++)
{
if(abs(i-a)+abs(j-b)<=d)
{
int tt=(a-1)*w+b;
A[t][tt]=1.0;
cnt++;
}
}
}//printf("cnt=%d\n",cnt);
A[t][num+1]=grid[i][j]*cnt;//printf("%lf\n",A[t][num+1]);
}
}
}
void Gauss()
{
int a,b;
for(a=1,b=1;a<=num&&b<=num;a++,b++)
{
int pos=a;
for(int i=a+1;i<=num;i++)
{
if(fabs(A[i][b])>fabs(A[pos][b]))
pos=i;
}
if(pos!=a)
{
for(int i=b;i<=num+1;i++)
swap(A[pos][i],A[a][i]);
}
if(A[a][b]==0)
{
a--;
continue;
}
for(int i=a+1;i<=num;i++)
{
double tmp=A[i][b]/A[a][b];
for(int j=b;j<=num+1;j++)
A[i][j]-=tmp*A[a][j];
}
}
for(int i=num;i>=1;i--)
{
for(int j=i+1;j<=num;j++)
A[i][num+1]-=(A[j][num+1]*A[i][j]);
A[i][num+1]/=A[i][i];
}
}
int main()
{
bool flag=0;
//fopen("input.txt","r",stdin);
//freopen("output.txt","w",stdout);
while(scanf("%d%d%d",&w,&h,&d),w||h||d)
{
if(!flag)
flag=1;
else
puts("");
for(int i=1;i<=h;i++)
{
for(int j=1;j<=w;j++)
scanf("%lf",&grid[i][j]);//printf("%lf\n",grid[i][j]);
}
Set();
Gauss();
for(int i=1;i<=num;i++)
{
if(i%w==0)
printf("%8.2lf\n",A[i][num+1]);
else
printf("%8.2lf",A[i][num+1]);
}
}
return 0;
}
hdu3976
给一个电路图，要求求节点1和节点n之间的等效电阻。
思路是看了博客才知道的，要用到电路基础中的基尔霍夫电流定律。流入和流出每个点的电流之和为0。因此可以以每个点的电势为未知数，根据每个节点列出方程组，即可用高斯消元解决。
用列主元消去法，一直会出现除0的情况，后来改为一般的高斯消元才不会。
为了方便计算，假设流入的电流为1，且我是以1节点的电势为高电势，n为低电势，那么两者的电势之差就是等效电阻。规定流出节点为正方向，那么对于1节点而言，1A的电流流出，所以常数项为-1，（移项），4节点也同理。
#include <bits/stdc++.h>
using namespace std;
int n,m;
double A[55][55];
void Set()
{
for(int i=1;i<=n;i++)
{
if(i==1)
A[i][n+1]=-1;
if(i==n)
A[i][n+1]=1;
}
}
void Gauss()
{
int a,b;
for(a=1,b=1;a<=n&&b<=n;a++,b++)
{
int pos=a;
for(int i=a+1;i<=n;i++)
{
if(fabs(A[i][b])>fabs(A[pos][b]))
pos=i;
}
if(a!=pos)
{
for(int i=b;i<=n+1;i++)
swap(A[a][i],A[pos][i]);
}
if(A[a][b]==0)
{//printf("****************%d %d\n",a,b);
a--;
continue;
}
for(int j=b+1;j<=n+1;j++)
A[a][j]/=A[a][b];
A[a][b]=1;
for(int i=a+1;i<=n;i++)
{
if(A[i][b])
{
for(int j=b+1;j<=n+1;j++)
A[i][j]-=A[a][j]*A[i][b];
A[i][b]=0;
}
}
//printf("###%lf\n",A[4][4]);
}
for(int i=n;i>=1;i--)
{//printf("%lf\n",A[i][i]);
for(int j=i+1;j<=n;j++)
A[i][n+1]-=A[i][j]*A[j][n+1];
}
printf("%.2lf\n",A[n][n+1]-A[1][n+1]);
}
int main()
{
int T,cas=0;
scanf("%d",&T);
while(T--)
{
scanf("%d%d",&n,&m);
memset(A,0,sizeof(A));
int a,b;
double c;
for(int i=1;i<=m;i++)
{
scanf("%d%d%lf",&a,&b,&c);
A[a][b]+=(-1.0/c);
A[b][a]+=(-1.0/c);
A[a][a]+=(1.0/c);
A[b][b]+=(1.0/c);
}
printf("Case #%d: ",++cas);
Set();
Gauss();
}
return 0;
}
京公网安备 11010502036488号