实数,取模,异或意义下的高斯消元法

今日讲线代的时候发现正好是去年这个时候过的高斯消元板子,但是听完之后意识到我去年对这东西一点理解都没有,于是 deep 科研了一下。

写代码的时候出现了多次一遍写对的情况,这可太爽了兄弟。

用矩阵表示方程组

一些疑似废话的引入 我们先从最接近日常的地方说起,要解一个二元一次方程组,你下一步会?

随便选一个方程,再随便选一个未知数,乘一个数让它和另一个方程的系数一样或者相反数,再加或减把这两个方程合起来变成一元一次方程。

把那个变量解出来,带回去,就有了两个变量的解。

你说得对但是我可能无法在数学试卷上正确模拟此过程

这时候我给你 \(n\) 元一次方程组,你可能就没有耐心模拟这个过程了,但是没关系计算机有。

高斯消元法就是模拟了这一个过程,可以得到这 \(n\) 个未知数的解,同时也可以判断无解、无限解的情况。

首先说下怎么处理输入的方程,我们当然不会用一个字符串把方程存起来,发现我们只关心每个方程中各个未知数的系数,当然不存在系数就是 \(0\),我们用一个二维数组存下所有的方程中各个未知数的系数。

具体来说,定义 \(f_{i,j}\) 表示第 \(i\) 个方程中未知数 \(a_j\) 的系数。特殊的,用 \(f_{i,n+1}\) 位置表示等于的常数 \(b\)

接下来想想我们可以给一个等式干什么,用以下三条来更严谨的描述小学学过的 “等式的基本性质”:

  • 交换两个等式,方程组的解不变。
  • 给其中一个等式的两边同时乘上一个非 \(0\) 实数,方程组的解不变。
  • 将其中两个等式相加组成一个新的等式并加入到方程组中,方程组的解不变。

通过这些变换,我们的目标是把原矩阵变成特殊的矩阵,特殊的矩阵有两种,我们称其为上三角矩阵和对角矩阵。

上三角矩阵就是不包括主对角线(左上到右下的对角线称为主对角线)下方的值全是 \(0\) 的矩阵,长这样:

\[\begin{bmatrix} a & b & c & x\\ 0 & d & e & y\\ 0 & 0 & f & z\\ \end{bmatrix} \]

当然 \(a,b,c,x,y\cdots\) 这些也有可能是 \(0\)

最右侧一列是上文提到的常数项。考虑这样的矩阵会对应怎么样的方程。

  • 对于第三行的方程,也就是 \(fx_3=z\),我们可以直接得到 \(x_3\) 的值。
  • 对于第二行的方程,也就是 \(dx_2+ex_3=y\),因为我们已知 \(x_3\),因此也可以得到 \(x_2\) 的值。
  • 经过不断往上代入,我们可以求出所有未知数的值。

还有一种矩阵是对角矩阵,对角矩阵是除了主对角线上外的值全是 \(0\) 的矩阵,长这样:

\[\begin{bmatrix} a & 0 & 0 & x\\ 0 & b & 0 & y\\ 0 & 0 & c & z\\ \end{bmatrix} \]

这种矩阵就更好办了,可以直接得出所有未知数的值。

你可能有疑问,如果 \(a,b,c\) 这些等于 \(0\) 会发生什么。

举个例子,如果 \(a=0\)\(x\ne 0\),此时原方程组是无解的。

如果 \(a=0\)\(x=0\),那么这个变量可以任意取值,我们称这样的变量为 自由元。一个方程组有自由元,意味着这个方程有无穷多个解。

实数意义下的高斯消元

知道了我们的目标,高斯消元就是利用上面三个性质,把原矩阵变成特殊矩阵的过程,显然第二种对角线矩阵是更方便的,我们下文的解说都是变换为对角线矩阵的算法。

设未知数总数为 \(n\),方程总数为 \(m\),过程如下:

首先枚举 \(i\in [1,n]\) (一般是从小到大)表示当前我们要消元的变量 \(a_i\),任取一个满足 \(j\in [1,m],f_{j,i} \ne 0\)\(j\),交换行 \(i,j\)

接下来进行消元,我们的目标就是除了 \(i,i\) 的位置 \(f_{i,i} \ne 0\),对于第 \(i\) 列的其它所有位置,都应该是 \(0\)

有个很简单的想法,就是枚举所有 \(j\in [1,m]\),让第 \(j\) 行的所有数(包括 \(m+1\) 的常数项)都乘上 \(\frac{f_{i,i}}{f_{j,i}}\),再用 \(i\) 行的每一位对应减去 \(j\) 行的每一位。当然,若 \(f_{j,i}\) 本来就是 \(0\),跳过这一行。

乘上 \(\frac{f_{i,i}}{f_{j,i}}\) 会让 \(f_{j,i}\) 变成 \(f_{i,i}\),再对应相减就会让 \(f_{j,i}\) 变成 \(0\),于是就达成了我们的目标。

由于我们事先交换了行,此时一定保证至少从左上角 \(i-1 \times i-1\) 的矩阵一定是一个对角线矩阵,这样子相乘或相减不会在之前消成 \(0\) 的位置造出不为 \(0\) 的数来,只会对对角线上的数造成影响,于是我们可以放心大胆的做操作,不会影响我们之前的结果。

处理完每一个位置,我们就可以统计答案了,也就是枚举 \(i\),去看 \(f_{i,i}\)\(f_{i,m+1}\) 的值,根据我们上文的讨论求解或判断其它情况。

来些模板:

P3389 【模板】高斯消元法

今年重构的代码真的比去年好理解了很多,而且一次就过了。

这题给定了 \(n\) 个方程,任意多个方程要变下写法,当然我们下面也有这样的模板。

code

Show me the code
#define rd read()
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll read(){
	ll x=0,f=1;
	char c=getchar();
	while(c>'9'||c<'0'){if(c=='-') f=-1;c=getchar();}
	while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c^48);c=getchar();}
	return x*f;
}
double x[105];
double f[114][114];
int n;
bool _swap(int b){
  int s=b;
  for(int i=b+1;i<=n;i++){
    if(fabs(f[i][b]-0)>=1e-9){
      for(int j=1;j<=n+1;j++){
        swap(f[s][j],f[i][j]);
      }
      return 0;
    }
  }
  return 1;
}
int main(){

	cin>>n;
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n+1;j++){
			cin>>f[i][j];
		}
	}
	for(int i=1;i<=n;i++){
    int tg=i;
    if(fabs(f[i][tg]-0)<=1e-9){
      if(_swap(i)){
        //换不动了,意味着这一位是 0,由于本题不区分无解和无限解的输出,因此直接输出即可。
        cout<<"No Solution"<<'\n';
			  return 0;
      }
    }
    for(int to=1;to<=n;to++){
      if(i==to)continue;
      double p=f[to][tg]/f[i][tg];
      for(int j=1;j<=n+1;j++){
        f[to][j]-=f[i][j]*p;
      }
    }
  }
	for(int i=1;i<=n;i++){
		printf("%.2lf\n",f[i][n+1]/f[i][i]);
	}

	return 0;
}

异或意义下的高斯消元

我们要求解异或方程组的解,首先把系数放在矩阵里,显然此时矩阵里的数尽可能是 \(0\)\(1\)

接下来依然是枚举 \(i\in [1,n]\) (一般是从小到大)表示当前我们要消元的变量 \(a_i\),任取一个满足 \(j\in [1,m],f_{j,i} \ne 0\)\(j\),交换行 \(i,j\)

枚举所有 \(j\in [1,m],f_{j,i} \ne 0\),让第 \(j\) 行 每一位异或上第 \(i\) 行的每一位,这样就可以让 \(f_{j,i} =0\),可以达成目标了。

经过此过程,可得到一个对角线矩阵,求值即可。

当然若 \(f_{i,i}=0\),也代表一个自由元。

来个例题:

ABC366G XOR Neighbors

模拟题意可得到 \(n\) 个异或方程组,但是能直接消元吗?

你会发现直接让所有点权为 \(0\) 就是一个解,但是题目不让你赋值 \(0\)

也就是说在这题中,如果异或消元消出来 \(f_{i,i} = 1\)\(i\) 行的其它位置都是 \(0\),意味着这个位置必须是 \(0\),也就是说是无解的。

别忘了还有自由元这个大手子,如果 \(f_{i,i} = 0\) 我们可以给自由元 \(x_i\) 赋值 \(2^{i}\),这样不会让两个变量互相影响,由于只有 \(\le 60\) 个点,是够的。

接下来考虑所有 \(f_{i,i}=1\),对于所有 \(k \in [1,n],k\ne i,f_{i,k}=1\),让 \(a_i\) 异或上 \(a_k\),即得 \(a_i\) 的答案。

由于存在自由元,此时矩阵一定不是一个对角线矩阵,这样的 \(k\) 一定是存在的。

这可太合理了。于是这题就做完了。

注意输出答案时在判断下答案是不是 \(0\)

code

Show me the code
#define psb push_back
#define mkp make_pair
#define ls p<<1
#define rs (p<<1)+1
#define rep(i,a,b) for( int i=(a); i<=(b); ++i)
#define per(i,a,b) for( int i=(a); i>=(b); --i)
#define rd read()
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll read(){
  ll x=0,f=1;
  char c=getchar();
  while(c>'9'||c<'0'){if(c=='-') f=-1;c=getchar();}
  while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c^48);c=getchar();}
  return x*f;
}
const int N=100;
int n,m;
vector<int> e[N];
bitset<N> bs[N];
bitset<N> vis;
void dfs(int u,int fa){
	vis[u]=1;
  for(int v:e[u]){
    bs[u][v]=1;
    if(v==fa)continue;
    if(!vis[v]){
    	vis[v]=1;
    	dfs(v,u);
		}
  }
  return ;
}
bool _swap(int b){
  int s=b;
  for(int i=b+1;i<=n;i++){
    if(bs[i][s]){
      for(int j=1;j<=n+1;j++){
        bool b1=bs[b][j],b2=bs[i][j];
        swap(b1,b2);
        bs[b][j]=b1;bs[i][j]=b2;
      }
      return 0;
    }
  }
  return 1;
}
ll ansl[N];
int main(){
  
  cin>>n>>m;
  for(int i=1;i<=m;i++){
    int u,v;cin>>u>>v;
    e[u].push_back(v);
    e[v].push_back(u);
  }
  dfs(1,0);
  // pr();
  for(int i=1;i<=n;i++){
    if(!bs[i][i]){
      if(_swap(i))continue;
      // pr();
    }
    for(int j=1;j<=n;j++){
      if(i==j)continue;
      // pr();
      if(bs[j][i]){
        for(int k=1;k<=n;k++){
          bs[j][k]=bs[j][k]^bs[i][k];
        }
      }
      // pr();
    }
  }
  ll tot=1;
  for(int i=1;i<=n;i++){
    if(bs[i][i]==0){
      bool exs=0;
      for(int j=i+1;j<=n;j++){
        exs|=bs[i][j];
      }
      if(exs){
        cout<<"No";return 0;
      }
      ansl[i]=(1ll<<tot)-1;
      tot++;
    }
    if(bs[i][i]){
      bool exs=0;
      for(int j=i+1;j<=n;j++){
        exs|=bs[i][j];
      }
      if(!exs){
        cout<<"No";return 0;
      }
    }
  }
  for(int i=n;i>=1;i--){
    if(!ansl[i]&&bs[i][i]){
      for(int j=i+1;j<=n;j++){
        ansl[i]=ansl[i]^(ansl[j]*bs[i][j]);
      }
    }
  }
  for(int i=1;i<=n;i++){
    if(ansl[i]==0){
      cout<<"No";return 0;
    }
  }
  cout<<"Yes"<<'\n';
  for(int i=1;i<=n;i++){
    cout<<ansl[i]<<' ';
  }
  
  return 0;
}

取模意义下的高斯消元

令模数为 \(p\)

我们对实数意义下的高斯消元稍作修改,对于选定的第 \(i\) 行和要消元的第 \(j\) 行。我们首先给第 \(i\) 行的所有数都乘上 \(\mathrm{inv}(f_{i,i})\),那么此时 \(f_{i,i}\) 就是 \(1\) 了,接下来把所有的 \(k\in [1,m]\),让 \(f_{j,k} \leftarrow f_{i,k} \times f_{j,k}\),就可以达成消元的目的了。

理解同实数高斯消元,这样子相乘或相减不会在之前消成 \(0\) 的位置造出不为 \(0\) 的数来,只会对对角线上的数造成影响,于是我们可以放心大胆的做操作,不会影响我们之前的结果。

于是就没了,其它操作同实数意义下的高斯消元,但是要注意处处取模。

这里,若模数为合数,可能存在一些数没有乘法逆元,但是题目应该是不会给合数模数情况的。

来个例题,顺便送给你取模意义下的高斯消元的板子:

CF1616F Tricolor Triangles

一个小观察:由于原题中限制了填充的颜色只能是 \(1,2,3\),那么 “三元环上的三条边颜色各不相同或全部相同” 完全等价于边颜色之和为 \(3\) 的倍数。

于是枚举出图中所有三元环,注意到这里点数很小,用邻接矩阵枚举即可。

然后可以构造出最多 \(O(m\sqrt{m})\) 个方程,可以证明。

于是解方程,判断有无解即可。

下面的代码是完全版的对于 \(n\) 个未知数,\(m\) 个方程的方程组,在取模 \(mod\) 意义下的高斯消元模板,这里让 \(mod=3\)

注释较详细。

code

Show me the code
#define psb push_back
#define mkp make_pair
#define ls p<<1
#define rs (p<<1)+1
#define rep(i,a,b) for( int i=(a); i<=(b); ++i)
#define per(i,a,b) for( int i=(a); i>=(b); --i)
#define rd read()
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll read(){
  ll x=0,f=1;
  char c=getchar();
  while(c>'9'||c<'0'){if(c=='-') f=-1;c=getchar();}
  while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c^48);c=getchar();}
  return x*f;
}
const int N=5000;
const int mod=3;
int expc=0;
int f[N][N];
int g[N][N];
int n,m;
ll ans[N];
ll ksm(ll n,ll k){
  ll res=1,f=n;
  while(k){
    if(k&1){
      res=res*f;
    }
    f=f*f;
    k>>=1;
  }
  return res%mod;
}
ll inv(ll a){
  return ksm(a,mod-2);
}
void solve(){
  cin>>n>>m;
  memset(f,0,sizeof f);
  memset(g,0,sizeof g);
  expc=0;
  memset(ans,0,sizeof ans);
  for(int i=1;i<=m;i++){
		int u,v,w;
		cin>>u>>v>>w;
		if (w!=-1) {
			expc++;
			f[expc][i]=1;
			f[expc][m+1]=w%3;
		}
		g[u][v]=g[v][u]=i;
	}
	for(int i=1;i<=n;i++)
		for(int j=i+1;j<=n;j++)
			for(int k=j+1;k<=n;k++)
				if(g[i][j]&&g[j][k]&&g[k][i]){
					expc++;
					f[expc][g[i][j]]=f[expc][g[j][k]]=f[expc][g[k][i]]=1;
				}
  
  int pro=0;
  for(int i=1;i<=m;i++){// 枚举消到了哪个元
    int tg=0;
    for(int j=pro+1;j<=expc;j++){// 任意找到一个不为 0 位置
      if(f[j][i]!=0){
        tg=j;
        break;
      }
    }
    if(!tg){// 如果没找到,为自由元,随便给一个值
      ans[i]=0;continue;
    }
    pro++;
    swap(f[tg],f[pro]);// 交换目标行到当前行
    if(f[pro][i]!=1){// 主元位置归一操作
      int inver=inv(f[pro][i]);// 先把这一行都乘上逆元
      for(int j=i;j<=m+1;j++){
        f[pro][j]=f[pro][j]*inver%mod;
      }
    }
    for(int j=1;j<=expc;j++){
      if(j==pro||!f[j][i])continue;
      int d=f[j][i]%mod;
      for(int k=1;k<=m+1;k++){// 对于主元位置,这样处理后一定是 0
        f[j][k]=(f[j][k]-f[pro][k]*d+mod*mod)%mod;
      }
    }
  }
  // 判断并处理答案即可
  for(int i=pro+1;i<=expc;i++){
    if(f[i][m+1]){
      cout<<-1<<'\n';return ;
    }
  }
  for(int i=1;i<=pro;i++){
    int p=1;
    while(f[i][p]==0)p++;
    ans[p]=f[i][m+1];
  }
  for(int i=1;i<=m;i++){
    if(ans[i]==0)cout<<3<<' ';
    else cout<<ans[i]<<' ';
  }
  cout<<'\n';
  return ;
}
int main(){
  
  int T;cin>>T;
  while(T--){
    solve();
  }
  
  return 0;
}

终于写完了。半个晚自习手敲废了(

posted @ 2025-07-26 21:33  hm2ns  阅读(96)  评论(0)    收藏  举报