[JSOI2008]球形空间产生器sphere

Description

[JSOI2008]球形空间产生器sphere

Solution

以二维的时候为例
球面上三个点分别为\((a_1,a_2)\)，\((b_1,b_2)\)，\((c_1,c_2)\)
设球心的点为\((x_1,x_2)\)，球的半径为\(T\)
列式得
\(\left\{ \begin{aligned} (a_1-x_1)^2+(a_2-x_2)^2&=(a_1)^2-2a_1x_1+(x_1)^2+(a_2)^2-2a_2x_2+(x_2)^2&=T\\ (b_1-x_1)^2+(b_2-x_2)^2&=(b_1)^2-2b_1x_1+(x_1)^2+(b_2)^2-2b_2x_2+(x_2)^2&=T\\ (c_1-x_1)^2+(c_2-x_2)^2&=(c_1)^2-2c_1x_1+(x_1)^2+(c_2)^2-2c_2x_2+(x_2)^2&=T \end{aligned} \right.\)
联立方程
\((a_1-x_1)^2+(a_2-x_2)^2=(b_1-x_1)^2+(b_2-x_2)^2=(c_1-x_1)^2+(c_2-x_2)^2\)
暴力展开
\((a_1)^2-2a_1x_1+(x_1)^2+(a_2)^2-2a_2x_2+(x_2)^2=(b_1)^2-2b_1x_1+(x_1)^2+(b_2)^2-2b_2x_2+(x_2)^2=(c_1)^2-2c_1x_1+(x_1)^2+(c_2)^2-2c_2x_2+(x_2)^2\)
整理得
\((a_1)^2-2a_1x_1+(a_2)^2-2a_2x_2=(b_1)^2-2b_1x_1+(b_2)^2-2b_2x_2=(c_1)^2-2c_1x_1+(c_2)^2-2c_2x_2\)
也可以写成这个样子
\(\left\{ \begin{aligned} (a_1)^2-2a_1x_1+(a_2)^2-2a_2x_2&=(b_1)^2-2b_1x_1+(b_2)^2-2b_2x_2\\ (b_1)^2-2b_1x_1+(b_2)^2-2b_2x_2&=(c_1)^2-2c_1x_1+(c_2)^2-2c_2x_2 \end{aligned} \right.\)
移来移去
\(\left\{ \begin{aligned} -2a_1x_1-2a_2x_2+2b_1x_1+2b_2x_2&=-(a_1)^2-(a_2)^2+(b_1)^2+(b_2)^2\\ -2b_1x_1-2b_2x_2+2c_1x_1+2c_2x_2&=-(b_1)^2-(b_2)^2+(c_1)^2+(c_2)^2 \end{aligned} \right.\)
移来移去
\(\left\{ \begin{aligned} (2b_1-2a_1)x_1+(2b_2-2a_2)x_2&=-(a_1)^2-(a_2)^2+(b_1)^2+(b_2)^2\\ (2c_1-2b_1)x_1+(2c_2-2b_2)x_2&=-(b_1)^2-(b_2)^2+(c_1)^2+(c_2)^2 \end{aligned} \right.\)
然后我们发现它成了一个\(n\)元一次方程组，于是我们来高斯消元就行了

Code

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define eps 1e-7
#define N 20
int n;
double c[N][N], a[N][N], ans[N];
inline int read() {
	int s = 0, w = 1;
	char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') w = -1;
	for (; isdigit(c); c = getchar()) s = (s << 1) + (s << 3) + (c ^ 48);
	return s * w;
}
int main() {
	n = read();
	for (register int i = 1; i <= n + 1; i++)
		for (register int j = 1; j <= n; j++)
			scanf("%lf", &c[i][j]);
	for (register int i = 1; i <= n; i++)
		for (register int j = 1; j <= n; j++)
			a[i][j] = - 2 * c[i][j] + 2 * c[i + 1][j],
			a[i][n + 1] += - c[i][j] * c[i][j] + c[i + 1][j] * c[i + 1][j];
	for (register int i = 1; i <= n; i++) {
		int r = i;
		for (register int j = i + 1; j <= n; j++)
			if (fabs(a[r][i]) < fabs(a[j][i]))
				r = j;
		if (fabs(a[r][i]) < eps) break;
		if (i != r) swap(a[i], a[r]);
		double div = a[i][i];
		for (register int j = i; j <= n + 1; j++)
			a[i][j] /= div;
		for (register int j = i + 1; j <= n; j++) {
			div = a[j][i];
			for (register int k = i; k <= n + 1; k++)
				a[j][k] -= a[i][k] * div;
		}
	}
	ans[n] = a[n][n + 1];
	for (register int i = n - 1; i; i--) {
		ans[i] = a[i][n + 1];
		for (register int j = i + 1; j <= n; j++)
			ans[i] -= a[i][j] * ans[j];
	}
	for (register int i = 1; i <= n; i++)
		printf("%.3lf ", ans[i]);
	return 0;
}