球形空间产生器
# 题意
n维度空间的球形，知道球面上的n+1个点的坐标，求球心
# 题解
设半径为r，球心为(x1 , x2 , ...... , xn)
那么会满足
(a11 - x1 )2 + (a12 -x2)2 + ...... + (a1n - xn)2 = r2
(a21 - x1 )2 + (a22 -x2)2 + ...... + (a2n - xn)2 = r2
......
(an1 - x1 )2 + (an2 -x2)2 + ...... + (ann - xn)2 = r2
(an+1 1 - x1 )2 + (an+1 2 -x2)2 + ...... + (an+1 n - xn)2 = r2
有n+1个点的坐标，某两个能够消去一个xi的二次项,同时右边的r2消成0，进行n次即可
单独看某两个的x1 项相减 ,(ai1 - x1)2 - (ai+1 1 - x1)2 = 0
展开化简得到 (ai12 - ai+1 12) + 2 · x1 · (ai+1 1 -ai1) =0
即(ai12 - ai+1 12) = 2 · (ai1 - ai+1 1 ) · x1
可以由此得到每两个方程合并后的方程为 ∑(j=1:n) 2 · (aij - ai+1,j) · xj = (aij2 - ai+1 j2)
然后高斯消元得圆心坐标即可
1 #include <bits/stdc++.h> 2 #define faststd ios::sync_with_stdio(0);cin.tie(0);cout.tie(0); 3 #define ll long long 4 using namespace std; 5 const int MAXN=110; 6 const double eps=1e-6; 7 int n; 8 double a[MAXN][MAXN],b[MAXN][MAXN];// 最后一列开始存的是方程右边常数，运算后解 9 // gauss 10 // O(n^3) 11 int gauss() 12 { 13 int c, r; 14 for (c = 0, r = 0; c < n; c ++ ) 15 { 16 int t = r; 17 for (int i = r; i < n; i ++ ) 18 if (fabs(a[i][c]) > fabs(a[t][c])) 19 t = i; 20 21 if (fabs(a[t][c]) < eps) continue; 22 23 for (int i = c; i < n + 1; i ++ ) swap(a[t][i], a[r][i]); 24 for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c]; 25 26 for (int i = r + 1; i < n; i ++ ) 27 if (fabs(a[i][c]) > eps) 28 for (int j = n; j >= c; j -- ) 29 a[i][j] -= a[r][j] * a[i][c]; 30 31 r ++ ; 32 } 33 34 if (r < n) 35 { 36 for (int i = r; i < n; i ++ ) 37 if (fabs(a[i][n]) > eps) 38 return 2; 39 return 1; 40 } 41 42 for (int i = n - 1; i >= 0; i -- ) 43 for (int j = i + 1; j < n; j ++ ) 44 a[i][n] -= a[j][n] * a[i][j]; 45 46 return 0; 47 } 48 int main(){ 49 faststd 50 cin>>n; 51 for(int i=0;i<=n;i++) { 52 for (int j = 0; j <n; j++) 53 cin >> b[i][j]; 54 } 55 for(int i=0;i<n;i++) { 56 double tmp=0; 57 for (int j = 0; j < n; j++) { 58 a[i][j]=2*(b[i][j]-b[i+1][j]); 59 tmp+=b[i][j]*b[i][j]-b[i+1][j]*b[i+1][j]; 60 } 61 a[i][n]=tmp; 62 } 63 int t =gauss(); 64 for(int i=0;i<n;i++) 65 printf("%.3lf ",a[i][n]); 66 return 0; 67 }
