高斯消元总结

https://zybuluo.com/ysner/note/1106886

本质

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

列方程

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

实现（加减方程组）

step 1：交换（将当前列系数绝对值最大的放到当前行）
step 2：将每行第一系数化为1
step 3：消成倒三角（在下面的方程中减去上面的方程）
step 4：回代求解

il void Gauss()
{
  fp(i,1,n)//枚举列
    {
      re int now=i;
      fp(j,i+1,n)//枚举行
    if(fabs(a[j][i])+eps>fabs(a[now][i])) now=j;
      if(now^i) swap(a[i],a[now]);
      fp(j,i+1,n)
    fq(k,n+1,i)
    a[j][k]-=a[i][k]*a[j][i]/a[i][i];
    }
  fq(i,n,1)
    {
      ans[i]=a[i][n+1];
      fq(j,n,i+1) ans[i]-=ans[j]*a[i][j];
      ans[i]/=a[i][i];
    }
}

最关键的就是\(a[j][k]-=a[i][k]*a[j][i]/a[i][i]\)了。
此时上面所有方程第一系数在\(/a[i][i]\)后已经被视为\(1\)。
想想平时怎么解方程的。
先把上面方程乘上当前方程第一系数（即\(*a[j][i]\))，再同位项相减即可。

实现（异或方程组)

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

il void Gauss()
{
  fp(i,1,n)
    {
      re int now=i;
      fp(j,i+1,n) if(a[j][i]>a[now][i]) now=j;
      if(now^i) swap(a[now],a[i]);
      if(!a[i][i]) continue;
      fp(j,i+1,n) if(a[j][i]) a[j]^=a[i];
    }
  fq(i,n,1)
    fq(j,i-1,1)
    if(a[j][i]) a[j]^=a[i];
}

特殊情况

无解：第一系数为0,常数非0
无数解：第一系数为0,常数为0
在开头交换后判一下系数最大值是否为0,有一个就记\(flag=1\)，回代时方程形式为\(x=b\),判b是否为0即可。
（注意先判无解，判完后若\(flag=1\),就是无穷解）

例题

T2 hihoCoder1196 高斯消元法二

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

bitset<50>a[50];
il void Gauss()
{
  fp(i,1,30)
    {
      re int now=i;
      fp(j,i+1,30)
	if(a[j][i]>a[now][i]) now=j;
      if(now^i) swap(a[now],a[i]);
      fp(j,i+1,30) if(a[j][i]) a[j]^=a[i];
    }
  fq(i,30,1)
    fq(j,i-1,1)
    if(a[j][i]) a[j]^=a[i];
}
int main()
{
  fp(i,1,5) scanf("%s",s[i]+1);
  fp(i,1,5)
    fp(j,1,6)
    {
      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);
      a[p][p]=1;a[p][31]=(s[i][j]-'0')^1;
      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;
    }
  Gauss();
  return 0;
}

T3Luogu2962 [USACO09NOV]灯Lights

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

显然高斯消元解异或方程组。
而且由于方程内值只有\(0/1\)，可以开\(bitset\)使复杂度降到\(O(\frac{n^3}{64})\)。
注意到这种方程解起来很容易得到形如\(0x=0\)的表达式（由于保证有解，不考虑\(0x=b\))。此时我们要枚举这个\(x\),才能继续往回代。

bitset<50>a[50];
int n,m,ans=1e9,tag[50];
il void Gauss()
{
  fp(i,1,n)
    {
      re int now=i;
      fp(j,i+1,n) if(a[j][i]>a[now][i]) now=j;
      if(now^i) swap(a[now],a[i]);
      if(!a[i][i]) continue;
      fp(j,i+1,n) if(a[j][i]) a[j]^=a[i];
    }
}
il void dfs(re int x,re int tot)
{
  if(tot>=ans) return;
  if(!x) {ans=tot;return;}
  if(a[x][x])
    {
      re int t=a[x][n+1];
      fp(i,x+1,n) if(a[x][i]) t^=tag[i];
      tag[x]=t;
      dfs(x-1,tot+t);
    }
  else
    {
      tag[x]=0;dfs(x-1,tot);
      tag[x]=1;dfs(x-1,tot+1);
    }
}
int main()
{
  n=gi();m=gi();
  fp(i,1,m)
    {
      re int u=gi(),v=gi();
      a[u][v]=a[v][u]=1;
    }
  fp(i,1,n) a[i][i]=a[i][n+1]=1;
  Gauss();
  dfs(n,0);
  printf("%d\n",ans);
  return 0;
}

T4Luogu2447 [SDOI2010]外星千足虫

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

太板子了。
注意至少要\(n\)个方程。

il int Gauss(re int m)
{
  fp(i,1,m) fp(j,1,n) a[i]=aa[i];
  fp(i,1,n)
    {
      re int now=i;
      fp(j,i+1,m) if(a[j][i]>a[now][i]) now=j;
      if(now^i) swap(a[now],a[i]);
      if(a[i][i]==0) return 0;
      fp(j,i+1,m) if(a[j][i]) a[j]^=a[i];
    }
  fq(i,m,1)
    fq(j,i-1,1)
    if(a[j][i]) a[j]^=a[i];
  if(mm==m) fp(i,1,n) ans[i]=a[i][n+1];
  return 1;
}
int main()
{
  n=gi();mm=m=gi();
  fp(i,1,m)
    {
      scanf("%s",s+1);re int len=strlen(s+1);
      fp(j,1,len)
      if(s[j]=='1') a[i][j]=aa[i][j]=1;
      a[i][n+1]=aa[i][n+1]=gi();
    }
  if(!Gauss(m)) {puts("Cannot Determine");return 0;}
  re int l=n,r=m;
  while(l<r)
    {
      re int mid=l+r>>1;
      if(Gauss(mid)) r=mid;
      else l=mid+1;
    }
  printf("%d\n",l);
  fp(i,1,n)
    if(ans[i]&1) puts("?y7M#");
    else puts("Earth");
  return 0;
}

T5Luogu2973 [USACO10HOL]赶小猪

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

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

il void Gauss()
{
  fp(i,1,n)//枚举列数
    {
      re int now=i;
      fp(j,i+1,n)//枚举行
    if(fabs(a[j][i])+eps>fabs(a[now][i])) now=j;
      if(now^i) swap(a[i],a[now]);
      fp(j,i+1,n)
    fq(k,n+1,i)
    a[j][k]-=a[i][k]*a[j][i]/a[i][i];
    }
  fq(i,n,1)
    {
      ans[i]=a[i][n+1];
      fq(j,n,i+1) ans[i]-=ans[j]*a[i][j];
      ans[i]/=a[i][i];
    }
}
int main()
{
  memset(h,-1,sizeof(h));
  n=gi();m=gi();p=gi();q=gi();pi=(double)1.0*p/q;
  fp(i,1,m)
    {
      re int u=gi(),v=gi();
      add(u,v);
    }
  a[1][n+1]=1;
  fp(u,1,n)
    {
      a[u][u]=1;
      for(re int i=h[u];i+1;i=e[i].next)
    {
      re int v=e[i].to;
      a[u][v]=-1.0*(1-pi)/in[v];
    }
    }
  Gauss();
  fp(i,1,n) printf("%.9lf\n",ans[i]*pi);
  return 0;
}

T6 bzoj3270博物馆

给定一个无向联通图，两个人分别在\(a\),\(b\)两点。在每一轮决策中，他们有\(pi\)的概率选择不动，有\(1-pi\)的概率等概率到任意一相邻房间。问在每个房间相遇的概率。
这题挺新奇的是吧。。。
要求每个房间的概率，我们就要求出每个状态的概率（即\(a\)在哪个点，\(b\)在哪个点）。
所以我们把一个状态的概率设为未知数。
方程？
\(P_{a在i，b在j}-\sum P_{a\in i的邻点，b\in i的邻点状态转移至此的概率}=0\)
（~~懒得详写分类讨论~~)
注意不要从已相遇状态转移过来就成。

void Gauss()
{
    fp(i,1,tot)
    {
        re int now=i;
        fp(j,i+1,tot) if(fabs(a[j][i])+eps>fabs(a[now][i])) now=j;
        if(now^i) swap(a[now],a[i]);
        fp(j,i+1,tot)
        fq(k,tot+1,i)
          a[j][k]-=a[i][k]*a[j][i]/a[i][i];
    }
    fq(i,tot,1)
    {
        fq(j,tot,i+1) a[i][tot+1]-=a[j][tot+1]*a[i][j];
        a[i][tot+1]/=a[i][i];
    }
}
int main()
{
  n=gi();m=gi();A=gi();B=gi();
  fp(i,1,m)
  {
      re int u=gi(),v=gi();
      d[u][v]=d[v][u]=1;++in[u];++in[v];
  }
  fp(i,1,n) scanf("%lf",&p[i]),d[i][i]=1;
  fp(i,1,n) fp(j,1,n) id[i][j]=++tot;
  a[id[A][B]][tot+1]=1;
  fp(i,1,n)
    fp(j,1,n)
    {
      a[id[i][j]][id[i][j]]=1;
      fp(ii,1,n)
        fp(jj,1,n)
          if(d[ii][i]&&d[jj][j]&&ii!=jj)//ii!=jj
          {
              re double p1,p2;
              if(i==ii) p1=p[ii];else p1=(1-p[ii])/in[ii];
              if(j==jj) p2=p[jj];else p2=(1-p[jj])/in[jj];
              a[id[i][j]][id[ii][jj]]+=-p1*p2;//+=???
          }
    }
  Gauss();
  fp(i,1,n) printf("%.6lf ",a[id[i][i]][tot+1]);
  return 0;
}

T7 Luogu3164 [CQOI2014]和谐矩阵

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

直接上方程啊。
啥？输出全为\(0\)？把方程里的自由元赋为\(1\)即可.

int main()
{
    n=gi();m=gi();tot=n*m;
    fp(i,1,n)
    fp(j,1,m)
        {
            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);
            a[p][p]=1;
            if (i>1) a[p][u]=1;
            if (j>1) a[p][l]=1;
            if (i<n) a[p][d]=1;
            if (j<m) a[p][r]=1;
        }
    fp(i,1,tot)
    {
        if (!a[i][i])
        {
            fp(j,i+1,tot)
                if (a[j][i]) {swap(a[i],a[j]);break;}
        }
          fp(j,i+1,tot)
            if (a[j][i]) a[j]^=a[i];
    }
    fq(i,tot,1)
    {
        ans[i]=a[i][tot+1];
        fq(j,tot,i+1)
            if (a[i][j]) ans[i]^=ans[j];
        if (!a[i][i]) ans[i]=1;
    }
    for (int i=1;i<=n;++i,puts(""))
        fp(j,1,m)
            printf("%d ",ans[pos(i,j)]);
    return 0;
}

T8 Luogu3265 [JLOI2015]装备购买

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

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

il int Insert(re int x)
{
    fp(i,1,m)
    {
        if(fabs(a[x].s[i])<=eps) continue;
        if(vis[i]) fq(j,m,i) a[x].s[j]-=p[i][j]*a[x].s[i]/p[i][i];
        else {vis[i]=1;fp(j,i,m) p[i][j]=a[x].s[j];return 1;}
    }
    return 0;
}
int main()
{
    n=gi();m=gi();
    fp(i,1,n) fp(j,1,m) a[i].s[j]=gi();
    fp(i,1,n) a[i].c=gi();
    sort(a+1,a+1+n);
    //fp(i,1,n) printf("%d %d %d\n",a[i].s[1],a[i].s[2],a[i].s[3]);
    fp(i,1,n)
    if(Insert(i)) ++tot,ans+=a[i].c;
    printf("%d %d\n",tot,ans);
    return 0;
}

posted @ 2018-04-12 10:08 小蒟蒻ysn 阅读(191) 评论(0) 编辑收藏举报

刷新页面返回顶部

ysner