模板 - 高斯约旦消元法

Posted 2021-10-25 yinku

真丶long double高斯约旦消元法

eps需要取得大一些，以免增加了矩阵的秩。
long double可能会慢一些但是无所谓，被卡精度太恶心了。

需要知道一些线代的知识（线代67说你呢！），比如秩、极大线性无关组（线性基）之类的。

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

namespace Gauss_Jordan_Elimination {
    const int MAXN=505;
    const int MAXM=505;
    long double a[MAXN][MAXM+1];
    double ans[MAXN];

    const long double eps=1e-6;

    //返回增广矩阵的秩，-1表示无解
    int gauss_jordan(int n,int m) {
        int r=0;
        for(int i=1; i<=m; i++) {
            //当前是第i列
            int mx=-1;
            //从当前秩的下一行开始往下数
            for(int j=r+1; j<=n; j++)
                if(mx==-1||fabs(a[j][i])>fabs(a[mx][i])){
                    mx=j;
                }
            if(mx==-1) {
                //该列全0，被跳过
                continue;
            }
            r++;
            //增加一个线性基，当前秩增加
            if(mx!=r) {
                //需要交换
                for(int j=1; j<=m+1; j++)
                    //m+1表示增广矩阵
                    swap(a[r][j],a[mx][j]);
                //交换行
            }

            for(int j=1; j<=n; j++) {
                //枚举每一行
                if(j!=r&&fabs(a[j][i])>eps) {
                    //消去除了r以外的其他行
                    long double tmp=a[j][i]/a[r][i];
                    //该行系数是当前列的tmp倍
                    for(int k=i; k<=m+1; k++) {
                        //把当前列对应行位置扩大tmp倍，然后用每个元素减去，由高斯约当的过程，左边的绝对是0
                        a[j][k]-=a[r][k]*tmp;
                    }
                }
            }
        }

        /*for(int i=1; i<=n; i++) {
            for(int j=1; j<=m; j++) {
                printf(" %.2f ",a[i][j]);
            }
            printf("\n");
        }*/

        /*当不需要关心方程右边时，删去以下部分*/
        int j=1;
        for(int i=1; i<=r; i++){
            //找第一个非0值，也就是阶梯型矩阵的第一个值，因为是高斯约当消元所以其他位置都是0
            while(fabs(a[i][j])<eps)
                j++;

            if(fabs(a[i][m+1])<eps){
                //在秩里面，居然有右侧为0的？
                return -1;
            }
            //消去x的系数，系数化为1
            ans[i]=(double)(a[i][m+1]/a[i][j]);

            j++;
            if(j>m){
                break;
            }
        }

        return r;
    }
}

using namespace Gauss_Jordan_Elimination;

int n,m;

int main() {
    scanf("%d",&n);
    m=n;
    for(int i=1; i<=n; i++) {
        for(int j=1; j<=m+1; j++) {
            double tmp;
            scanf("%lf",&tmp);
            a[i][j]=tmp;
        }
    }

    int r=gauss_jordan(n,m);

    if(r!=n){
        puts("No Solution");
    }
    else{
        for(int i=1;i<=r;i++)
            printf("%.2f\n",ans[i]);
    }
}

https://www.luogu.org/problemnew/show/P3265

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

/**
 *某个线性方程组的矩阵的秩是确定的
 *用高斯约旦消元可以求出矩阵的秩
 *明显可以贪心，每次选择最便宜的非零向量？
 *当矩阵为零矩阵时，可以买最便宜的零向量
 */

int price[505];

int r=0;

namespace Gauss_Jordan_Elimination {
    const int MAXN=505;
    long double a[MAXN][MAXN];
    long double ans[MAXN];

    const long double eps=1e-6;

    bool gauss_jordan(int n,int m) {
        r=0;
        for(int i=1; i<=m; i++) {
            //当前是第i列
            int mx=-1;
            //从当前秩的下一行开始往下数
            for(int j=r+1; j<=n; j++)
                if(fabs(a[j][i])>eps){
                    if(mx==-1||price[j]<price[mx]){
                        mx=j;
                    }
                }
            if(mx==-1) {
                //该列全0，被跳过
                continue;
            }
            r++;
            //增加一个线性基，当前秩增加
            if(mx!=r) {
                //最小花费的行需要交换
                for(int j=1; j<=m; j++)
                    swap(a[r][j],a[mx][j]);
                //交换行
                swap(price[r],price[mx]);
                //这两个装备交换了位置，价格当然也变了
            }

            for(int j=1; j<=n; j++) {
                //枚举每一行
                if(j!=r&&fabs(a[j][i])>eps) {
                    //消去除了r以外的其他行
                    long double tmp=a[j][i]/a[r][i];
                    //该行系数是当前列的tmp倍
                    for(int k=i; k<=m; k++) {
                        //把当前列对应行位置扩大tmp倍，然后用每个元素减去，由高斯约当的过程，左边的绝对是0
                        a[j][k]-=a[r][k]*tmp;
                    }
                }
            }
        }

        /*for(int i=1; i<=n; i++) {
            for(int j=1; j<=m; j++) {
                printf(" %.2f ",a[i][j]);
            }
            printf("\n");
        }*/

        /*for(int i=1; i<=r; i++)
            //xi的系数化为1
            ans[i]=a[i][m+1]/a[i][i];*/

        return true;
    }
}

using namespace Gauss_Jordan_Elimination;

int n,m;

int ansmoney=0x3f3f3f3f;
int anscount=0;

int tx[505];

bool iszero(int id) {
    int j=id;
    while(j<=m) {
        if(fabs(a[id][j])>eps) {
            return false;
        } else {
            j++;
        }
    }
    if(j>m) {
        return true;
    }
    return false;
}

void dfs(int dep,int curcount,int curmoney) {
    if(curcount<=anscount&&curmoney>=ansmoney)
        //花更多的钱买了更少的装备，脑子有洞
        return;
    if(dep==0) {
        /*printf("cnt=%d mon=%d\n",curcount,curmoney);
        for(int i=1; i<=n; i++)
            printf(" %d",tx[i]);
        printf("\n-------\n");*/

        /*for(int i=r+1; i<=n; i++) {
            //自由变量能否被表示
            double tmp=0.0;
            for(int j=r+1; i<=n; i++) {
                if(tx[j]>0&&fabs(a[dep][j])>eps) {
                    //该自由变量存在，且属性与这个装备有关
                    tmp+=a[dep][j];
                }
                if(fabs(tmp)<eps) {
                    //该变量能被线性表示
                    tx[dep]=0;
                    dfs(dep-1,curcount,curmoney);
                } else {
                    cout<<"dep="<<dep<<" can't be e"<<endl;
                    tx[dep]=1;
                    dfs(dep-1,curcount+1,curmoney+price[dep]);
                }
            }*/

        if(curcount>anscount) {
            anscount=curcount;
            ansmoney=curmoney;
        } else if(curcount==anscount) {
            ansmoney=min(ansmoney,curmoney);
        }
        return;
    }
    if(!iszero(dep)) {
        //cout<<"dep="<<dep<<" can't be e"<<endl;
        tx[dep]=1;
        dfs(dep-1,curcount+1,curmoney+price[dep]);
    } else {
        //cout<<"dep="<<dep<<" "<<"iszero"<<endl;
        tx[dep]=0;
        dfs(dep-1,curcount,curmoney);
    }
}

int main() {
    scanf("%d%d",&n,&m);
    int minprice=0;
    for(int i=1; i<=n; i++) {
        for(int j=1; j<=m; j++) {
            double tmp;
            scanf("%lf",&tmp);
            a[i][j]=tmp;
        }
    }

    for(int i=1; i<=n; i++) {
        scanf("%d",&price[i]);
        minprice=min(minprice,price[i]);
    }

    gauss_jordan(n,m);

    //已经求矩阵的秩
    //cout<<r<<endl;

    if(r==0){
        printf("0 %d\n",minprice);
    }
    else{
        //给自由变量赋值
        //dfs(n,0,0);

        //直接得到sumprice
        int sumprice=0;
        for(int i=1;i<=r;i++)
            sumprice+=price[i];

        printf("%d %d\n",r,sumprice);
    }
}

double的n*n高斯消元法

虽然经过模板验证，但是不明白为什么可以中途退出。假设我们用来插值，第一个是常数，当某列全为0，那么不是只跟当前变量没有关系吗？建议使用模意义的模板，有数学证明。（虽然这个交上去没有翻过车）

namespace Gauss_Jordan_Elimination{
    const int MAXN=105;
    double a[MAXN][MAXN],del;
    double ans[MAXN];

    const double eps=1e-10;

    bool gauss_jordan(int n){
        for(int i=1;i<=n;i++){
            int mx=i;
            for(int j=i+1;j<=n;j++)
                if(fabs(a[j][i])>fabs(a[mx][i]))
                    mx=j;
            for(int j=1;j<=n+1;j++)
                swap(a[i][j],a[mx][j]);
            if(fabs(a[i][i])<eps)
                return false;
            for(int j=1;j<=n+1;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++)
            ans[i]=a[i][n+1]/a[i][i];

        return true;
    }
}


using namespace Gauss_Jordan_Elimination;

模意义下高斯消元法

这里是模质数的高斯消元，假如不是质数，乘法逆元都不一定存在。
https://codeforces.com/contest/1155/problem/E

顺便也发现了一种插值的策略，直接高斯消元然后求出多项式。
这里的高斯消元好像和上面有点不同，但是上面好像是书上类似的写法。
顺带一提，交互题可以写个jury函数，通过传入一个询问来返回一个回答。

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

const ll p=1000003;

ll inv[1000005];

inline ll mp_init(ll x) {
    if(x>=0&&x<p)
        return x;
    x%=p;
    if(x<0)
        x+=p;
    return x;
}

inline ll mp_add(ll x,ll y) {
    return mp_init(x+y);
}

inline ll mp_sub(ll x,ll y) {
    return mp_init(x-y);
}

inline ll mp_mul(ll x,ll y) {
    return (ll)x*y%p;
}

inline ll qpow(ll x,ll n) {
    ll res=1%p;
    while(n) {
        if(n&1)
            res=res*x%p;
        x=x*x%p;
        n>>=1;
    }
    return res;
}

inline ll get_inv(ll x) {
    return qpow(x,p-2);
}

inline ll mp_div(ll x,ll y) {
    return (ll)x*inv[y]%p;
}


namespace Gauss_Jordan_Elimination_Modulo_Prime {
    const int MAXN=55;
    ll a[MAXN][MAXN];
    ll ans[MAXN];

    bool gauss_jordan(int n) {
        //化为模意义下的数
        for(int i=1; i<=n; i++)
            for(int j=1; j<=n+1; j++)
                a[i][j]=mp_init(a[i][j]);

        for(int i=1; i<=n; i++) {
            //当前是第i列
            int mx=i;
            //当前列的最大变量系数所在的行
            for(int j=i+1; j<=n; j++)
                if(a[j][i]>a[mx][i])
                    mx=j;
            if(a[mx][i]==0) {
                //当当前列最大的行的变量系数都为0，则不影响
                continue;
            }
            if(mx!=i) {
                for(int j=1; j<=n+1; j++)
                    swap(a[i][j],a[mx][j]);
                //交换行，使得当前位置系数最大
            }

            //第j行
            for(int j=1; j<=n+1; j++) {
                if(j!=i) {
                    //消去其他行
                    ll tmp=mp_div(a[j][i],a[i][i]);
                    //该行系数是当前列的tmp倍
                    for(int k=i+1; k<=n+1; k++) {
                        //把这里的k从1开始，可以看见完全消除的结果，但是没必要
                        //把当前列对应行位置扩大tmp倍，然后用每个元素减去
                        a[j][k]=mp_sub(a[j][k],mp_mul(a[i][k],tmp));
                    }
                }
            }
        }


        for(int i=1; i<=n; i++)
            //xi的系数化为1
            ans[i]=mp_div(a[i][n+1],a[i][i]);

        return true;
    }
}


using namespace Gauss_Jordan_Elimination_Modulo_Prime;


ll jury(ll x) {
    ll a[11]= {2,9,1,4,8,10,0,2,8,7,6};
    ll ret=0;
    ll tx=1;
    for(int i=0; i<=10; i++) {
        ret=(ret+tx*a[i]%p)%p;
        tx=tx*x%p;
    }
    return ret;
}

int query(int x) {
    printf("? %d\n",x);
    fflush(stdout);
    scanf("%d",&x);
    //x=jury(x);
    return x;
}

void set_matrix(int id,ll x,ll y) {
    a[id][1]=1;
    for(int i=2; i<=11; i++) {
        a[id][i]=mp_mul(a[id][i-1],x);
    }
    a[id][12]=mp_init(y);
}

ll calc(ll x) {
    ll ret=0;
    ll tx=1;
    for(int i=1; i<=11; i++) {
        ret=(ret+tx*ans[i]%p)%p;
        tx=tx*x%p;
    }
    return ret;
}


void answer(int x) {
    printf("! %d\n",x);
}


int main() {
    for(int i=1;i<1000003;i++){
        inv[i]=get_inv(i);
    }

    for(int i=0; i<=10; i++) {
        int res=query(i);
        if(res==0) {
            /*answer(i);
            exit(0);*/
        }
        set_matrix(i+1,i,res);
    }

    int res=gauss_jordan(11);

    if(res==0) {
        answer(-1);
    } else {

        /*for(int i=1;i<=11;i++){
            printf("%lld ",ans[i]);
        }
        printf("\n");*/

        for(int i=0; i<1e6+3; i++) {
            int r=(int)calc(i);
            if(r==0) {
                answer(i);
                exit(0);
            }
        }
        answer(-1);
    }
}

模2意义下高斯消元法

模2的话，加减法变成异或，乘法变成与，没有除法了（除法按道理说，除以1就是与1，保持不变，除以0……说明这个是自由变量）。
顺便给一个带自由变量的高斯消元，高斯求出来的只是其中一个解，我们可以通过高斯消元消出来的矩阵简化计算，在这种运用中，建议把高斯消元的起点变为1，直观看见消元的结果。

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


namespace Gauss_Jordan_Elimination_Modulo_2 {
    const int MAXN=40;
    int a[MAXN][MAXN];
    int ans[MAXN];

    bool gauss_jordan(int n) {
        //化为模意义下的数
        /*for(int i=1; i<=n; i++)
            for(int j=1; j<=n+1; j++)
                a[i][j]=a[i][j]&1;*/

        /*for(int i=1; i<=n; i++) {
            for(int j=1; j<=n+1; j++) {
                printf("%d ",a[i][j]);
            }
            printf("\n");
        }
        printf("---\n");*/

        for(int i=1; i<=n; i++) {
            //当前是第i列
            int mx=i;
            for(int j=i; j<=n; j++)
                if(a[j][i]) {
                    mx=j;
                    break;
                }
            if(a[mx][i]==0) {
                //当当前列最大的行的变量系数都为0，则不影响
                continue;
            }
            if(mx!=i) {
                for(int j=1; j<=n+1; j++)
                    swap(a[i][j],a[mx][j]);
                //交换行，使得当前位置系数最大
            }

            //第j行
            for(int j=1; j<=n+1; j++) {
                if(j!=i&&a[j][i]) {
                    //消去其他行
                    for(int k=/*i+*/1; k<=n+1; k++) {
                        a[j][k]^=a[i][k];
                    }
                }
            }
        }


        /*for(int i=1; i<=n; i++) {
            for(int j=1; j<=n+1; j++) {
                printf("%d ",a[i][j]);
            }
            printf("\n");
        }
        printf("---\n");*/


        for(int i=1; i<=n; i++)
            //xi的系数化为1
            ans[i]=a[i][n+1]&a[i][i];



        /*for(int i=1; i<=n; i++)
            printf("%d ",ans[i]);
        printf("\n");*/
        return true;
    }
}


using namespace Gauss_Jordan_Elimination_Modulo_2;

int n,m,res=0;

int tx[40];

void dfs(int x,int cur) {
    //cout<<"x="<<x<<" res="<<cur<<endl;
    if(cur>=res)
        //最优性剪枝
        return;
    if(x==0) {
        //经过若干调整到达终点，需要检验
        //cout<<"cur="<<cur<<endl;
        for(int i=1; i<=n; i++) {
            int y=0;
            for(int j=1; j<=n; j++) {
                y^=a[i][j]*tx[j];
            }
            //如果这个灯没有前往他要去的目标状态
            if(y!=a[i][n+1]) {
                /*cout<<"fail on r " <<i<<endl;
                for(int i=1; i<=n; i++)
                    printf("%d ",tx[i]);
                printf("\n");*/
                return;
            }
        }
        res=min(res,cur);
        /*printf("res=%d\n",res);
        for(int i=1; i<=n; i++)
            printf("%d ",tx[i]);
        printf("\n");*/
        return;
    }
    if(a[x][x]) {
        //不是自由变量，取值由左下角已经确定了的值来确定
        int y=a[x][n+1];
        //该变量被右端函数值以及前面以及“确定”的值约束
        for(int i=x+1; i<=n; i++) {
            //其实就是y-=a[i]x[i]，把已经确定了的x的值乘上他的贡献因数然后减掉
            y^=a[x][i]&tx[i];//第i个变量当前取值是tx，那么现在的x要取什么才会等于y？
            //y=0说明其他变量凑够了
            //cout<<"y="<<y<<endl;
        }
        int ttx=tx[x];
        tx[x]=y;
        dfs(x-1,cur+y);
        tx[x]=ttx;
    } else {
        //当前是自由变量，取值是任意的
        tx[x]=0;
        dfs(x-1,cur);
        //原本的高斯消元后，自由变量的取值都被归一化步骤赋值0，现在赋值1，不一定满足右侧全1，最后需要检验
        tx[x]=1;
        dfs(x-1,cur+1);
        tx[x]=0;
    }
}

int main() {
    scanf("%d%d",&n,&m);
    for(int i=1; i<=m; i++) {
        int u,v;
        scanf("%d%d",&u,&v);
        a[u][v]=1;
        a[v][u]=1;
    }
    for(int i=1; i<=n; i++)
        a[i][i]=a[i][n+1]=1;
    gauss_jordan(n);

    res=0;
    for(int i=1; i<=n; i++) {
        res+=ans[i];
    }

    dfs(n,0);

    printf("%d\n",res);
}