P4035 [JSOI2008]球形空间产生器 （向量，高斯消元）
题面
有一个 n n n 维球，给定 n + 1 n+1 n+1 个在球面上的点，求球心坐标。
n ≤ 10 n\leq 10 n≤10 。
题解
好久以前的题了，昨天首 A 。
n n n 太小了！明明可以开 100！
看到题解里有两种做法：
- 模拟退火
- 距离半径法高斯消元
后面那种，以前一直搞不懂，不知道怎么解二次的。于是很久都没做出来。
现在，我有一种更好的方法：向量高斯消元法。
我们知道，对于球上两个点 A , B A,B A,B 和球心 P P P ，满足： P P P 到 A B AB AB 中点的连线与 A B AB AB 垂直。
于是，我们就可以列出这样的等式：
A
B
⟶
⋅
(
A
P
⟶
+
B
P
⟶
)
=
0
\stackrel{\longrightarrow}{AB}\cdot(\stackrel{\longrightarrow}{AP}+\stackrel{\longrightarrow}{BP})=0\\
AB⟶⋅(AP⟶+BP⟶)=0
即
(
O
B
⟶
−
O
A
⟶
)
⋅
(
2
O
P
⟶
−
(
O
B
⟶
+
O
A
⟶
)
)
=
0
(
O
B
⟶
−
O
A
⟶
)
⋅
2
O
P
⟶
=
(
O
B
⟶
−
O
A
⟶
)
⋅
(
O
B
⟶
+
O
A
⟶
)
=
O
B
⟶
2
−
O
A
⟶
2
(\stackrel{\longrightarrow}{OB}-\stackrel{\longrightarrow}{OA})\cdot(2\stackrel{\longrightarrow}{OP}-(\stackrel{\longrightarrow}{OB}+\stackrel{\longrightarrow}{OA}))=0\\ (\stackrel{\longrightarrow}{OB}-\stackrel{\longrightarrow}{OA})\cdot2\stackrel{\longrightarrow}{OP}=(\stackrel{\longrightarrow}{OB}-\stackrel{\longrightarrow}{OA})\cdot(\stackrel{\longrightarrow}{OB}+\stackrel{\longrightarrow}{OA})=\stackrel{\longrightarrow}{OB}^2-\stackrel{\longrightarrow}{OA}^2
(OB⟶−OA⟶)⋅(2OP⟶−(OB⟶+OA⟶))=0(OB⟶−OA⟶)⋅2OP⟶=(OB⟶−OA⟶)⋅(OB⟶+OA⟶)=OB⟶2−OA⟶2
从坐标来看，就是
∑
i
=
1
n
2
(
X
B
[
i
]
−
X
A
[
i
]
)
⋅
X
P
[
i
]
‾
=
∣
O
B
∣
2
−
∣
O
A
∣
2
\sum_{i=1}^n2(X_B[i]-X_A[i])\cdot \underline{~X_P[i]~}=|OB|^2-|OA|^2
i=1∑n2(XB[i]−XA[i])⋅ XP[i] =∣OB∣2−∣OA∣2
其中画横线的就是方程未知量，右边是两个距离的平方相减。
我们把 n + 1 n+1 n+1 个点中每两个相邻的点列一个这样的方程。
下一步拉板子。
CODE
一定有解，所以不考虑这么多了。
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<bitset>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define MAXN 205
#define LL long long
#define DB double
#define lowbit(x) (-(x) & (x))
#define ENDL putchar('\n')
#define FI first
#define SE second
#define SQ 450
#define eps 1e-9
LL read() {
LL f=1,x=0;int s = getchar(); if(s<0) return -1;
while(s < '0' || s > '9') {if(s=='-')f = -f;s = getchar();}
while(s >= '0' && s <= '9') {x=x*10+(s^48);s = getchar();}
return f*x;
}
void putpos(LL x) {if(!x)return ;putpos(x/10);putchar('0'+(x%10));}
void putnum(LL x) {
if(!x) {putchar('0');return ;}
if(x<0) {putchar('-');x = -x;}
return putpos(x);
}
void AIput(LL x,int c) {putnum(x);putchar(c);}
int MOD = 1;
int n,m,s,o,k;
DB pt[105][105];
DB a[105][105];
DB Abs(DB x) {return x<0 ? -x:x;}
void Gauss(int n) {
for(int i = 1;i <= n;i ++) {
for(int j = i+1;j <= n;j ++) {
if(Abs(a[j][i]) > Abs(a[i][i])) {
swap(a[j],a[i]); break;
}
}
for(int j = i+1;j <= n;j ++) {
DB nm = a[j][i] / a[i][i];
for(int k = n+1;k >= i;k --) {
a[j][k] -= a[i][k] * nm;
}
}
}
for(int i = n;i > 0;i --) {
for(int j = n;j > i;j --) {
a[i][n+1] -= a[i][j] * a[j][n+1];
}
a[i][n+1] /= a[i][i];
}return ;
}
int main() {
n = read();
for(int i = 1;i <= n+1;i ++) {
for(int j = 1;j <= n;j ++) {
scanf("%lf",&pt[i][j]);
}
}
for(int i = 1;i <= n;i ++) {
for(int j = 1;j <= n;j ++) {
a[i][j] = 2.0*(pt[i+1][j]-pt[i][j]);
a[i][n+1] += pt[i+1][j]*pt[i+1][j] - pt[i][j]*pt[i][j];
}
}
Gauss(n);
for(int i = 1;i <= n;i ++) {
printf("%.3f",a[i][n+1]);
putchar(i==n ? '\n':' ');
}
return 0;
}
