[2018.12.18]BZOJ1013 [JSOI2008]球形空间产生器sphere

设中心点坐标\(B=(b_1,b_2,b_3,...,b_n)\),球面上的点的坐标A=\((a_1,a_2,a_3,...,a_n)\)

则

\(dist_{A,B}=\sqrt{(a_1-b_1)^2+(a_2-b_2)^2+(a_3-b_3)^2+...+(a_n-b_n)^2}\)

\((a_1-b_1)^2+(a_2-b_2)^2+(a_3-b_3)^2+...+(a_n-b_n)^2=dist_{A,B}^2\)

\(a_1^2-2a_1b_1+b_1^2+a_2^2-2a_2b_2+b_2^2+a_3^2-2a_3b_3+b_3^2+...+a_n^2-2a_nb_n+b_n^2=dist_{A,B}^2\)

\(-2a_1b_1-2a_2b_2-2a_3b_3-...-2a_nb_n=dist_{A,B}^2-(a_1^2+a_2^2+a_3^2+...+a_n^2)-(b_1^2+b_2^2+b_3^2+...+b_n^2)\)

设\(dist_{A,B}^2-(b_1^2+b_2^2+b_3^2+...+b_n^2)=k\)

得\(2a_1b_1+2a_2b_2+2a_3b_3+...+2a_nb_n+k=a_1^2+a_2^2+a_3^2+...+a_n^2\)

由于\(a_1,a_2,a_3,...,a_n\)都是常数,所以这实际上是一个\(n+1\)元一次方程。

我们有\(n+1\)个这样的方程。

比如样例,可得以下方程:

\(0b_1+0b_2+k=0\)

\(-2b_1+2b_2+k=2\)

\(2b_1+0b_2+k=1\)

解得\(b_1=\frac{1}{2},b_2=\frac{3}{2},k=0\)

至于算法,高斯消元求解即可。

code:

#include<bits/stdc++.h>
using namespace std;
struct equation{
	double n[15],x;
}e[15];
int n;
double ans[15];
void Printeq(){//debug
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++)cout<<e[i].n[j]<<" ";
		cout<<"=>"<<e[i].x<<"\n";
	}
}
void Work(int x,int y){
	double v=-e[y].n[x]/e[x].n[x];
	for(int i=1;i<=n;i++)e[y].n[i]+=e[x].n[i]*v;
	e[y].x+=e[x].x*v;
}
void Solve(){
	for(int i=1;i<n;i++){
		for(int j=i;j<=n;j++){
			if(e[j].n[i]){
				swap(e[j],e[i]);
				break;
			}
		}
		for(int j=i+1;j<=n;j++)Work(i,j);
	}
	ans[n]=e[n].x/e[n].n[n];
	double k=0;
	for(int i=n-1;i>=1;i--,k=0){
		for(int j=i+1;j<=n;j++)k+=ans[j]*e[i].n[j];
		ans[i]=(e[i].x-k)/e[i].n[i];
	}
}
int main(){
	scanf("%d",&n);
	n++;
	for(int i=1;i<=n;i++){
		for(int j=1;j<n;j++)scanf("%lf",&e[i].n[j]),e[i].x+=e[i].n[j]*e[i].n[j],e[i].n[j]*=2;
		e[i].n[n]=1.0;
	}
	Solve();
	for(int i=1;i<n;i++)printf("%.3lf ",ans[i]);
	return 0;
}