[关闭]
@ysner 2018-04-12T10:06:04.000000Z 字数 5973 阅读 2048

高斯消元总结

高斯消元


本质

模拟加减消元,先消成上三角再代入求解。
据说是小学内容???

列方程

有这个量的影响就赋系数。

实现(加减方程组)

  1. il void Gauss()
  2. {
  3. fp(i,1,n)//枚举列
  4. {
  5. re int now=i;
  6. fp(j,i+1,n)//枚举行
  7. if(fabs(a[j][i])+eps>fabs(a[now][i])) now=j;
  8. if(now^i) swap(a[i],a[now]);
  9. fp(j,i+1,n)
  10. fq(k,n+1,i)
  11. a[j][k]-=a[i][k]*a[j][i]/a[i][i];
  12. }
  13. fq(i,n,1)
  14. {
  15. ans[i]=a[i][n+1];
  16. fq(j,n,i+1) ans[i]-=ans[j]*a[i][j];
  17. ans[i]/=a[i][i];
  18. }
  19. }

最关键的就是了。
此时上面所有方程第一系数在后已经被视为
想想平时怎么解方程的。
先把上面方程乘上当前方程第一系数(即),再同位项相减即可。

Update:只要有无解情况就需要交换。

实现(异或方程组)

解这种方程不难,和解加减方程差不多,原理是不断消去倒三角左边的系数。一开始消下面的系数,往回代时消上面的。
除了continue,代码思路一致。

  1. il void Gauss()
  2. {
  3. fp(i,1,n)
  4. {
  5. re int now=i;
  6. fp(j,i+1,n) if(a[j][i]>a[now][i]) now=j;
  7. if(now^i) swap(a[now],a[i]);
  8. if(!a[i][i]) continue;
  9. fp(j,i+1,n) if(a[j][i]) a[j]^=a[i];
  10. }
  11. fq(i,n,1)
  12. fq(j,i-1,1)
  13. if(a[j][i]) a[j]^=a[i];
  14. }

特殊情况

例题

T1 Luogu3389【模板】高斯消元法

见上

T2 hihoCoder1196 高斯消元法二

据说异或方程组高斯消元和普通高斯消元一模一样,就是把减改成异或。
若只有01建议用bitset压常数

  1. bitset<50>a[50];
  2. il void Gauss()
  3. {
  4. fp(i,1,30)
  5. {
  6. re int now=i;
  7. fp(j,i+1,30)
  8. if(a[j][i]>a[now][i]) now=j;
  9. if(now^i) swap(a[now],a[i]);
  10. fp(j,i+1,30) if(a[j][i]) a[j]^=a[i];
  11. }
  12. fq(i,30,1)
  13. fq(j,i-1,1)
  14. if(a[j][i]) a[j]^=a[i];
  15. }
  16. int main()
  17. {
  18. fp(i,1,5) scanf("%s",s[i]+1);
  19. fp(i,1,5)
  20. fp(j,1,6)
  21. {
  22. re int p=pos(i,j),l=pos(i,j-1),r=pos(i,j+1),u=pos(i-1,j),d=pos(i+1,j);
  23. a[p][p]=1;a[p][31]=(s[i][j]-'0')^1;
  24. if(i>1) a[p][u]=1;if(i<5) a[p][d]=1;if(j>1) a[p][l]=1;if(j<6) a[p][r]=1;
  25. }
  26. Gauss();
  27. return 0;
  28. }

T3Luogu2962 [USACO09NOV]灯Lights

给定一个无向联通图,开始点值全为0。你可以选择将一个点和与其相邻的点同时取反,问最少进行多少次这样的操作可以使点值全为1。

显然高斯消元解异或方程组。
而且由于方程内值只有,可以开使复杂度降到
注意到这种方程解起来很容易得到形如的表达式(由于保证有解,不考虑)。此时我们要枚举这个,才能继续往回代。

  1. bitset<50>a[50];
  2. int n,m,ans=1e9,tag[50];
  3. il void Gauss()
  4. {
  5. fp(i,1,n)
  6. {
  7. re int now=i;
  8. fp(j,i+1,n) if(a[j][i]>a[now][i]) now=j;
  9. if(now^i) swap(a[now],a[i]);
  10. if(!a[i][i]) continue;
  11. fp(j,i+1,n) if(a[j][i]) a[j]^=a[i];
  12. }
  13. }
  14. il void dfs(re int x,re int tot)
  15. {
  16. if(tot>=ans) return;
  17. if(!x) {ans=tot;return;}
  18. if(a[x][x])
  19. {
  20. re int t=a[x][n+1];
  21. fp(i,x+1,n) if(a[x][i]) t^=tag[i];
  22. tag[x]=t;
  23. dfs(x-1,tot+t);
  24. }
  25. else
  26. {
  27. tag[x]=0;dfs(x-1,tot);
  28. tag[x]=1;dfs(x-1,tot+1);
  29. }
  30. }
  31. int main()
  32. {
  33. n=gi();m=gi();
  34. fp(i,1,m)
  35. {
  36. re int u=gi(),v=gi();
  37. a[u][v]=a[v][u]=1;
  38. }
  39. fp(i,1,n) a[i][i]=a[i][n+1]=1;
  40. Gauss();
  41. dfs(n,0);
  42. printf("%d\n",ans);
  43. return 0;
  44. }

T4Luogu2447 [SDOI2010]外星千足虫

给m个方程,最多n个未知数,解出方程并询问这最少需要要前几个方程组。

太板子了。
注意至少要个方程。

  1. il int Gauss(re int m)
  2. {
  3. fp(i,1,m) fp(j,1,n) a[i]=aa[i];
  4. fp(i,1,n)
  5. {
  6. re int now=i;
  7. fp(j,i+1,m) if(a[j][i]>a[now][i]) now=j;
  8. if(now^i) swap(a[now],a[i]);
  9. if(a[i][i]==0) return 0;
  10. fp(j,i+1,m) if(a[j][i]) a[j]^=a[i];
  11. }
  12. fq(i,m,1)
  13. fq(j,i-1,1)
  14. if(a[j][i]) a[j]^=a[i];
  15. if(mm==m) fp(i,1,n) ans[i]=a[i][n+1];
  16. return 1;
  17. }
  18. int main()
  19. {
  20. n=gi();mm=m=gi();
  21. fp(i,1,m)
  22. {
  23. scanf("%s",s+1);re int len=strlen(s+1);
  24. fp(j,1,len)
  25. if(s[j]=='1') a[i][j]=aa[i][j]=1;
  26. a[i][n+1]=aa[i][n+1]=gi();
  27. }
  28. if(!Gauss(m)) {puts("Cannot Determine");return 0;}
  29. re int l=n,r=m;
  30. while(l<r)
  31. {
  32. re int mid=l+r>>1;
  33. if(Gauss(mid)) r=mid;
  34. else l=mid+1;
  35. }
  36. printf("%d\n",l);
  37. fp(i,1,n)
  38. if(ans[i]&1) puts("?y7M#");
  39. else puts("Earth");
  40. return 0;
  41. }

T5Luogu2973 [USACO10HOL]赶小猪

一个无向图,节点1有一个炸弹,在每个单位时间内,有p/q的概率在这个节点炸掉,有1-p/q的概率随机选择一条出去的路到其他的节点上。问最终炸弹在每个节点上爆炸的概率。

看到题就想起了[HNOI2013]游走
好像思路方程都是一样的咦。
不过高斯消元解的应该是炸弹到某点的概率,炸弹爆炸概率最后再乘。

  1. il void Gauss()
  2. {
  3. fp(i,1,n)//枚举列数
  4. {
  5. re int now=i;
  6. fp(j,i+1,n)//枚举行
  7. if(fabs(a[j][i])+eps>fabs(a[now][i])) now=j;
  8. if(now^i) swap(a[i],a[now]);
  9. fp(j,i+1,n)
  10. fq(k,n+1,i)
  11. a[j][k]-=a[i][k]*a[j][i]/a[i][i];
  12. }
  13. fq(i,n,1)
  14. {
  15. ans[i]=a[i][n+1];
  16. fq(j,n,i+1) ans[i]-=ans[j]*a[i][j];
  17. ans[i]/=a[i][i];
  18. }
  19. }
  20. int main()
  21. {
  22. memset(h,-1,sizeof(h));
  23. n=gi();m=gi();p=gi();q=gi();pi=(double)1.0*p/q;
  24. fp(i,1,m)
  25. {
  26. re int u=gi(),v=gi();
  27. add(u,v);
  28. }
  29. a[1][n+1]=1;
  30. fp(u,1,n)
  31. {
  32. a[u][u]=1;
  33. for(re int i=h[u];i+1;i=e[i].next)
  34. {
  35. re int v=e[i].to;
  36. a[u][v]=-1.0*(1-pi)/in[v];
  37. }
  38. }
  39. Gauss();
  40. fp(i,1,n) printf("%.9lf\n",ans[i]*pi);
  41. return 0;
  42. }

T6 bzoj3270博物馆

给定一个无向联通图,两个人分别在,两点。在每一轮决策中,他们有的概率选择不动,有的概率等概率到任意一相邻房间。问在每个房间相遇的概率。
这题挺新奇的是吧。。。
要求每个房间的概率,我们就要求出每个状态的概率(即在哪个点,在哪个点)。
所以我们把一个状态的概率设为未知数。
方程?

懒得详写分类讨论)
注意不要从已相遇状态转移过来就成。

  1. void Gauss()
  2. {
  3. fp(i,1,tot)
  4. {
  5. re int now=i;
  6. fp(j,i+1,tot) if(fabs(a[j][i])+eps>fabs(a[now][i])) now=j;
  7. if(now^i) swap(a[now],a[i]);
  8. fp(j,i+1,tot)
  9. fq(k,tot+1,i)
  10. a[j][k]-=a[i][k]*a[j][i]/a[i][i];
  11. }
  12. fq(i,tot,1)
  13. {
  14. fq(j,tot,i+1) a[i][tot+1]-=a[j][tot+1]*a[i][j];
  15. a[i][tot+1]/=a[i][i];
  16. }
  17. }
  18. int main()
  19. {
  20. n=gi();m=gi();A=gi();B=gi();
  21. fp(i,1,m)
  22. {
  23. re int u=gi(),v=gi();
  24. d[u][v]=d[v][u]=1;++in[u];++in[v];
  25. }
  26. fp(i,1,n) scanf("%lf",&p[i]),d[i][i]=1;
  27. fp(i,1,n) fp(j,1,n) id[i][j]=++tot;
  28. a[id[A][B]][tot+1]=1;
  29. fp(i,1,n)
  30. fp(j,1,n)
  31. {
  32. a[id[i][j]][id[i][j]]=1;
  33. fp(ii,1,n)
  34. fp(jj,1,n)
  35. if(d[ii][i]&&d[jj][j]&&ii!=jj)//ii!=jj
  36. {
  37. re double p1,p2;
  38. if(i==ii) p1=p[ii];else p1=(1-p[ii])/in[ii];
  39. if(j==jj) p2=p[jj];else p2=(1-p[jj])/in[jj];
  40. a[id[i][j]][id[ii][jj]]+=-p1*p2;//+=???
  41. }
  42. }
  43. Gauss();
  44. fp(i,1,n) printf("%.6lf ",a[id[i][i]][tot+1]);
  45. return 0;
  46. }

T7 Luogu3164 [CQOI2014]和谐矩阵

   我们称一个由0和1组成的矩阵是和谐的,当且仅当每个元素都有偶数个相邻的1。一个元素相邻的元素包括它本身,及他上下左右的4个元素(如果存在)。给定矩阵的行数和列数,请计算并输出一个和谐的矩阵。注意:所有元素为0的矩阵是不允许的。

直接上方程啊。
啥?输出全为?把方程里的自由元赋为即可.

  1. int main()
  2. {
  3. n=gi();m=gi();tot=n*m;
  4. fp(i,1,n)
  5. fp(j,1,m)
  6. {
  7. re int p=pos(i,j),l=pos(i,j-1),r=pos(i,j+1),u=pos(i-1,j),d=pos(i+1,j);
  8. a[p][p]=1;
  9. if (i>1) a[p][u]=1;
  10. if (j>1) a[p][l]=1;
  11. if (i<n) a[p][d]=1;
  12. if (j<m) a[p][r]=1;
  13. }
  14. fp(i,1,tot)
  15. {
  16. if (!a[i][i])
  17. {
  18. fp(j,i+1,tot)
  19. if (a[j][i]) {swap(a[i],a[j]);break;}
  20. }
  21. fp(j,i+1,tot)
  22. if (a[j][i]) a[j]^=a[i];
  23. }
  24. fq(i,tot,1)
  25. {
  26. ans[i]=a[i][tot+1];
  27. fq(j,tot,i+1)
  28. if (a[i][j]) ans[i]^=ans[j];
  29. if (!a[i][i]) ans[i]=1;
  30. }
  31. for (int i=1;i<=n;++i,puts(""))
  32. fp(j,1,m)
  33. printf("%d ",ans[pos(i,j)]);
  34. return 0;
  35. }

T8 Luogu3265 [JLOI2015]装备购买

个装备,每个装备个属性,每个装备还有个价格。如果手里有的装备的每一项属性为它们分配系数(实数)后可以相加得到某件装备,则不必要买这件装备。求最多装备下的最小花费。

看到这道题可以想起线性基的性质:其中没有一个元素能被其它元素异或和求得。
那这题也类似?可以拿线性基的的方法套?
线性基里存维,拿一个表达式进去时,若该表达式对应维系数为0就换维,如果该维没访问过就把这个式子整个塞进去,否则就把这个式子当前维消掉(高斯消元)。
于是,式子只有两个结果,一个是被塞进去,另一个是被消成蛋。
(竟然卡精度)

  1. il int Insert(re int x)
  2. {
  3. fp(i,1,m)
  4. {
  5. if(fabs(a[x].s[i])<=eps) continue;
  6. if(vis[i]) fq(j,m,i) a[x].s[j]-=p[i][j]*a[x].s[i]/p[i][i];
  7. else {vis[i]=1;fp(j,i,m) p[i][j]=a[x].s[j];return 1;}
  8. }
  9. return 0;
  10. }
  11. int main()
  12. {
  13. n=gi();m=gi();
  14. fp(i,1,n) fp(j,1,m) a[i].s[j]=gi();
  15. fp(i,1,n) a[i].c=gi();
  16. sort(a+1,a+1+n);
  17. //fp(i,1,n) printf("%d %d %d\n",a[i].s[1],a[i].s[2],a[i].s[3]);
  18. fp(i,1,n)
  19. if(Insert(i)) ++tot,ans+=a[i].c;
  20. printf("%d %d\n",tot,ans);
  21. return 0;
  22. }
添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注