板子
#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;
}