BZOJ 1013 球形空间产生器sphere 高斯消元

题目链接：

https://www.lydsy.com/JudgeOnline/problem.php?id=1013

题目大意：

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

思路：

每两个可以构成一个n个变量的式子，可以构造出n个不同的式子，进行高斯消元求解。

高斯消元模板：

 1 typedef double Matrix[maxn][maxn];
 2 void gauss(Matrix A, int n)
 3 {
 4     for(int i = 0; i < n; i++)
 5     {
 6         int r = i;
 7         for(int j = i + 1; j < n; j++)
 8             if(fabs(A[j][i]) > fabs(A[r][i]))r = j;
 9         if(r != i)for(int j = 0; j <= n; j++)swap(A[r][j], A[i][j]);
10         for(int k = i + 1; k < n; k++)
11         {
12             double f = A[k][i] / A[i][i];
13             for(int j = i; j <= n; j++)A[k][j] -= f * A[i][j];
14         }
15     }
16  
17     for(int i = n - 1; i >= 0; i--)
18     {
19         for(int j = i + 1; j < n; j++)
20         {
21             A[i][n] -= A[j][n] * A[i][j];
22         }
23         A[i][n] /= A[i][i];
24     }
25 }

 1 #include<bits/stdc++.h>
 2 #define IOS ios::sync_with_stdio(false);//不可再使用scanf printf
 3 #define Max(a, b) ((a) > (b) ? (a) : (b))//禁用于函数，会超时
 4 #define Min(a, b) ((a) < (b) ? (a) : (b))
 5 #define Mem(a) memset(a, 0, sizeof(a))
 6 #define Dis(x, y, x1, y1) ((x - x1) * (x - x1) + (y - y1) * (y - y1))
 7 #define MID(l, r) ((l) + ((r) - (l)) / 2)
 8 #define lson ((o)<<1)
 9 #define rson ((o)<<1|1)
10 #pragma comment(linker, "/STACK:102400000,102400000")//栈外挂
11 using namespace std;
12 inline int read()
13 {
14     int x=0,f=1;char ch=getchar();
15     while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
16     while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
17     return x*f;
18 }
19  
20 typedef long long ll;
21 const int maxn = 10 + 10;
22 const int mod = 100003;//const引用更快，宏定义也更快
23 const int INF = 1e9;
24 typedef double Matrix[maxn][maxn];
25 double a[maxn][maxn];
26 double tmp[maxn][maxn];
27 void gauss(Matrix A, int n)
28 {
29     for(int i = 0; i < n; i++)
30     {
31         int r = i;
32         for(int j = i + 1; j < n; j++)
33             if(fabs(A[j][i]) > fabs(A[r][i]))r = j;
34         if(r != i)for(int j = 0; j <= n; j++)swap(A[r][j], A[i][j]);
35         for(int k = i + 1; k < n; k++)
36         {
37             double f = A[k][i] / A[i][i];
38             for(int j = i; j <= n; j++)A[k][j] -= f * A[i][j];
39         }
40     }
41  
42     for(int i = n - 1; i >= 0; i--)
43     {
44         for(int j = i + 1; j < n; j++)
45         {
46             A[i][n] -= A[j][n] * A[i][j];
47         }
48         A[i][n] /= A[i][i];
49     }
50 }
51 int main()
52 {
53     int n;
54     scanf("%d", &n);
55     for(int i = 0; i <= n; i++)
56         for(int j = 0; j < n; j++)scanf("%lf", &a[i][j]);
57     for(int i = 1; i <= n; i++)
58     {
59         for(int j = 0; j < n; j++)
60         {
61             tmp[i - 1][n] -= a[0][j] * a[0][j];
62             tmp[i - 1][n] += a[i][j] * a[i][j];
63         }
64         for(int j = 0; j < n; j++)
65             tmp[i - 1][j] = 2.0 * (a[i][j] - a[0][j]);
66     }
67     gauss(tmp, n);
68     for(int i = 0; i < n; i++)
69     {
70         if(i == 0)printf("%.3lf", tmp[i][n]);
71         else printf(" %.3lf", tmp[i][n]);
72     }
73     puts("");
74     return 0;
75 }