「模板」高斯消元

目录

• 模板

模板

因为 \(\text{gauss}\) 消元法没有什么好讲的。

主要知识使用了小学的消元技巧，主要注意细节。

但是小学知识大家都应该知道吧...

所以直接上板子吧...为封装

#include<cstdio>
#include<cstring>

#define rep(i,__l,__r) for(register int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(register int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define uint unsigned int
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long
// #define int unsigned

#ifdef FILEOI
    #define MAXBUFFERSIZE 500000
    inline char fgetc(){
        static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
        return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
    }
    #undef MAXBUFFERSIZE
    #define cg (c=fgetc())
#else
    #define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
    char c;bool f=0;
    while(cg<'0'||'9'<c)f|=(c=='-');
    for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
    if(f)x=-x;
}
inline int qread(){
    int x=0;char c;bool f=0;
    while(cg<'0'||'9'<c)f|=(c=='-');
    for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
    return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
    inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
    if(x<0)return (void)(putchar('-'),fwrit(-x));
    if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
    return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}

int n;

namespace gauss{
    const int MAXN=100;
    const double eps=1e-20;
    struct Equation{
        double dat[MAXN+5];
        inline Equation(){
            memset(dat,0,sizeof dat);
        }
        inline double& operator [](const int __index){
            return dat[__index];
        }
        inline Equation operator *(const double __x)const{
            Equation ret=*this;
            rep(i,1,n+1)ret[i]*=__x;
            return ret;
        }
        inline Equation operator /(const double __x)const{
            Equation ret=*this;
            rep(i,1,n+1)ret[i]/=__x;
            return ret;
        }
        inline Equation operator -(Equation x)const{
            Equation ret=*this;
            rep(i,1,n+1)ret[i]-=x[i];
            return ret;
        }
        inline Equation operator *=(const double __x){return (*this)=(*this)*__x;}
        inline Equation operator /=(const double __x){return (*this)=(*this)/__x;}
        inline Equation operator -=(Equation x){return (*this)=(*this)-x;}
    }data[MAXN+5];
    int ans[MAXN+5];
    bool vis[MAXN+5];
    inline int get(const int x){
        //得到 x_i 的最大的系数的那个等式
        double maxx=0;int ret;
        rep(i,1,n)if(!vis[i] && maxx<fab(data[i][x]))
            maxx=fab(data[i][x]),ret=i;
        vis[ret]=true;
        return ret;
    }
    inline void launchGauss(){
        //高斯消元的主过程
        rep(i,1,n){
            int now=get(i);
            const double a=data[now][i];
            rep(j,1,n)if(j!=now){
                if(fab(data[j][i])<eps){
                //避免除 0
                    printf("No Solution\n");
                    return;
                }
                data[j]*=(a/data[j][i]);
                data[j]-=data[now];
                //此处用高斯消元消掉第 i 个未知数
            }
            ans[i]=now;
        }
        rep(i,1,n)if(fab(data[ans[i]][i])<eps){
        //检测是否有几个方程有线性相关性, 将我们要求的未知数消掉了
            printf("No Solution\n");
            return;
        }
        rep(i,1,n){
            //因为只输出部分小数, 所以有可能会出现输出 -0.00 的情况
            if(fab(data[ans[i]][n+1]/data[ans[i]][i])<eps)printf("0.00");
            else printf("%.2lf\n",data[ans[i]][n+1]/data[ans[i]][i]);
        }
    }
};

using namespace gauss;

signed main(){
#ifdef FILEOI
    freopen("file.in","r",stdin);
    freopen("file.out","w",stdout);
#endif
    qread(n);
    rep(i,1,n)rep(j,1,n+1)scanf("%lf",&data[i][j]);
    launchGauss();
    return 0;
}