207. 球形空间产生器
题目链接
207. 球形空间产生器
有一个球形空间产生器能够在 \(n\) 维空间中产生一个坚硬的球体。
现在，你被困在了这个 \(n\) 维球体中，你只知道球面上 \(n+1\) 个点的坐标，你需要以最快的速度确定这个 \(n\) 维球体的球心坐标，以便于摧毁这个球形空间产生器。
注意： 数据保证有唯一解。
输入格式
第一行是一个整数 \(n\)。
接下来的 \(n+1\) 行，每行有 \(n\) 个实数，表示球面上一点的 \(n\) 维坐标。
每一个实数精确到小数点后 \(6\) 位，且其绝对值都不超过 \(20000\)。
输出格式
有且只有一行，依次给出球心的 \(n\) 维坐标（\(n\) 个实数），两个实数之间用一个空格隔开。
每个实数精确到小数点后 \(3\) 位。
数据范围
\(1≤n≤10\)
输入样例：
2
0.0 0.0
-1.0 1.0
1.0 0.0
输出样例：
0.500 1.500
解题思路
高斯消元
设球心坐标为 \((x_1,x_2,\dots,x_n)\)，依题意，有如下方程组（\(n+1\) 个等式，\(n+1\) 个未知数）：
\[\left\{\begin{array}{l}
\left(a_{01}-x_{1}\right)^{2}+\left(a_{02}-x_{2}\right)^{2}+\cdots+\left(a_{0 n}-x_{n}\right)^{2}=R^{2} \\
\left(a_{11}-x_{1}\right)^{2}+\left(a_{12}-x_{2}\right)^{2}+\cdots+\left(a_{1 n}-x_{n}\right)^{2}=R^{2} \\
\cdots \\
\left(a_{n1}-x_{1}\right)^{2}+\left(a_{n2}-x_{2}\right)^{2}+\cdots+\left(a_{n n}-x_{n}\right)^{2}=R^{2}
\end{array}\right.
\]
相邻两个等式作差，得：
\[\left\{\begin{array}{l}
2\left(a_{0\,1}-a_{1\,1}\right) x_{1}+2\left(a_{0\,2}-a_{1\,2}\right) x_{2}+\cdots+2\left(a_{0\,n}-a_{1\, n}\right) x_{n}=(a_{0\,1}^{2}-a_{1\,1}^{2})+(a_{0\,2}^{2}-a_{1\,2}^{2}) +\dots+(a_{0\,n}^{2}-a_{1\,n}^{2})\\
2\left(a_{1\,1}-a_{2\,1}\right) x_{1}+2\left(a_{1\,2}-a_{2\,2}\right) x_{2}+\cdots+2\left(a_{1\, n}-a_{2\, n}\right) x_{n}=(a_{1\,1}^{2}-a_{2\,1}^{2})+(a_{1\,2}^{2}-a_{2\,2}^{2}) +\dots+(a_{1\,n}^{2}-a_{2\,n}^{2})\\
\cdots \\
2\left(a_{{n-1}\,{1}}-a_{n\,1}\right) x_{1}+2\left(a_{n-1\,2}-a_{n\,2}\right) x_{2}+\cdots+2\left(a_{n-1\, n}-a_{n\, n}\right) x_{n}=(a_{n-1\,1}^{2}-a_{n\,1}^{2})+(a_{n-1\,2}^{2}-a_{n\,2}^{2}) +\dots+(a_{n-1\,n}^{2}-a_{n\,n}^{2})\\
\end{array}\right.
\]
此即为 \(n\) 元一次方程组，可先将其转化为 \(n\) 行 \(n+1\) 列的矩阵形式，再进行高斯消元：
- 先构造成如下形式（上三角矩阵）：
\( \begin{bmatrix} 1 & & & \\ 0 & 1 & & \\ 0&0 & 1 & \end{bmatrix} \)
按对角线 \((r,c)\) 枚举，为了使除数不为 \(0\)，找绝对值最大的 \(b[i][c]\)，将对应的行换到 \(r\) 行，将 \(b[r][c]\) 化成 \(1\)，每次按列将 \((r,c)\) 下面列消为 \(0\)，
- 最后构造如下形式（简化阶梯型矩阵）：
\( \begin{bmatrix} 1 &0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0& 1& \end{bmatrix} \)
按列消，从最后一列往前，最终 \(n+1\) 列即为答案
- 时间复杂度：\(O(n^3)\)
爬山法
即找到一个点，使该点到其他各点之间的距离相等，不妨先考虑一维两个点的情况：显然中点即为答案，此时中点到两个点的距离最小，引申到 \(n\) 维空间，即寻找一个点使其到已知点的距离之和最小，对一维来说，\(dx=(x-x0)^2\) 是凸函数，引申到 \(n\) 维，由于凸函数的和仍为凸函数，所以，最终形成的函数是一个单峰，可用爬山法来求解：不用求解函数具体值，只要确定一个方向，每次往峰顶走
本题可先将中心作为备选答案，求出各个点到备选点的距离，以及各个距离的平均值，如果距离大于平均值的话，对应维度需要变化，且距离相对平均值越大的话变化就越大
- 时间复杂度：\(O(随机)\)
模拟退火
一般情况，模拟退火可以取代爬山法
爬山法关键在于确定方向，即通过本身及相关性质确定方向，而模拟退火关键在于跳点，即在答案范围内随机生成一个点，如果生成的点在相关性质下要比当前点的性质要强烈，且越强烈跳到该点的概率越大
本题球心要求备选点到其他各点的距离相等，可以将这些距离求出来，即要求这些距离最小，相当于方差最小，方差为 \(0\) 时备选点即为答案
- 时间复杂度：\(O(随机)\)
代码
- 高斯消元
// Problem: 球形空间产生器
// Contest: AcWing
// URL: https://www.acwing.com/problem/content/description/209/
// Memory Limit: 64 MB
// Time Limit: 1000 ms
//
// Powered by CP Editor (https://cpeditor.org)
// %%%Skyqwq
#include <bits/stdc++.h>
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
template <typename T> void inline read(T &x) {
int f = 1; x = 0; char s = getchar();
while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
x *= f;
}
const int N=15;
double a[N][N],b[N][N];
int n;
void Gauss()
{
for(int r=1,c=1;r<=n;r++,c++)
{
int t=r;
for(int i=r+1;i<=n;i++)
if(fabs(b[i][c])>fabs(b[t][c]))t=i;
for(int i=c;i<=n+1;i++)swap(b[r][i],b[t][i]);
for(int i=n+1;i>=c;i--)b[r][i]/=b[r][c];
for(int i=r+1;i<=n;i++)
for(int j=n+1;j>=c;j--)b[i][j]-=b[i][c]*b[r][j];
}
for(int j=n;j>1;j--)
for(int i=j-1;i>=1;i--)
{
b[i][n+1]-=b[j][n+1]*b[i][j];
b[i][j]=0;
}
}
int main()
{
scanf("%d",&n);
for(int i=0;i<=n;i++)
for(int j=1;j<=n;j++)scanf("%lf",&a[i][j]);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
{
b[i][j]=2*(a[i-1][j]-a[i][j]);
b[i][n+1]+=a[i-1][j]*a[i-1][j]-a[i][j]*a[i][j];
}
Gauss();
for(int i=1;i<=n;i++)printf("%.3lf ",b[i][n+1]);
return 0;
}
- 爬山法
// Problem: 球形空间产生器
// Contest: AcWing
// URL: https://www.acwing.com/problem/content/description/209/
// Memory Limit: 64 MB
// Time Limit: 1000 ms
//
// Powered by CP Editor (https://cpeditor.org)
// %%%Skyqwq
#include <bits/stdc++.h>
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
template <typename T> void inline read(T &x) {
int f = 1; x = 0; char s = getchar();
while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
x *= f;
}
const int N=15;
double d[N][N],dist[N],res[N],delta[N];
int n;
void cal()
{
double avg=0;
for(int i=0;i<=n;i++)
{
dist[i]=delta[i]=0;
for(int j=1;j<=n;j++)dist[i]+=(d[i][j]-res[j])*(d[i][j]-res[j]);
dist[i]=sqrt(dist[i]);
avg+=dist[i]/(n+1);
}
for(int i=0;i<=n;i++)
for(int j=1;j<=n;j++)delta[j]+=(d[i][j]-res[j])*(dist[i]-avg)/avg;
}
int main()
{
scanf("%d",&n);
for(int i=0;i<=n;i++)
for(int j=1;j<=n;j++)scanf("%lf",&d[i][j]),res[j]+=d[i][j]/(n+1);
for(double T=1e4;T>1e-6;T*=0.99995)
{
cal();
for(int i=1;i<=n;i++)res[i]+=delta[i]*T;
}
for(int i=1;i<=n;i++)printf("%.3lf ",res[i]);
return 0;
}
- 模拟退火
// Problem: 球形空间产生器
// Contest: AcWing
// URL: https://www.acwing.com/problem/content/description/209/
// Memory Limit: 64 MB
// Time Limit: 1000 ms
//
// Powered by CP Editor (https://cpeditor.org)
// %%%Skyqwq
#include <bits/stdc++.h>
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
template <typename T> void inline read(T &x) {
int f = 1; x = 0; char s = getchar();
while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
x *= f;
}
const int N=15;
double d[N][N],dist[N],res[N],cur[N],t[N],min_dist=1e9;
int n;
double rand(double l,double r)
{
return (double)rand()/RAND_MAX*(r-l)+l;
}
double cal(double a[])
{
double avg=0;
for(int i=0;i<=n;i++)
{
dist[i]=0;
for(int j=1;j<=n;j++)dist[i]+=(d[i][j]-a[j])*(d[i][j]-a[j]);
avg+=dist[i]/(n+1);
}
double ret=0;
for(int i=0;i<=n;i++)ret+=(dist[i]-avg)*(dist[i]-avg);
if(min_dist>ret)
{
min_dist=ret;
memcpy(res,cur,sizeof cur);
}
return ret;
}
void simulate_anneal()
{
for(double T=1e4;T>1e-6;T*=0.9999)
{
for(int i=1;i<=n;i++)t[i]=cur[i]+rand(-T,+T);
double delta=cal(t)-cal(cur);
if(exp(-delta/T)>rand(0,1))memcpy(cur,t,sizeof t);
}
}
int main()
{
scanf("%d",&n);
srand(time(NULL));
for(int i=0;i<=n;i++)
for(int j=1;j<=n;j++)
{
scanf("%lf",&d[i][j]);
cur[j]+=d[i][j]/(n+1);
}
double s=clock();
while(((double)clock()-s)/CLOCKS_PER_SEC<0.7)simulate_anneal();
for(int i=1;i<=n;i++)printf("%.3lf ",res[i]);
return 0;
}
