bzoj1013/luogu4035 球形空间生成器 (高斯消元)

时间:2023-03-08 23:27:36
bzoj1013/luogu4035 球形空间生成器 (高斯消元)

根据各点到圆心的距离相等,可以列出有N个等号的方程

假设圆心坐标是(x,y,z,...)的话,$x^2,y^2,z^2$等是可以消掉的

于是整理一下,就变成了N元1次方程组,有N个方程、而且保证是相容的

高斯消元的话,就是拿着第一式去把剩下的第一项都消了,再拿第二式把剩下的第二项都消了,...到最后方程组呈一个阶梯状,然后从最后一点点往回带就能解出来

复杂度$O(n^3)$

 #include<bits/stdc++.h>
#define pa pair<int,int>
#define CLR(a,x) memset(a,x,sizeof(a))
using namespace std;
typedef long long ll;
const int maxn=; inline ll rd(){
ll x=;char c=getchar();int neg=;
while(c<''||c>''){if(c=='-') neg=-;c=getchar();}
while(c>=''&&c<='') x=x*+c-'',c=getchar();
return x*neg;
} int N;
double pos[maxn][maxn],a[maxn][maxn],b[maxn]; int main(){
//freopen("","r",stdin);
int i,j,k;
N=rd();
for(i=;i<=N+;i++){
for(j=;j<=N;j++){
scanf("%lf",&pos[i][j]);
}
if(i>){
for(j=;j<=N;j++){
a[i-][j]=*(pos[i][j]-pos[i-][j]);
a[i-][N+]+=pos[i][j]*pos[i][j]-pos[i-][j]*pos[i-][j];
}
}
}
for(i=;i<=N;i++){
for(j=i+;j<=N;j++){
double p=a[j][i]/a[i][i];
for(k=i;k<=N+;k++){
a[j][k]-=p*a[i][k];
}
}
}
for(i=N;i;i--){
double r=a[i][N+];
for(j=N;j>i;j--){
r-=b[j]*a[i][j];
}
b[i]=r/a[i][i];
}
for(i=;i<=N;i++)
printf("%.3lf ",b[i]);
return ;
}