球
型
空
间
产
生
器
球型空间产生器
球型空间产生器
D
e
s
c
r
i
p
t
i
o
n
\mathcal{Description}
Description
有一个球形空间产生器能够在
n
n
n 维空间中产生一个坚硬的球体。现在,你被困在了这个 n 维球体中,你只知道球面上
n
+
1
n+1
n+1 个点的坐标,你需要以最快的速度确定这个
n
n
n 维球体的球心坐标,以便于摧毁这个球形空间产生器。
n
<
=
10
n <= 10
n<=10
S
o
l
u
t
i
o
n
\mathcal{Solution}
Solution
题意: 求 n n n维点 O O O( a 1 , a 2 , . . . , a n a_1,a_2,...,a_n a1,a2,...,an), 使 O O O 距离 球面上 n + 1 n+1 n+1 个点的 距离相等,
使用 式子 表示出来, 即
(
a
1
−
b
1
,
1
)
2
+
(
a
2
−
b
1
,
2
)
2
+
.
.
.
+
(
a
n
−
b
1
,
n
)
2
=
D
2
(
a
1
−
b
2
,
1
)
2
+
(
a
2
−
b
2
,
2
)
2
+
.
.
.
+
(
a
n
−
b
2
,
n
)
2
=
D
2
.
.
.
(
a
1
−
b
n
+
1
,
1
)
2
+
(
a
2
−
b
n
+
1
,
2
)
2
+
.
.
.
+
(
a
n
−
b
n
+
1
,
n
)
2
=
D
2
\ \\ (a_1-b_{1,1})^2+(a_2-b_{1,2})^2+...+(a_n-b_{1,n})^2=D^2\\ \ \\ (a_1-b_{2,1})^2+(a_2-b_{2,2})^2+...+(a_n-b_{2,n})^2=D^2\\ . \\ . \\ . \\ \ \\ (a_1-b_{n+1,1})^2+(a_2-b_{n+1,2})^2+...+(a_n-b_{n+1,n})^2=D^2 \\ \ \\
(a1−b1,1)2+(a2−b1,2)2+...+(an−b1,n)2=D2 (a1−b2,1)2+(a2−b2,2)2+...+(an−b2,n)2=D2... (a1−bn+1,1)2+(a2−bn+1,2)2+...+(an−bn+1,n)2=D2
相邻式子相减消去
D
2
D^2
D2, 得到
(
2
a
1
−
b
1
,
1
−
b
2
,
1
)
(
b
2
,
1
−
b
1
,
1
)
+
.
.
.
+
(
2
a
n
−
b
1
,
n
−
b
2
,
n
)
(
b
2
,
n
−
b
1
,
n
)
=
0
.
.
略
\ \\ (2a_1-b_{1,1}-b_{2,1})(b_{2,1}-b_{1,1})+...+(2a_n-b_{1,n}-b_{2,n})(b_{2,n}-b_{1,n})=0\\ . \\ . \\ 略 \\ \ \\
(2a1−b1,1−b2,1)(b2,1−b1,1)+...+(2an−b1,n−b2,n)(b2,n−b1,n)=0..略
发现已经是一个 线性方程组 了,
系数 A [ i , j ] = 2 ∗ ( b [ i + 1 , j ] − b [ i , j ] ) A[i,j]=2*(b[i+1,j]-b[i,j]) A[i,j]=2∗(b[i+1,j]−b[i,j])
常数 A [ i , n + 1 ] = ∑ j = 1 n ( b [ i , j ] + b [ i + 1 , j ] ) ∗ ( b [ i + 1 , j ] − b [ i , j ] ) A[i,n+1]=\sum_{j=1}^n \ (b[i,j] + b[i+1,j])\ *\ (b[i+1,j]-b[i,j]) A[i,n+1]=j=1∑n (b[i,j]+b[i+1,j]) ∗ (b[i+1,j]−b[i,j])
于是 高斯消元 O ( N 3 ) O(N^3) O(N3) 求解 .
C
o
d
e
\mathcal{Code}
Code
#include<bits/stdc++.h>
#define reg register
const int maxn = 15;
int N;
double A[maxn][maxn];
double B[maxn][maxn];
int main(){
scanf("%d", &N);
for(reg int i = 1; i <= N+1; i ++)
for(reg int j = 1; j <= N; j ++) scanf("%lf", &B[i][j]);
for(reg int i = 1; i <= N; i ++){
for(reg int j = 1; j <= N; j ++) A[i][j] = 2.0*(B[i+1][j] - B[i][j]);
A[i][N+1] = 0;
for(reg int j = 1; j <= N; j ++)
A[i][N+1] += (B[i][j]+B[i+1][j])*(B[i+1][j] - B[i][j]);
}
for(reg int i = 1; i <= N; i ++){
int max_id = i;
for(reg int j = i+1; j <= N; j ++)
if(fabs(A[max_id][i]) < fabs(A[j][i])) max_id = j;
std::swap(A[i], A[max_id]);
double tmp = A[i][i];
for(reg int j = i; j <= N+1; j ++) A[i][j] /= tmp;
for(reg int j = i+1; j <= N; j ++){
tmp = A[j][i];
for(reg int k = i; k <= N+1; k ++)
A[j][k] -= A[i][k] * tmp;
}
}
for(reg int i = N; i >= 1; i --)
for(reg int j = i+1; j <= N; j ++) A[i][N+1] -= A[i][j]*A[j][N+1];
for(reg int i = 1; i <= N; i ++) printf("%.3lf ", A[i][N+1]);
return 0;
}