BZOJ P1013 [JSOI2008]球形空间产生器

发表于 | 更新于 | 阅读次数：

题目描述

有一个球形空间产生器能够在n维空间中产生一个坚硬的球体。现在，你被困在了这个n维球体中，你只知道球面上n+1个点的坐标，你需要以最快的速度确定这个n维球体的球心坐标，以便于摧毁这个球形空间产生器。

输入格式

第一行是一个整数n。

接下来的n+1行，每行有n个实数，表示球面上一点的n维坐标。

输出格式

有且只有一行，依次给出球心的n维坐标（n个实数），两个实数之间用一个空格隔开。每个实数精确到小数点后3位。

数据保证有解。你的答案必须和标准输出一模一样才能够得分。

样例输入

1
2
3
4
2
0.0 0.0
-1.0 1.0
1.0 0.0

样例输出

1
0.500 1.500

数据规模与约定

对于40%的数据，$1<=n<=3$。

对于100%的数据，$1<=n<=10$。

1、 球心：到球面上任意一点距离都相等的点。

2、 距离：设两个n为空间上的点A, B的坐标为$(a_1,a_2,…,a_n),(b_1,b_2,…,b_n)，则AB的距离定义为：

$dist=sqrt((a_1−b_1)^2+(a_2−b_2)^2+…+(a_n−b_n)^2)$

数据中所有数据绝对值小于 10000 。

  • 因为是一个球，所以每个点到球心的距离都相等，

我们设这个半径为R，球心坐标为$O(x_1,x_2,….,x_n)$;

那么对于每一个点$P(ai1,ai2,…,ain)$：我们易得

$sqrt ( ( ai1 - x1 ) ^ 2 + ( ai2 - x2 ) ^ 2 + … + ( ain - xn ) ^2 ) = R$;

  • 考虑构造方程组

将上式两侧平方再展开，得

$-2 ( ai1 x1 + ai2 x2 + … + ain * xn )+( ai1 ^ 2 + ai2 ^ 2 + … + ain ^ 2 )+( x1 ^ 2 + x2 ^ 2 + … + xn ^ 2 )$

= $R ^ 2$;

给出n+1个点，n个点就可以构造该方程组，那多给的一个点是用来（设选第一个点为这个点） 消掉重复出现的部分$( x1 ^ 2 + x2 ^ 2 + … + xn ^ 2 )$和$R ^ 2$，

即令其他所有的点的方程都减掉多的一个点的方程，整理得到其他方程格式为：

$2 ( ai1 x1 + ai2 x2 + … + ain * xn )$

=$ ( ai1 ^ 2 + ai2 ^ 2 + … + ain ^ 2 ) - ( a11 ^ 2 + a12 ^ 2 + … + a1n ^ 2 ) - 2 ( a11 x1 + a12 x2 + … + a1n * xn ) $;

右侧是常数，左侧展开就是一个愉快的高斯消元方程组.

Code…

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
#include<bits/stdc++.h>

const double eps=1e-8;
const int maxn=101;
double det;
/*Namespace's template*/
template<class type>type _max(type a, type b){return a > b ? a : b;}
template<class type>type _min(type a, type b){return a < b ? a : b;}
template<class type>type _abs(type a){return a < 0 ? -a : a;}
template<class type>type _gcd(type a, type b){return (!b) ? a : gcd(b, a % b);}
template<class type>type _mod(type a, type p){if(a >= 0 && a < p)return a;a %= p;if(a < 0)a += p;return a;}
template<class type>type _qpow(type a, type b, type p){type ans = 1;for(; b; b >>= 1, a = _mod(a * a, p))if(b & 1)ans = _mod(ans * a, p);return ans;}
char buf[1<<21],*p1=buf,*p2=buf;
inline int getc() {
    return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;
}
inline int r() {
	register int x = 0; register char ch = getchar();
	for(; ch < '0' || ch > '9'; ch = getchar());
	for(; ch >= '0' && ch <= '9'; ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ 48);
	return x;
}
inline double r_db () {
    double s = 0.0; int f = 1; char c = getchar();
    while (c < '0' || c > '9') {if (c == '-') f = -1;c = getchar();}
    while (c >= '0' && c <= '9') {s = s * 10 + (c ^ 48);c = getchar();}
    if (c == '.') {int k = 10;c = getchar();while (c >= '0' && c <= '9') {s += (double)(c ^ 48) / k;k *= 10;c = getchar();}}
    return s;
}
/*end template*/
inline int Gauss(double a[][maxn],bool l[],double ans[],const int &n,const int &m){
    int res=0,r=0;
    for(int i=0;i<n;++i) l[i]=false;
    for(int i=0;i<n;++i){
        for(int j=r;j<m;++j)
            if(fabs(a[j][i])>eps){
                for(int k=i;k<=n;++k)
                    std::swap(a[j][k],a[r][k]);
                break;
            }
        if(fabs(a[r][i])<eps){
            ++res;
            continue;
        }
        for(int j=0;j<m;++j)
            if(j!=r && fabs(a[j][i])>eps){
                double tmp=a[j][i]/a[r][i];
                for(int k=i;k<=n;++k)
                    a[j][k]-=tmp*a[r][k];
            }
        l[i]=true,++r;
    }
    det=1;
    for(int i=0;i<n;++i)
        det*=a[i][i];
    //检查是否无解
    for(int i=n-res;i<m;++i){
        if(fabs(a[i][n])>eps) return -1;
    }
    for(int i=0;i<n;++i){
        if(l[i])//不是自由变元
            for(int j=0;j<n;++j)
                if(fabs(a[j][i])>eps)
                        ans[i]=a[j][n]/a[j][i];
    }
    return res;//返回自由变元数
}

int main() {
    double a[maxn][maxn];
    bool l[maxn];
    double ans[maxn];
    double arr[maxn];
    double brr[maxn];
    int n,m,res;
    memset(a,0,sizeof(a));
    n = r();
    m=n;
    for(int i=0;i<n;++i) scanf("%lf", &arr[i]);
    for(int i=0;i<n;++i){
        for(int j=0;j<n;++j) scanf("%lf", &brr[j]);
        for(int j=0;j<n;++j) a[i][j]=2*(brr[j]-arr[j]);
        for(int j=0;j<n;++j) a[i][n]+=(brr[j]*brr[j]-arr[j]*arr[j]);
        for(int j=0;j<n;++j) arr[j]=brr[j];
    }
    res=Gauss(a,l,ans,n,m);
    for(int i=0;i<n;++i)//
        if(i==n-1) printf("%.3lf\n",ans[i]);
        else printf("%.3lf ",ans[i]);
    return 0;
}