实数,取模,异或意义下的高斯消元法
今日讲线代的时候发现正好是去年这个时候过的高斯消元板子,但是听完之后意识到我去年对这东西一点理解都没有,于是 deep 科研了一下。
写代码的时候出现了多次一遍写对的情况,这可太爽了兄弟。
用矩阵表示方程组
一些疑似废话的引入
我们先从最接近日常的地方说起,要解一个二元一次方程组,你下一步会?随便选一个方程,再随便选一个未知数,乘一个数让它和另一个方程的系数一样或者相反数,再加或减把这两个方程合起来变成一元一次方程。
把那个变量解出来,带回去,就有了两个变量的解。
你说得对但是我可能无法在数学试卷上正确模拟此过程
这时候我给你 \(n\) 元一次方程组,你可能就没有耐心模拟这个过程了,但是没关系计算机有。
高斯消元法就是模拟了这一个过程,可以得到这 \(n\) 个未知数的解,同时也可以判断无解、无限解的情况。
首先说下怎么处理输入的方程,我们当然不会用一个字符串把方程存起来,发现我们只关心每个方程中各个未知数的系数,当然不存在系数就是 \(0\),我们用一个二维数组存下所有的方程中各个未知数的系数。
具体来说,定义 \(f_{i,j}\) 表示第 \(i\) 个方程中未知数 \(a_j\) 的系数。特殊的,用 \(f_{i,n+1}\) 位置表示等于的常数 \(b\)。
接下来想想我们可以给一个等式干什么,用以下三条来更严谨的描述小学学过的 “等式的基本性质”:
- 交换两个等式,方程组的解不变。
- 给其中一个等式的两边同时乘上一个非 \(0\) 实数,方程组的解不变。
- 将其中两个等式相加组成一个新的等式并加入到方程组中,方程组的解不变。
通过这些变换,我们的目标是把原矩阵变成特殊的矩阵,特殊的矩阵有两种,我们称其为上三角矩阵和对角矩阵。
上三角矩阵就是不包括主对角线(左上到右下的对角线称为主对角线)下方的值全是 \(0\) 的矩阵,长这样:
当然 \(a,b,c,x,y\cdots\) 这些也有可能是 \(0\)。
最右侧一列是上文提到的常数项。考虑这样的矩阵会对应怎么样的方程。
- 对于第三行的方程,也就是 \(fx_3=z\),我们可以直接得到 \(x_3\) 的值。
- 对于第二行的方程,也就是 \(dx_2+ex_3=y\),因为我们已知 \(x_3\),因此也可以得到 \(x_2\) 的值。
- 经过不断往上代入,我们可以求出所有未知数的值。
还有一种矩阵是对角矩阵,对角矩阵是除了主对角线上外的值全是 \(0\) 的矩阵,长这样:
这种矩阵就更好办了,可以直接得出所有未知数的值。
你可能有疑问,如果 \(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;
}
终于写完了。半个晚自习手敲废了(

浙公网安备 33010602011771号