球形空间产生器

球形空间产生器

洛谷&天立OI

思路

奇妙的好题啊！！

推导

首先我们学习孟德尔的遗传定律发现的经验，从一般到特殊。

假设 \(n=2\)，三个点的坐标 ：$ (a_1,b_1),(a_2,b_2),(a_3,b_3)$ 。

设球心为\((x_1,x_2)\)，半径为 \(R\)。

易得到以下式子：

\[\begin{align}\notag &(x_1-a_1)^2+(x_2-b_2)^2=R^2 \iff x_1^2+x_2^2-R^2+a_1^2+b_1^2=2a_1x_1+2b_1x_2 \end{align} \]

所以有：

\[\begin{align} \tag1 x_1^2+x_2^2-R^2+a_1^2+b_1^2=2a_1x_1+2b_1x_2\\ \tag2 x_1^2+x_2^2-R^2+a_2^2+b_2^2=2a_2x_1+2b_2x_2\\ \tag3 x_1^2+x_2^2-R^2+a_3^2+b_3^2=2a_3x_1+2b_3x_2 \end{align} \]

可见这是一个关于\(x\)的二元二次方程组，但是高斯消元是一次的怎么办？

不妨可以这样：

\[令 c_n=a_n^2+b_n^2\\ \begin{align} \tag4 (1)-(3) \iff c_1-c_3=(2a_1-2a_3)x_1+(2b_1-2b_3)x_2\\ \tag5 (2)-(3) \iff c_2-c_3=(2a_2-2a_3)x_1+(2b_2-2b_3)x_2 \end{align} \]

这样我们就可以得到一个二元一次方程组：

\[\begin{cases} \notag c_1-c_3=(2a_1-2a_3)x_1+(2b_1-2b_3)x_2\\ c_2-c_3=(2a_2-2a_3)x_1+(2b_2-2b_3)x_2 \end{cases} \]

接下来就简单啦！直接高斯约旦消元就做出来了！！

但如果\(n!=2\)呢？

其实也很简单，我们稍微想一下也可以得到一个\(n\)元\(n\)次方程组：

\[\notag 令 c_n=a_n^2+b_n^2+c_n^2+......\\ 令 圆心坐标为(x_1,x_2,x_3,x_4,,,,,x_n)\\ 令 给定球面点坐标为(a_n,b_n,c_n,,,,A_n)\\ \begin{cases} c_1-c_{n+1}=(2a_1-2a_{n+1})x_1+(2b_1-2b_{n+1})x_2+(2c_1-2c_{n+1})x_3+.....+(2A_1-2A_{n+1})x_n\\ c_2-c_{n+1}=(2a_2-2a_{n+1})x_1+(2b_2-2b_{n+1})x_2+(2c_2-2c_{n+1})x_3+.....+(2A_2-2A_{n+1})x_n\\ .\\ .\\ .\\ c_n-c_{n+1}=(2a_n-2a_{n+1})x_1+(2b_n-2b_{n+1})x_2+(2c_n-2c_{n+1})x_3+.....+(2A_n-2A_{n+1})x_n\\ \end{cases} \]

这这就是一个高斯消元即可！！！！

CODE

#include<bits/stdc++.h>
#define ll long long
#define ld long double
using namespace std;
const ll maxn = 10+3;
ll ll_maxn=1e18;
inline ll read_int(){
    ll a=0,f=0,g=getchar();
    while(g<'0'||g>'9'){if(g=='-') f=1;g=getchar();}
    while('0'<=g&&g<='9') a=a*10+g-'0',g=getchar();
    return f ? -a : a;
}

inline void write(ll s,bool f){
    ll top=0,a[40];
    if(s<0) s=-s,putchar('-');
    while(s) a[++top]=s%10,s/=10;
    if(top==0) a[++top]=0;
    while(top) putchar(a[top]+'0'),top--;
    if(f) putchar('\n');
}

int n;
ld a[maxn][maxn];
ld b[maxn][maxn];
ld eps=1e-6;

inline void gaosi(){
	for(int i=1;i<=n;i++){
		int max_set=i;
		for(int e=i+1;e<=n;e++){
			if(fabs(a[e][i])>fabs(a[max_set][i])) max_set=e;
		}
		if(max_set^i) for(int e=i;e<=n+1;e++) swap(a[max_set][e],a[i][e]);
		for(int e=n+1;e>=i;e--) a[i][e]/=a[i][i];
		for(int e=1;e<=n;e++){
			if(e==i) continue;
			for(int k=n+1;k>=i;k--){
				a[e][k]-=a[e][i]*a[i][k];
			}
		}
	}
}

inline void PRINTF(){
	for(int i=1;i<=n;i++) printf("%.3Lf ",a[i][n+1]);
}

inline void read(){
	n=read_int();
	for(int i=1;i<=n+1;i++){
		for(int e=1;e<=n;e++){
			scanf("%Lf",&b[i][e]);
			b[i][n+1]+=b[i][e]*b[i][e];
			b[i][e]*=2;
		}
	}
	for(int i=1;i<=n;i++){
		for(int e=1;e<=n+1;e++){
			a[i][e]=b[i][e]-b[n+1][e];
		}
	}
	gaosi();
	PRINTF();
}

int main (){
	read();
}