• 模板


    真丶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("
    ");
            }*/
    
            /*当不需要关心方程右边时,删去以下部分*/
            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的?要是是常数项就完了
                    ;
                }
                //消去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
    ",ans[i]);
        }
    }
    
    
    

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

    ```cpp #include #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("
    ");
        }*/
    
        /*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 ",curcount,curmoney);
    for(int i=1; i<=n; i++)
    printf(" %d",tx[i]);
    printf(" ------- ");
    /

        /*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
    ",minprice);
    }
    else{
        //给自由变量赋值
        //dfs(n,0,0);
    
        //直接得到sumprice
        int sumprice=0;
        for(int i=1;i<=r;i++)
            sumprice+=price[i];
    
        printf("%d %d
    ",r,sumprice);
    }
    

    }

    </details>
    
    ##double的n*n高斯消元法
    
    虽然经过模板验证,但是不明白为什么可以中途退出。假设我们用来插值,第一个是常数,当某列全为0,那么不是只跟当前变量没有关系吗?建议使用模意义的模板,有数学证明。(虽然这个交上去没有翻过车)
    
    ```cpp
    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
    ",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
    ",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("
    ");*/
    
            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("
    ");
            }
            printf("---
    ");*/
    
            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("
    ");
            }
            printf("---
    ");*/
    
    
            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("
    ");*/
            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("
    ");*/
                    return;
                }
            }
            res=min(res,cur);
            /*printf("res=%d
    ",res);
            for(int i=1; i<=n; i++)
                printf("%d ",tx[i]);
            printf("
    ");*/
            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
    ",res);
    }
    
    
  • 相关阅读:
    Coursera Algorithms Programming Assignment 2: Deque and Randomized Queue (100分)
    Coursera Algorithms week1 查并集 练习测验:3 Successor with delete
    Coursera Algorithms week1 查并集 练习测验:2 Union-find with specific canonical element
    项目选题报告答辩总结
    项目UML设计(团队)
    第七次作业--项目需求分析(团队)
    结对项目--第二次作业
    【软件工程实践 · 团队项目】 第二次作业
    第五次作业--原型设计(结对)
    团队项目系列博客 —— 在路上(之wampserver 修改根目录以及配置多站点以及修改端口号)
  • 原文地址:https://www.cnblogs.com/Yinku/p/10705687.html
Copyright © 2020-2023  润新知