作用：
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;
}
![CF1394C 計算幾何[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)