2024-9-10 高斯消元

最新推荐文章于 2024-10-18 22:48:33 发布
最新推荐文章于 2024-10-18 22:48:33 发布
阅读量835 收藏 15
点赞数 25
文章标签： c++ 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/jyy0806/article/details/142260052
版权

说在前面

  • 一直听说高斯消元这个算法,但是由于自己对数学的抵触情绪,一直没有敢于来学习。直到看见了这道题，和完这道题后的罪恶感让我下定了学习高斯消元的决心，学完感觉并没有想象中那样难，而且这种人机思维的转换十分有趣，值得一学

高斯消元

这是什么

  • 高斯消元，一种在约 O ( n 3 ) O(n^3) O(n3)的时间内求解n元一次方程组的算法，通常作为程序的一部分，可以与期望概率、位运算、dp等配合使用，可以用于求解受后继状态影响的dp方程

引入

  • 考虑我们解方程的步骤：（以二元一次方程为例）

    • 我们知道，普通人解决这个问题通常运用两种方法：加减消元和代入消元

    • 对于一个二元一次方程组：

    { x + 2 y = 3    ( 1 ) 4 x + 4 y = 5    ( 2 ) \begin{cases}x+2y=3 \ \ (1) \\4x+4y=5\ \ (2) \end{cases} {x+2y=3  (1)4x+4y=5  (2)​

    • ( 2 ) − 2 × ( 1 ) (2)-2\times(1) (2)−2×(1)， x = − 1 2 x=-\frac{1}{2} x=−21​，将x带入y中，解得 y = 7 4 y=\frac{7}{4} y=47​

  • 像这样，一个多元一次方程组就可以在进行多次加减消元和带入消元后被求解，高斯消元利用的就是这种思路

实现（以P3389 【模板】高斯消元法 为例）

  • 考虑把这种加减和代入的思路转换为计算机可以处理的方式

前置知识：矩阵的初等行变换

  • 用非零的数a乘以某一整行

  • 把一行的a倍加到另一行

  • 互换任意两行的位置

  • 利用这种变换可以实现加减消元
  1. 将方程组转化为矩阵

  • 我们发现，在解方程过程中真正关注的其实是方程系数，所以把方程的系数单独拿出来，形成一个矩阵，即系数矩阵A，满足 A i , j = 第 i 个方程中的第 j 个未知数的系数 A_{i,j}=第i个方程中的第j个未知数的系数 Ai,j​=第i个方程中的第j个未知数的系数（在每个方程中未知数顺序一样）

  • 再把右侧的常数拿出来，形成一个矩阵，称为常数矩阵B

  • 则整个方程组可表示为增广矩阵A|B，即矩阵A和B拼起来的矩阵：

[ 1 2 3 4 4 5 ] \left[\begin{array}{cc|c} 1&2&3\\ 4&4&5 \end{array}\right] [14​24​35​]
2. 对增广矩阵进行求解

  • 为了方便，我们让第i个方程保留第i个未知数，这样消元完成后的矩阵就是一个只有右上三角有数（普通高斯消元）或只有对角线有数的矩阵（高斯-约旦消元）

  • 普通高斯消元：

  • 枚举未知数i（1->n）进行消元（将矩阵通过矩阵初等变换变换为三角矩阵）

    主元：主元(pivot element)，一种变元。指在消去过程中起主导作用的元素。
    主元就是在矩阵消去过程中,每列的要保留的非零元素,用它可以把该列其他消去
    在阶梯型矩阵中,主元就是每个非零行第一个非零元素就是主元

    • 首先，在整个矩阵中寻找第i个未知量的系数不为零的方程，然后将它们进行互换（如果找不到不为0的方程，则第i个位置数的值可以取任意值，该方程组有多解）

    • 在这一步中，为了减少误差，通常要找系数最大的一个作为主元

    • 接着枚举i之后的方程，把它们中的第i个未知数消掉

  • 消元完成后，对矩阵进行回带操作

    • 倒序枚举行，通过带入消元将后面仍然未被消掉的未知数消掉
inline void gauss(){
    for(int i=1;i<=n;i++){
        int maxx=i;
        for(int j=i+1;j<=n;j++){//寻找主元
            if(fabs(a[j][i])>fabs(a[maxx][i])) maxx=j;//记得加绝对值
        }
        if(fabs(a[maxx][i])<=eps){//记得加绝对值
            cout<<"No Solution";
            exit(0);
        }
        swap(a[maxx],a[i]);
        for(int j=i+1;j<=n;j++){//加减消元
            double tmp=a[j][i]/a[i][i];
            for(int k=1;k<=n+1;k++){
                a[j][k]-=tmp*a[i][k];
            }
        }
    }
    for(int i=n;i>=1;i--){//回带
        for(int j=i+1;j<=n;j++){
            a[i][n+1]-=a[i][j]*res[j];
        }
        res[i]=a[i][n+1]/a[i][i];
    }
    return;
}

  • 高斯-约旦消元： 在加减消元时对整个方程组都进行一次消元，可以避免回带操作， a i , n + 1 a i , i \frac{a_{i,n+1}}{a_{i,i}} ai,i​ai,n+1​​即为答案

    inline void gauss(){
    for(int i=1;i<=n;i++){
        int maxx=i;
        for(int j=i+1;j<=n;j++){//寻找主元
            if(fabs(a[j][i])>fabs(a[maxx][i])) maxx=j;//记得加绝对值
        }
        if(fabs(a[maxx][i])<=eps){//记得加绝对值
            cout<<"No Solution";
            exit(0);
        }
        swap(a[maxx],a[i]);
        for(int j=1;j<=n;j++){//加减消元
            if(i==j) continue;
            double tmp=a[j][i]/a[i][i];
            for(int k=1;k<=n+1;k++){
                a[j][k]-=tmp*a[i][k];
            }
        }
    }
    return ;
}
  1. 关于解的数量

  • 显而易见的是，当消元结束后，如果主元为0，则存在无解和有无数解的情况

  • 当常数项为0时，因为 0 × x = 0 0\times x=0 0×x=0恒成立，所以有无数解

  • 当常数项不为0时，无论 x x x取何值，等式都不会成立，所以方程无解

  • 在本题中，因为不进行区分，所以只需在过程中判断系数是否为0即可

  1. 注意事项

  • 小数运算存在精度问题，判零时不能判是否为0，而是要判断其是否小于一个 e p s eps eps，一般情况下， e p s eps eps取 1 0 − 7 ∼ 9 10^{-7\sim9} 10−7∼9

  • 判断大小是否为零和选择主元时的数字比较要记得打绝对值，并且double的绝对值时fabs()

  1. 完整代码
  • 普通高斯消元
#include<bits/stdc++.h>
using namespace std;

const int N = 105;
const double eps=1e-9;
int n;
double a[N][N],res[N];
inline void gauss(){
    for(int i=1;i<=n;i++){
        int maxx=i;
        for(int j=i+1;j<=n;j++){//寻找主元
            if(fabs(a[j][i])>fabs(a[maxx][i])) maxx=j;//记得加绝对值
        }
        if(fabs(a[maxx][i])<=eps){//记得加绝对值
            cout<<"No Solution";
            exit(0);
        }
        swap(a[maxx],a[i]);
        for(int j=i+1;j<=n;j++){//加减消元
            double tmp=a[j][i]/a[i][i];
            for(int k=1;k<=n+1;k++){
                a[j][k]-=tmp*a[i][k];
            }
        }
    }
    for(int i=n;i>=1;i--){//回带
        for(int j=i+1;j<=n;j++){
            a[i][n+1]-=a[i][j]*res[j];
        }
        res[i]=a[i][n+1]/a[i][i];
    }
    return;
}

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n+1;j++){
            scanf("%lf",&a[i][j]);
        }
    }
    gauss();
    for(int i=1;i<=n;i++){
        printf("%.2lf\n",res[i]?res[i]:0);
    }
    return 0;
}
  • 高斯-约旦消元
#include<bits/stdc++.h>
using namespace std;

const int N = 105;
const double eps=1e-9;
int n;
double a[N][N];
inline void gauss(){
    for(int i=1;i<=n;i++){
        int maxx=i;
        for(int j=i+1;j<=n;j++){//寻找主元
            if(fabs(a[j][i])>fabs(a[maxx][i])) maxx=j;//记得加绝对值
        }
        if(fabs(a[maxx][i])<=eps){//记得加绝对值
            cout<<"No Solution";
            exit(0);
        }
        swap(a[maxx],a[i]);
        for(int j=1;j<=n;j++){//加减消元
            if(i==j) continue;
            double tmp=a[j][i]/a[i][i];
            for(int k=1;k<=n+1;k++){
                a[j][k]-=tmp*a[i][k];
            }
        }
    }
    return ;
}

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n+1;j++){
            scanf("%lf",&a[i][j]);
        }
    }
    gauss();
    for(int i=1;i<=n;i++){
        double tmp=a[i][n+1]/a[i][i];
        printf("%.2lf\n",tmp?tmp:0);
    }
    return 0;
}
  1. 两种算法的优劣：好像没有什么太多的区别，但是高斯-约旦消元稍微短小好写一些

例题

例题1：P2455 [SDOI2006] 线性方程组

  • 本质也是板子题，不过有一些细节需要注意

  • 由于主播写的有些累了，这里引用一下题解

    普通的高斯-约旦消元存在漏洞：以 hehe_zhou 大佬的hack数据为例：

    2
0 2 3
0 0 0

    答案显然是0，但是我们的程序会说这是-1。

    我们发现，我们可以根据第一个式子来求出第二元，但是程序会用它来算第一元，并且以后算第二元的时候不会再考虑该式，因为我们认为它已经对应第一元了！

    也就是说，我们的答案受到消元顺序影响。

    题解里部分选手当场召唤玄学势力，把消元顺序随机乱改，大概率能出正解。

    这里提供一种新的做法。

    注意到，我们的症结在于把还有用的式子放到了用不上它的地方，并且以后也不来看它是否可用。

    那么，我们就来看它是否可用。

    原本的高斯消元中，我们认为只有 i 之后的式子是可用的，因为我们不用管无解还是无数解，系数为0直接判掉。

    但这里，我们有可能会在系数为0之后继续做下去。这就是受到消元顺序影响的原因。

    那么，可用的就不仅仅是之后的式子，还有之前系数为0的式子。

    于是解决了。

  • 手动膜拜大佬Rui_R orz orz，大佬不要告我侵权拜谢拜谢

#include<bits/stdc++.h>
using namespace std;

const int N = 55;
const double eps=1e-7;
int n;
double a[N][N];

inline double Abs(double a){
    return a<0?-a:a;
} 

inline void gauss(){
    for(int i=1;i<=n;i++){//未知数 
        int maxx=i;
        for(int j=1;j<=n;j++){//枚举所有行,寻找最大的系数当做主元 
            if(Abs(a[j][j])>eps&&j<i) continue;//排除掉在i之前且系数不为0的（已确定） 
            if(Abs(a[j][i])>Abs(a[maxx][i])) maxx=j;
        }
        for(int j=1;j<=n+1;j++){//将选出主元所在的行换到第i行 
            swap(a[i][j],a[maxx][j]);
        } 
        if(Abs(a[i][i])<=eps) continue;//系数为0,先跳过
        for(int j=1;j<=n;j++){//对剩下的每一行进行加减消元 
            if(i==j) continue;
            double tmp=a[j][i]/a[i][i];
            for(int k=i;k<=n+1;k++){
                a[j][k]-=tmp*a[i][k];
            }
        } 
    }
}

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n+1;j++){
            scanf("%lf",&a[i][j]);
        }
    }
    gauss();
    int flag=1;
    for(int i=1;i<=n;i++){//枚举每一个未知数的系数,判断无解和无数解 
        if(Abs(a[i][i])<=eps){
            if(Abs(a[i][n+1])>eps) flag=-1;//系数为0,常数不为零,无解,且必须所有都符合才无解 
            else if(flag!=-1) flag=0;
        }
    }
    if(flag!=1) printf("%d\n",flag);
    else for(int i=1;i<=n;i++){
        //printf("%.10lf %d\n",Abs(a[i][n+1]/a[i][i]),Abs(a[i][n+1]/a[i][i])<=eps);
        printf("x%d=%.2lf\n",i,Abs(a[i][n+1]/a[i][i])<=eps?0:a[i][n+1]/a[i][i]);//用常数除以第i个未知数的最终系数 
    }

    return 0;
}

例题2：P4035 [JSOI2008] 球形空间产生器

  1. 思路：

  • 注意到这里要求球心坐标，很容易想到圆的轨迹方程 (whk乱入)，又发现已经给出了 n + 1 n+1 n+1个点，所以通过半径相等，我们可以列出n个方程，求解n个未知数，高斯消元即可

  • 注意到，高斯消元只能求解标准形式的n元一次方程，所以我们需要帮忙进行一个化简：

    令已知的 n + 1 n+1 n+1个点分别为 P i ( x i , y i , z i , …   ) , i ∈ { i ∈ R ∣ 1 ≤ i ≤ n + 1 } P_i(x_i,y_i,z_i,\dots),i\in\{i\in R|1\le i\le n+1\} Pi​(xi​,yi​,zi​,…),i∈{i∈R∣1≤i≤n+1}

    设圆心 O ( x , y , z , …   ) O(x,y,z,\dots) O(x,y,z,…)

    则有 ( x i − x ) 2 + ( y i − y ) 2 + ( z i − z ) 2 + ⋯ = ( x i − 1 − x ) 2 + ( y i − 1 − y ) 2 + ( z i − 1 − z ) 2 + … (x_i-x)^2+(y_i-y)^2+(z_i-z)^2+\dots=(x_{i-1}-x)^2+(y_{i-1}-y)^2+(z_{i-1}-z)^2+\dots (xi​−x)2+(yi​−y)2+(zi​−z)2+⋯=(xi−1​−x)2+(yi−1​−y)2+(zi−1​−z)2+…

    化简后，得到

    x i 2 + x 2 − 2 x i x + ⋯ = x i − 1 2 + x 2 − 2 x i − 1 x + … x_i^2+x^2-2x_ix+\dots=x_{i-1}^2+x^2-2x_{i-1}x+\dots xi2​+x2−2xi​x+⋯=xi−12​+x2−2xi−1​x+…

    移项后，得到

    x i 2 − x i − 1 2 + y i 2 − y i − 1 2 ⋯ = 2 ( x i x − x i − 1 ) x + 2 ( y i − y i − 1 ) y … x_i^2-x_{i-1}^2+y_i^2-y_{i-1}^2\dots=2(x_ix-x_{i-1})x+2(y_i-y_{i-1})y\dots xi2​−xi−12​+yi2​−yi−12​⋯=2(xi​x−xi−1​)x+2(yi​−yi−1​)y…

  • 这样我们就可以构造出增广矩阵了：

        //a里存的是原始数据,mat是增广矩阵
    for(int i=2;i<=n+1;i++){
        for(int j=1;j<=n;j++){
            mat[i-1][j]=2*(a[i-1][j]-a[i][j]);
            mat[i-1][n+1]+=(a[i-1][j]*a[i-1][j]-a[i][j]*a[i][j]);
        }
    }
  1. 代码：
#include<bits/stdc++.h>
using namespace std;

const int N = 15;
const double eps = 1e-7;

double a[N][N];
double mat[N][N],ans[N];
int n;

inline double Abs(double a){
    return a<0?-a:a;
}

inline void gauss(){
    for(int i=1;i<=n;i++){
        int maxx=i;
        for(int j=i+1;j<=n;j++){
            if(Abs(mat[j][i])>Abs(mat[maxx][i])) maxx=j;
        }
        if(maxx!=i){
            for(int j=1;j<=n+1;j++){
                swap(mat[i][j],mat[maxx][j]);
            }
        }
        for(int j=i+1;j<=n;j++){
            double tmp=mat[j][i]/mat[i][i];
            for(int k=i;k<=n+1;k++){
                mat[j][k]-=mat[i][k]*tmp;
            }
        }
    }
    for(int i=n;i>=1;i--){
        double tmp=mat[i][n+1];
        for(int j=n;j>i;j--){
            tmp-=ans[j]*mat[i][j]; //代入消元,减掉后面数的贡献 
        }
        ans[i]=tmp/mat[i][i];//最后除掉系数即为答案 
    }
    return ;
}

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n+1;i++){
        for(int j=1;j<=n;j++){
            scanf("%lf",&a[i][j]);
        }
    }
    for(int i=2;i<=n+1;i++){
        for(int j=1;j<=n;j++){
            mat[i-1][j]=2*(a[i-1][j]-a[i][j]);
            //cerr<<mat[i-1][j]<<' '; 
            mat[i-1][n+1]+=(a[i-1][j]*a[i-1][j]-a[i][j]*a[i][j]);
        }
        //cerr<<mat[i-1][n+1]<<'\n';
    }
    gauss();
    for(int i=1;i<=n;i++){
        printf("%.3lf",ans[i]);
        if(i==n) puts("");
        else putchar(' ');
    }
    return 0;
}

经过了上述的努力学习，我们终于可以来学点高级的了

高斯消元+期望

例题3: P3232 [HNOI2013] 游走
  1. 说在这里的前面
  • 联系说在前面，这道题就是促使我来学习高斯消元的那道题
  1. 思路：

  • 我们知道，对于每一条边，它都有一个期望经过次数，而 ∑ \sum ∑这个概率 × \times ×其边权即为分值和的期望

  • 题意要求分值和期望的最小值，存在一个贪心思路，即给访问期望次数最大的边赋上最小的边权值

  • 这样，我们要做的事情就成为了求出所有边的期望经过次数，最后排序赋值即可

  • 对于一条边u->v，它被访问的期望次数为 访问 u 的期望次数 × 1 d e g [ u ] + 访问 v 的期望次数 × 1 d e g [ v ] 访问u的期望次数\times \frac{1}{deg[u]}+访问v的期望次数\times \frac{1}{deg[v]} 访问u的期望次数×deg[u]1​+访问v的期望次数×deg[v]1​

  • 所以问题又转换为求访问某个点的期望次数

  • 我们发现，每个点的期望访问次数会由与它相连的点转移而来，即 f u = ∑ v ∈ t o [ u ] f v d e g v f_u=\sum_{v\in to[u]} \frac{f_v}{deg_v} fu​=∑v∈to[u]​degv​fv​​

  • 特别的，当这个点为1时，还要加上它初始的一次访问

  • 则有：

    f i = ∑ v ∈ t o [ u ] f v d e g v   ( 1 < i < n ) f 1 = ∑ v ∈ t o [ 1 ] f v d e g v + 1 f_i=\sum_{v\in to[u]}\frac{f_v}{deg_v} \ (1<i<n)\\ f_1=\sum_{v\in to[1]}\frac{f_v}{deg_v}+1 fi​=v∈to[u]∑​degv​fv​​ (1<i<n)f1​=v∈to[1]∑​degv​fv​​+1

  • 进一步：

    f i − ∑ v ∈ t o [ u ] f v d e g v = 0   ( 1 < i < n ) f 1 − ∑ v ∈ t o [ 1 ] f v d e g v = 1 f_i-\sum_{v\in to[u]}\frac{f_v}{deg_v}=0 \ (1<i<n)\\ f_1-\sum_{v\in to[1]}\frac{f_v}{deg_v}=1 fi​−v∈to[u]∑​degv​fv​​=0 (1<i<n)f1​−v∈to[1]∑​degv​fv​​=1

  • 我们发现，这个转移会受到后继状态的影响，所以需要设 f i f_i fi​为未知数， 1 d e g i \frac{1}{deg_i} degi​1​为系数（不要算n点，因为访问到n点后就结束了，没有贡献），列出n-1个方程进行求解

        --n;//避免第n个点被统计 
    a[1][n+1]=-1.0;//为了求和方便，将方程都乘上-1，使Σ前系数为1
    for(int u=1;u<=n;u++){//初始化，用于高斯消元求解a数组 
        a[u][u]=-1;
        for(int i=fi[u];i;i=e[i].nxt){
            int v=e[i].to;
            if(v!=n+1) a[u][v]=1.0/d[v];
        }
    }
  1. 代码：
#include<bits/stdc++.h>
using namespace std;

const int N = 510,M = 5e5+10;
const double eps=1e-9;

int n,m;
double f[M],a[N][N],ans;
int d[N],u[M],v[M];
struct edge{
    int to,nxt;
}e[M<<1];
int ecnt,fi[N];

inline void add(int x,int y){
    e[++ecnt].to=y;
    e[ecnt].nxt=fi[x];
    fi[x]=ecnt;
    d[x]++;

    e[++ecnt].to=x;
    e[ecnt].nxt=fi[y];
    fi[y]=ecnt;
    d[y]++;
}

inline void gauss(){
    for(int i=1;i<=n;i++){
        int maxx=i;
        for(int j=i+1;j<=n;j++){
            if(fabs(a[j][i])>fabs(a[maxx][i])) 
                swap(maxx,j);
        }
        if(maxx!=i){
            for(int j=1;j<=n+1;j++){
                swap(a[i][j],a[maxx][j]);
            }
        }
        if(fabs(a[i][i])<eps) continue;
        for(int j=1;j<=n;j++){
            if(j!=i){
                double tmp=a[j][i]/a[i][i];
                for(int k=i+1;k<=n+1;k++)
                    a[j][k]-=a[i][k]*tmp;
            }
        }
    }
    for(int i=1;i<=n;i++){
        a[i][n+1]/=a[i][i];
    }
}

int main(){
    scanf("%d%d",&n,&m);
    for(int i=1;i<=m;i++){
        scanf("%d%d",&u[i],&v[i]);
        add(u[i],v[i]);
    }
    --n;//避免第n个点被统计 
    a[1][n+1]=-1.0;
    for(int u=1;u<=n;u++){//初始化，用于高斯消元求解a数组 
        a[u][u]=-1;
        for(int i=fi[u];i;i=e[i].nxt){
            int v=e[i].to;
            if(v!=n+1) a[u][v]=1.0/d[v];
        }
    }
    gauss();
    for(int i=1;i<=m;i++){//来到第i条边的期望次数为从左端点来的期望次数+从右端点来的期望次数 
        f[i]=a[u[i]][n+1]/d[u[i]]+a[v[i]][n+1]/d[v[i]];
    } 
    sort(f+1,f+1+m);//贪心求解答案，给期望次数大的放小编号 
    for(int i=1;i<=m;i++){
        ans+=f[i]*(m-i+1);
    }
    printf("%.3lf",ans);
    return 0;
}
例题4：P3211 [HNOI2011] XOR和路径——高斯消元+期望+位运算
  1. 思路：

  • 看到xor，多少和位运算有点关系

  • 考虑到最终答案中每一位是否为1 只与计算过程中对应位是否为1有关，并且期望具有可加性： E ( a + b ) = E ( a ) + E ( b ) E(a+b)=E(a)+E(b) E(a+b)=E(a)+E(b)，所以我们对每一位分别进行计算

  • 设 f u f_u fu​为从u走到n边权异或和为1的概率

  • 则有转移：

    f u = ∑ ( u , v ) ∈ E , w ( u , v ) = 1 ( 1 − f v ) d e g u + ∑ ( u , v ) ∈ E , w ( u , v = 0 ) f v d e g u ( 1 d e g u 为从 u 点走到 v 的概率 ) f_u=\sum_{(u,v)\in E,w(u,v)=1}\frac{(1-f_v)}{deg_u}+\sum_{(u,v)\in E,w(u,v=0)} \frac{f_v}{deg_u}\\ (\frac{1}{deg_u}为从u点走到v的概率) fu​=(u,v)∈E,w(u,v)=1∑​degu​(1−fv​)​+(u,v)∈E,w(u,v=0)∑​degu​fv​​(degu​1​为从u点走到v的概率)

  • 进一步，有

    f u d e g u = ∑ ( u , v ) ∈ E , w ( u , v ) = 1 ( 1 − f v ) + ∑ ( u , v ) ∈ E , w ( u , v ) = 0 f v f u d e g u + ∑ ( u , v ) ∈ E , w ( u , v ) = 1 f v − ∑ ( u , v ) ∈ E , w ( u , v ) = 0 f v = ∑ ( u , v ) ∈ E , w ( u , v ) = 1 1 f_udeg_u=\sum_{(u,v)\in E,w(u,v)=1}(1-f_v)+\sum_{(u,v)\in E,w(u,v)=0 }f_v\\ f_udeg_u+\sum_{(u,v)\in E,w(u,v)=1}f_v-\sum_{(u,v)\in E,w(u,v)=0 }f_v=\sum_{(u,v)\in E,w(u,v)=1}1 fu​degu​=(u,v)∈E,w(u,v)=1∑​(1−fv​)+(u,v)∈E,w(u,v)=0∑​fv​fu​degu​+(u,v)∈E,w(u,v)=1∑​fv​−(u,v)∈E,w(u,v)=0∑​fv​=(u,v)∈E,w(u,v)=1∑​1

            a[n][n]-=1.0;
        for(int u=1;u<n;u++){
            a[u][u]=d[u];
            for(int i=fi[u];i;i=e[i].nxt){
                int v=e[i].to;
                int w=(e[i].w>>k)&1;//看一下边权在这一位上是否为1
                if(w){
                   // cout<<u<<' '<<n+1<<' '<<a[u][n+1]<<'\n';
                    a[u][n+1]+=1.0;
                    a[u][v]+=1.0;
                    //cout<<d[u]<<' '<<a[u][v]<<' '<<a[u][n+1]<<'\n';
                }else{
                    a[u][v]-=1.0;
                   // cout<<a[u][v]<<'\n';
                }
            }
       }

    消元后，将每一位的 f 1 f_1 f1​乘上这一位为1时对应权值（ 2 i 2^i 2i)，总和即为答案

  1. 代码
#include<bits/stdc++.h>
using namespace std;

const int N = 120;
const int M = 20005;

int n,m;
double a[N][N],f[N];

struct edge{
    int to,nxt,w;
}e[M];
int fi[N],ecnt=0,d[N];
inline void add(int u,int v,int w){
    e[++ecnt].to=v;
    e[ecnt].nxt=fi[u];
    fi[u]=ecnt;
    e[ecnt].w=w;
    d[u]++;
}

inline void gauss(){
    for(int i=1;i<n;i++){//f[n]不用单独求(为0)
        for(int j=1;j<=n;j++){
            if(i==j) continue;
            double tmp=a[j][i]/a[i][i];
            for(int k=1;k<=n+1;k++){
                a[j][k]-=tmp*a[i][k];
            }
        }
    }
}

int main(){
    scanf("%d%d",&n,&m);
    for(int i=1;i<=m;i++){
        int u,v,w;
        cin>>u>>v>>w;
        add(u,v,w);
        if(u!=v) add(v,u,w);
    }
    double ans=0;
    for(int k=0;k<=30;k++){
        memset(a,0,sizeof(a));
        a[n][n]-=1.0;
        for(int u=1;u<n;u++){
            a[u][u]=d[u];
            for(int i=fi[u];i;i=e[i].nxt){
                int v=e[i].to;
                int w=(e[i].w>>k)&1;//看一下边权在这一位上是否为1
                if(w){
                   // cout<<u<<' '<<n+1<<' '<<a[u][n+1]<<'\n';
                    a[u][n+1]+=1.0;
                    a[u][v]+=1.0;
                    //cout<<d[u]<<' '<<a[u][v]<<' '<<a[u][n+1]<<'\n';
                }else{
                    a[u][v]-=1.0;
                   // cout<<a[u][v]<<'\n';
                }
            }
        }
        gauss();//注意位置!!
        ans+=(a[1][n+1]/a[1][1])*(1<<k);//记得要把系数除掉
    }
    printf("%.3lf\n",ans);
    return 0;
}

高斯消元求解异或方程组

例题5：P2447 [SDOI2010] 外星千足虫
  1. 思路

  • 题意很显然可以抽象为一个系数为1/0的异或方程组

  • 求解时，由于异或性质，我们可以直接用n个bitset存储整个方程组，直接对bitset进行xor来消元即可

    bitset

    b.set(pos,val)把b的第pos个二进制位赋值为val

    b.test(pos)询问b的第pos个二进制位是否位0

  1. 代码
#include<bits/stdc++.h>
using namespace std;

const int N = 2e3+10;

int n,m,ans=0;
string s;
bitset<1010>mat[N];

inline void gauss(){
    ans=-1;
    for(int i=1;i<=n;i++){
        int cur=i;
        while(cur<=m&&!mat[cur].test(i)){
            cur++;//寻找不为0系数
        }
        if(cur>m){
            ans=0;
            return ;//有多解
        }
        ans=max(ans,cur);
        if(cur!=i){
            swap(mat[cur],mat[i]);
        }
        for(int j=1;j<=m;j++){
            if(j!=i&&mat[j].test(i)!=0){
                mat[j]^=mat[i];
            }
        }
    }
    return ;
}

int main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin>>n>>m;
    int res;
    for(int i=1;i<=m;i++){
        cin>>s;
        cin>>res;
        for(int j=0;j<n;j++){
            mat[i].set(j+1,s[j]=='1');
        }
        mat[i].set(0,res);//系数存在第0位
    }
    gauss();
    if(ans){
        cout<<ans<<'\n';
        for(int i=1;i<=n;i++){
            if(mat[i].test(0)){//由于系数只能是1/0,所以不用除系数
                cout<<"?y7M#\n";
            }else{
                cout<<"Earth\n";
            }
        }
    }else{
        cout<<"Cannot Determine";
    }
    return 0;
}
jyy0806
关注
高斯消元:线性方程组: 高斯消元-matlab开发
05-29
高斯消元，也称为行约简，是线性代数中的一种算法，用于求解线性方程组。 通常理解为对相应系数矩阵执行的一系列操作。 该方法还可用于求矩阵的秩、计算矩阵的行列式以及计算可逆方阵的逆。 该方法以 Carl Friedrich...
gaussian_eliminatio​n:高斯消元和 Gaussian-Jordan 消元-matlab开发
05-29
总结来说，高斯消元和高斯-约旦消元是线性代数中的基础工具，MATLAB通过简洁的代码实现了这些方法，使得线性方程组的求解变得简单易行。通过理解和掌握这些算法，不仅能够提高MATLAB编程能力，也能加深对线性代数的...
【数学】高斯-约旦消元
OI博客
05-03 1309
讲解基本高斯-约旦消元
笔记---高斯消元
2302_79440616的博客
01-31 603
但是考虑到本题作为一道模板题，考察点并不在于此，在此处卡住大多同学的代码没有太大意义，故增加 SPJ，对输出。之后继续枚举，但是被换到最上面的那一行就不要动了，直至最后。在数学中，一般没有正零或负零的概念，所以严格来说应当输出。通过高斯消元，来使方程变为上三角形式，从而可以解得方程（大学线性代数）输入一个包含 n 个方程 n 个未知数的线性方程组。如果消过之后，是完美的上三角形，那么会得到唯一解。个未知数的解，结果保留两位小数。个实数，表示一个方程的。（非零）的方程，那么无解。的方程，那么就有无穷多解。
数论 - 高斯消元算法
林深不见鹿 的博客
07-27 2168
1.高斯消元 (1)定义 高斯消元法是求解线性方阵组的一种算法，它也可用来求矩阵的秩，以及求可逆方阵的逆矩阵。它通过逐步消除未知数来将原始线性系统转化为另一个更简单的等价的系统。它的实质是通过初等行变化将线性方程组的增广矩阵转化为行阶梯矩阵. (2)初等行变换 初等行变换 定义：所谓数域P上矩阵的初等行变换是指下列3种变换： (1）以P中一个非零的数乘矩阵的某一行 (2）把矩阵的某一行的c倍加到另一行，这里c是P中的任意一个数 (3）互换矩阵中两行的位置 一般来说，一个矩阵经过初等行变换后就变成了另一个矩
数论 - 高斯消元
2203_75327677的博客
02-20 690
数学上，高斯消元法（或译：高斯消去法），是线性代数规划中的一个算法，可用来为线性方程组求解
数论----高斯消元
Lin_disguiser的博客
08-05 907
学了一天的高斯消元，又退了两天，才接着补坑，唉~~自己为什么这么不争气~~ 主要的学习高斯消元的来源还是论文---何江舟的《高斯消元解线性方程组》 注意几点： 1.equ和var分别代表方程数和未知数 2.在代码中注意k和col的实时变化，k循环结束后表示消元后的最后一个行，col循环后的为第var-1列（0~var） 3.注意无穷解的代码段 4.temp=a[i][var];
计算机编程 高斯消元,高斯-若尔当消元法
weixin_39896617的博客
07-21 2940
高斯－若尔当消元法(英语：Gauss-Jordan Elimination)，或译为高斯-约旦消元法，简称G-J消元法，是数学中的一个算法，是高斯消元法的另一个版本。它在线性代数中用来找出线性方程组的解，其方法与高斯消去法相同。唯一相异之处就是这算法产生出来的矩阵是一个简化行梯阵式，而不是高斯消元法中的行梯阵式。相比起高斯消元法，此算法的效率比较低，却可把方程组的解用矩阵一次过表示出来。中文名高斯...
高斯消元-线性方程组求解
w75779的博客
09-24 559
高斯消元
数论算法-高斯消元
m0_53792316的博客
04-13 1022
编程中高斯消元的学习及使用。
高斯消元-异或版
qq_40329819的博客
11-07 578
题目 2020icpc济南-A-Matrix Equation 思路 矩阵 A 与答案每列乘积异或和互不干扰，故对每列分别求解，等式右边化为 0 找到当前处理的未知数为 1 的行，交换，对之后该未知数为 1 的行整行异或，化为上三角求出解 代码 #include<bits/stdc++.h> #define LL long long #define mod 998244353 using namespace std; const int N=209; int n; LL ans=1LL; LL
线性方程组：高斯消元：线性方程组：高斯消元-matlab开发
05-30
下面将详细介绍如何使用MATLAB进行高斯消元求解线性方程组的步骤： 1. **创建系数矩阵与常数向量**： 首先，我们需要构建线性方程组的系数矩阵A和常数向量b。例如，如果我们有以下线性方程组： \[ \begin{align*}...
高斯消元java源码-Lights-Out:熄灯
06-06
高斯消元java源码熄灯 Lights Out 的游戏从一个 nxn 的按钮板开始，它们可以兼作灯泡，其中 n ≥ 2。在我们的实现中，我们从打开所有灯开始，目标是关闭它们。 问题是当一个按钮被切换时，它的紧邻的邻居也会被切换...
C语言[函数实现模块化程序设计]
2401_87108776的博客
10-09 458
函数的名字函数的类型函数参数的名字及类型函数的功能。
【C++刷题】力扣-#121-买卖股票的最佳时机
会写代码的饭桶
10-16 323
给定一个数组 prices，其中 prices[i] 表示第 i 天的股票价格。假设你可以在第 i 天买入并在第 j 天卖出股票（i ≤ j），设计一个算法来计算你所能获取的最大利润。注意你只能持有一股股票，并且你不能同时参与多笔交易（即在再次买入前必须卖出股票）。
C++中的vector介绍（常用函数）
zxybf_的博客
10-15 752
vector的文档介绍vector是序列容器，表示可以改变大小的数组。(constructor)构造函数声明接口说明vector()（重点）无参构造构造并初始化n个val（重点）拷贝构造使用迭代器进行初始化构造int main()//无参的构造//n个val的构造//拷贝构造//迭代器区间构造return 0;
[C++ 核心编程]笔记 5 函数提高
Alex的博客
10-11 446
作用:函数名可以相同，提高复用性同一个作用域下函数名称相同函数参数类型不同 或者 个数不同 或者 顺序不同注意: 函数的返回值不可以作为函数重载的条件//原函数cout
C++中的CRTP
最新发布
moter的专栏
10-18 568
我们知道，要定义派生类肯定需要知道基类的类型，而在使用类模板来实例化一个模板类时又得需要知道模板的实参类型，在这里实例化基类时的实参类型又是派生类的类型。CRTP基类是类模板，基类实例化之后，和派生类是一一对应的，即一个基类只能派生一个子类，而模板参数化this指针的实现方式是基类只有一个，而且是非模板类，它的派生类可以有多个，但是每定义一个派生类，在该基类中都会使用派生类类型实例化出一个对应的interface()函数，也就是在基类中会有多个interface()成员函数和它的派生类一一对应。
C++：继承
WG_17的博客
10-13 878
下⾯我们看到Person是基类，也称作⽗类。Student是派⽣类，也称作⼦类。(因为翻译的原因，所以 既叫基类/派⽣类，也叫⽗类/⼦类)
高斯消元与高斯-约旦算法教程
高斯-约当消元法可以直观地得到矩阵的逆，因为通过此算法将增广矩阵（原矩阵与单位矩阵拼接形成的矩阵）的原矩阵部分消元到单位矩阵，单位矩阵部分自然就变成了原矩阵的逆。 在本教程中，可能会包含以下几个方面的...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值