高斯消元

板子

#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#define maxn 110
using namespace std;
template<typename T>
inline void read(T &x){
    x=0;bool flag=0;char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') flag=1;
    for(;isdigit(c);c=getchar()) x=x*10+(c^48);
    if(flag) x=-x;
}

double eps=1e-8;
int n;
double a[maxn][maxn];
bool no,many;

void gauss(){
    for(int i=1;i<=n;i++){
        int maxi=i;
        for(int j=i+1;j<=n;j++){
            if(fabs(a[j][i])>fabs(a[maxi][i])) maxi=j;
        }
        for(int j=1;j<=n+1;j++) swap(a[maxi][j],a[i][j]);
        if(fabs(a[i][i])<eps) continue;
        double tmp=a[i][i];//后面a[i][i]会变的，所以不能直接除以a[i][i]，必须拿tmp存一下 
        for(int j=1;j<=n+1;j++) a[i][j]/=tmp;
        for(int j=1;j<=n;j++){
            if(i==j) continue;
            double tmp=a[j][i];
            for(int k=1;k<=n+1;k++) a[j][k]-=a[i][k]*tmp;
        }
    }
    for(int i=1;i<=n;i++){
        int cnt=0;
        for(int j=1;j<=n+1;j++) if(fabs(a[i][j])<eps) cnt++;
        if(cnt==n&&fabs(a[i][n+1])>=eps) no=1;
        if(cnt==n+1) many=1;
        
    }
	if(no||many) return ;
    for(int i=n-1;i;i--)
        for(int j=i+1;j<=n;j++)
            a[i][n+1]-=a[j][n+1]*a[i][j];
}

int main(){
    read(n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            read(a[i][j]);
    gauss();
    if(no){
        cout<<"No Solution"<<endl;
        return 0;
    }
    if(many){
        cout<<"No Solution"<<endl;
        return 0;
    }
    for(int i=1;i<=n;i++) printf("%.2f
",fabs(a[i][n+1])<eps?0:a[i][n+1]);
    return 0;
}

例题

Luogu P2455线性方程组

题目

#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#define maxn 60
using namespace std;
template<typename T>
inline void read(T &x){
    x=0;bool flag=0;char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') flag=1;
    for(;isdigit(c);c=getchar()) x=x*10+(c^48);
    if(flag) x=-x;
}

double eps=1e-8;
int n;
double a[maxn][maxn];
bool no,many;
bool no_flag,inf_flag;

void gauss(){
    for(int i=1;i<=n;i++){
        int maxi=i;
        for(int j=i+1;j<=n;j++){
            if(fabs(a[j][i])>fabs(a[maxi][i])) maxi=j;
        }
        for(int j=1;j<=n+1;j++) swap(a[maxi][j],a[i][j]);
        if(fabs(a[i][i])<eps) continue;
        double tmp=a[i][i];//后面a[i][i]会变的，所以不能直接除以a[i][i]，必须拿tmp存一下 
        for(int j=1;j<=n+1;j++) a[i][j]/=tmp;
        for(int j=1;j<=n;j++){
            if(i==j) continue;
            double tmp=a[j][i];
            for(int k=1;k<=n+1;k++) a[j][k]-=a[i][k]*tmp;
        }
    }
    for(int i=1;i<=n;i++){
        int cnt=0;
        for(int j=1;j<=n+1;j++) if(fabs(a[i][j])<eps) cnt++;
        if(cnt==n&&fabs(a[i][n+1])>=eps) no=1;
        if(cnt==n+1) many=1;
        
    }
	if(no||many) return ;
    	for(int i=n-1;i;i--)
            for(int j=i+1;j<=n;j++)
                a[i][n+1]-=a[j][n+1]*a[i][j];
}

int main(){
    read(n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            read(a[i][j]);
    gauss();
    if(no){
        cout<<-1<<endl;
        return 0;
    }
    if(many){
        cout<<0<<endl;
        return 0;
    }
    for(int i=1;i<=n;i++) printf("x%d=%.2f
",i,fabs(a[i][n+1])<eps?0:a[i][n+1]);
    return 0;
}

Luogu P4035球形空间产生器

题目

#include<iostream>
#include<cmath> 
#include<cstdio>
#include<cstdlib>
#define maxn 20
#define ll long long
using namespace std;
template<typename T>
inline void read(T &x){
    x=0;bool flag=0;char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') flag=1;
    for(;isdigit(c);c=getchar()) x=x*10+(c^48);
    if(flag) x=-x;
}

const int eps=1e-8;
int n;
double a[maxn][maxn],b[maxn][maxn];
bool no,many;

void gauss(){
    for(int i=1;i<=n;i++){
        int maxi=i;
        for(int j=i+1;j<=n;j++){
            if(fabs(a[j][i])>fabs(a[maxi][i])) maxi=j;
        }
        for(int j=1;j<=n+1;j++) swap(a[maxi][j],a[i][j]);
        if(fabs(a[i][i])<eps) continue;
        double tmp=a[i][i];//后面a[i][i]会变的，所以不能直接除以a[i][i]，必须拿tmp存一下 
        for(int j=1;j<=n+1;j++) a[i][j]/=tmp;
        for(int j=1;j<=n;j++){
            if(i==j) continue;
            double tmp=a[j][i];
            for(int k=1;k<=n+1;k++) a[j][k]-=a[i][k]*tmp;
        }
    }
    for(int i=1;i<=n;i++){
        int cnt=0;
        for(int j=1;j<=n+1;j++) if(fabs(a[i][j])<eps) cnt++;
        if(cnt==n&&fabs(a[i][n+1])>=eps) no=1;
        if(cnt==n+1) many=1;
        
    }
	if(no||many) return ;
    for(int i=n-1;i;i--)
        for(int j=i+1;j<=n;j++)
            a[i][n+1]-=a[j][n+1]*a[i][j];
}

int main(){
	read(n);
	for(int i=1;i<=n+1;i++)
		for(int j=1;j<=n;j++)
			cin>>b[i][j];
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++) {
			a[i][j]=2*(b[i][j]-b[i+1][j]);
			a[i][n+1]+=b[i][j]*b[i][j]-b[i+1][j]*b[i+1][j];
		}
	}
	gauss();
	for(int i=1;i<=n;i++) printf("%.3f ",fabs(a[i][n+1])<eps?0:a[i][n+1]);
	printf("
");
	return 0;
}

Luogu P4111小Z的房间

题目
矩阵树定理，关于这个的blog先咕着吧~

#include<iostream>
#include<cstdio>
#include<cstdlib>
#define maxn 110
#define int long long
using namespace std;
template<typename T>
inline void read(T &x){
    x=0;bool flag=0;char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') flag=1;
    for(;isdigit(c);c=getchar()) x=x*10+(c^48);
    if(flag) x=-x;
}

const int mod=1e9;
int n,m,b[maxn][maxn],mp[maxn][maxn],cnt,ans=1;
char a[maxn][maxn];

void add(int x,int y){
	mp[x][x]++;mp[y][y]++;mp[x][y]--;mp[y][x]--;
}

void gauss(){
	for(int i=1;i<=cnt;i++){
		for(int j=i+1;j<=cnt;j++){
			while(mp[j][i]){
				int tmp=mp[i][i]/mp[j][i];
				for(int k=i;k<=cnt;k++) mp[i][k]=(mp[i][k]-(mp[j][k]*tmp+mod)%mod+mod)%mod;
				for(int k=1;k<=cnt;k++) swap(mp[i][k],mp[j][k]);
				ans*=-1;
			}
		}
		(ans*=mp[i][i])%=mod;
	}
}

signed main(){
	read(n),read(m);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=m;j++){
			cin>>a[i][j];
			if(a[i][j]=='.') b[i][j]=++cnt;
		}
	cnt--;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=m;j++){
			if(a[i][j]=='.'){
				if(a[i+1][j]=='.') add(b[i][j],b[i+1][j]);
				if(a[i][j+1]=='.') add(b[i][j],b[i][j+1]);
			}
		}		
	gauss();
	cout<<(ans+mod)%mod<<endl;
	return 0;
}

Luogu P3317重建

题目
仍然事矩阵树定理，咕~

#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#define maxn 210
using namespace std;
template<typename T>
inline void read(T &x){
    x=0;bool flag=0;char c=getchar();
    for(;!isdigit(c);c=getchar()) if(c=='-') flag=1;
    for(;isdigit(c);c=getchar()) x=x*10+(c^48);
    if(flag) x=-x;
}
//
//double eps=1e-6;
//int n,cnt=100;
//double g[maxn][maxn],mp[maxn][maxn],ans=1;

//void pre(int x,int y,double p){
//	mp[x][y]-=p/(1-p);
//	mp[x][x]+=p/(1-p);
//}

//void gauss(){
//	for(int i=1;i<=cnt;i++){
//		int maxi=i;
//		for(int j=i+1;j<=cnt;j++){
//			if(fabs(mp[j][i])>fabs(mp[maxi][i])) maxi=j;
//		}
//		if(maxi!=i){
//			for(int j=1;j<=cnt;j++) swap(mp[maxi][j],mp[i][j]);
//			ans*=-1; 
//		}
////		for(int j=i+1;j<=cnt;j++){
////			if(fabs(mp[j][i])>fabs(mp[i][i])){
////				for(int k=1;k<=cnt;k++)  swap(mp[i][k],mp[j][k]);
////				ans*=-1;
////			}
////		}
//		if(fabs(mp[i][i]<eps)){//
//			ans=0;return ;//
//		}//
//		for(int k=i+1;k<=cnt;k++){
//			double t=mp[k][i]/mp[i][i];
//			for(int j=i;j<=cnt;j++) mp[k][j]-=mp[i][j]*t;
//		}
//		ans*=mp[i][i];
//	}
////	for(int i=1;i<=cnt;i++) ans*=mp[i][i];
//}

//int main(){
//	read(n);
//	for(int i=1;i<=n;i++){
//		for(int j=1;j<=n;j++){
//			cin>>g[i][j];
//			g[i][j]+=eps;
//			pre(i,j,g[i][j]);
//			if(i<j) ans*=(1-g[i][j]);
//		}	
//	}
////	cout<<"*"<<ans<<endl;
//	cnt=n-1;
//	gauss();
//	printf("%.3Lf
",fabs(ans));
//	return 0;
//}

double eps=1e-6;
long double z[111][111];
int nh=100;
double f=1;

double gauss(){
	for(int i=1;i<=nh;i++){
		for(int j=i;j<=nh;j++){
			if(fabs(z[j][i])>fabs(z[i][i])){
				f=-f;
				for(int k=1;k<=nh;k++)swap(z[j][k],z[i][k]);
			}
		}
		if(fabs(z[i][i])<=eps)return 0;
		for(int j=i+1;j<=nh;j++)	{
			double x=z[j][i]/z[i][i];
			for(int k=1;k<=nh;k++){
				z[j][k]-=z[i][k]*x;
			}
		}
	}
	for(int i=1;i<=nh;i++){
		f*=z[i][i];
	}
	return f;
}

void qadx(int f,int t,long double pb){
    z[f][t]-=pb/(1-pb);
    z[f][f]+=pb/(1-pb);
}

signed main(){
	int n;
	cin>>n;
	long double tmp;
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			cin>>tmp;
			tmp+=eps;
			qadx(i,j,tmp);
			if(i<j)f*=(1-tmp);
		}
	}
	nh=n-1;
	cout<<gauss()<<endl;
	return 0;
}