球形空间产生器

思路

奇妙的好题啊！！

推导

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

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

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

易得到以下式子：

\begin{aligned} (x_{1} - a_{1})^{2} + (x_{2} - b_{2})^{2} = R^{2} ⟺ x_{1}^{2} + x_{2}^{2} - R^{2} + a_{1}^{2} + b_{1}^{2} = 2 a_{1} x_{1} + 2 b_{1} x_{2} \end{aligned}

$\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{aligned} (1) & x_{1}^{2} + x_{2}^{2} - R^{2} + a_{1}^{2} + b_{1}^{2} = 2 a_{1} x_{1} + 2 b_{1} x_{2} \\ (2) & x_{1}^{2} + x_{2}^{2} - R^{2} + a_{2}^{2} + b_{2}^{2} = 2 a_{2} x_{1} + 2 b_{2} x_{2} \\ (3) & x_{1}^{2} + x_{2}^{2} - R^{2} + a_{3}^{2} + b_{3}^{2} = 2 a_{3} x_{1} + 2 b_{3} x_{2} \end{aligned}

$\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{aligned} (4) & (1) - (3) ⟺ c_{1} - c_{3} = (2 a_{1} - 2 a_{3}) x_{1} + (2 b_{1} - 2 b_{3}) x_{2} \\ (5) & (2) - (3) ⟺ c_{2} - c_{3} = (2 a_{2} - 2 a_{3}) x_{1} + (2 b_{2} - 2 b_{3}) x_{2} \end{aligned}

$令 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} c_{1} - c_{3} = (2 a_{1} - 2 a_{3}) x_{1} + (2 b_{1} - 2 b_{3}) x_{2} \\ c_{2} - c_{3} = (2 a_{2} - 2 a_{3}) x_{1} + (2 b_{2} - 2 b_{3}) x_{2} \end{cases}

$\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$ 次方程组：

令 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} = (2 a_{1} - 2 a_{n + 1}) x_{1} + (2 b_{1} - 2 b_{n + 1}) x_{2} + (2 c_{1} - 2 c_{n + 1}) x_{3} + . . . . . + (2 A_{1} - 2 A_{n + 1}) x_{n} \\ c_{2} - c_{n + 1} = (2 a_{2} - 2 a_{n + 1}) x_{1} + (2 b_{2} - 2 b_{n + 1}) x_{2} + (2 c_{2} - 2 c_{n + 1}) x_{3} + . . . . . + (2 A_{2} - 2 A_{n + 1}) x_{n} \\ . \\ . \\ . \\ c_{n} - c_{n + 1} = (2 a_{n} - 2 a_{n + 1}) x_{1} + (2 b_{n} - 2 b_{n + 1}) x_{2} + (2 c_{n} - 2 c_{n + 1}) x_{3} + . . . . . + (2 A_{n} - 2 A_{n + 1}) x_{n} \end{cases}

$\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();
}