[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<<"
";
	}
}
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;
}