bzoj1013 [JSOI2008]球形空间产生器sphere

1013: [JSOI2008]球形空间产生器sphere

Time Limit: 1 Sec  Memory Limit: 162 MB
Submit: 2981  Solved: 1565
[Submit][Status][Discuss]

Description

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

Input

第一行是一个整数，n。接下来的n+1行，每行有n个实数，表示球面上一点的n维坐标。每一个实数精确到小数点后6位，且其绝对值都不超过20000。

Output

有且只有一行，依次给出球心的n维坐标（n个实数），两个实数之间用一个空格隔开。每个实数精确到小数点后3位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。

Sample Input

2
0.0 0.0
-1.0 1.0
1.0 0.0

Sample Output

0.500 1.500

HINT

 

数据规模：

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

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

提示：给出两个定义：

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

2、 距离：设两个n为空间上的点A, B的坐标为(a1, a2, …, an), (b1, b2, …, bn)，则AB的距离定义为：dist = sqrt( (a1-b1)^2 + (a2-b2)^2 + … + (an-bn)^2 )

 

Source

 

题意：给出n维球的球面上的n+1个点，求出球心

分析：显然，我们会想到用模拟退火，但是题目说要与答案一样才得分，而模拟退火什么的很难保证绝对精确，

所以此题只能用高斯消元

根据提示及球心的定义，设球心为（x1, x2, x3, ....,  xn）：

得n+1个式子

(a1-x1)^2+(a2-x2)^2+......+(an-xn)^2 = r^2
(b1-x1)^2+(b2-x2)^2+......+(bn-xn)^2 = r^2
(c1-x1)^2+(c2-x2)^2+......+(cn-xn)^2 = r^2

.........

变形即得n+1个式子

a1^2+a2^2+.....+an^2+   x1^2+x2^2+.....+xn^2   =   (2*a1*x1+2*a2*x2+.......+2*an*xn) 
b1^2+b2^2+.....+bn^2+  x1^2+x2^2+.....+xn^2   =   (2*b1*x1+2*b2*x2+.......+2*bn*xn) 
c1^2+c2^2+.....+cn^2+  x1^2+x2^2+.....+xn^2   =   (2*c1*x1+2*c2*x2+.......+2*cn*xn) 

........

每个式子都与第一个相减，得到n个式子

2*(a1-b1)*x1+2*(a2-b2)*x2+......+2*(an-bn)*xn = (a1^2+a2^2+....+an^2) - (b1^2+b2^2+......+bn^2)
2*(a1-c1)*x1+2*(a2-c2)*x2+......+2*(an-cn)*xn = (a1^2+a2^2+....+an^2) - (c1^2+c2^2+......+cn^2)

..........

其中a1、a2、a3、.......、an，b1、b2、........、bn，.......，都是已知的

所以共有n个式子，n个未知数，用高斯消元即可

 

综上所述，本题得解

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <deque>
#include <vector>
#include <queue>
#include <iostream>
#include <algorithm>
#include <map>
#include <set>
#include <ctime>
using namespace std;
typedef long long LL;
typedef double DB;
#define For(i, s, t) for(int i = (s); i <= (t); i++)
#define Ford(i, s, t) for(int i = (s); i >= (t); i--)
#define MIT (2147483647)
#define INF (1000000001)
#define MLL (1000000000000000001LL)
#define sz(x) ((int) (x).size())
#define clr(x, y) memset(x, y, sizeof(x))
#define puf push_front
#define pub push_back
#define pof pop_front
#define pob pop_back
#define ft first
#define sd second
#define mk make_pair
inline void SetIO(string Name) {
    string Input = Name+".in",
    Output = Name+".out";
    freopen(Input.c_str(), "r", stdin),
    freopen(Output.c_str(), "w", stdout);
}

const int N = 20;
int n;
DB Dat[N][N];
DB K[N][N], Ans[N];

inline void Input() {
    scanf("%d", &n);
    For(i, 0, n)
        For(j, 1, n) scanf("%lf", &Dat[i][j]);
}

inline DB Sqr(DB x) {
    return x*x;
}

inline void Solve() {
    For(i, 1, n) {
        For(j, 1, n)
            K[i][j] = 2.0*(Dat[0][j]-Dat[i][j]);
        K[i][n+1] = 0.0;
        For(j, 1, n)
            K[i][n+1] += (Sqr(Dat[0][j])-Sqr(Dat[i][j]));
    }
    
    For(i, 1, n) {
        DB Max = -1.0*MLL;
        int x;
        For(j, i, n)
            if(Max < K[j][i]) {
                Max = K[j][i];
                x = j;
            }
        For(j, i, n+1) swap(K[i][j], K[x][j]);
        For(j, i+1, n) {
            DB Temp = K[j][i]/K[i][i];
            For(w, i, n+1) K[j][w] -= Temp*K[i][w];
        }
    }
    
    Ford(i, n, 1) {
        Ford(j, n, i+1) K[i][n+1] -= K[i][j]*Ans[j];
        Ans[i] = K[i][n+1]/K[i][i];
    }
    
    For(i, 1, n-1) printf("%.3lf ", Ans[i]);
    printf("%.3lf\n", Ans[n]);
}

int main() {
    SetIO("1013");
    Input();
    Solve();
    return 0;
}