Time Limit: 1 Sec Memory Limit: 162 MB
洛谷:
Description
有一个球形空间产生器能够在n维空间中产生一个坚硬的球体。现在,你被困在了这个n维球体中,你只知道球
面上n+1个点的坐标,你需要以最快的速度确定这个n维球体的球心坐标,以便于摧毁这个球形空间产生器。
Input
第一行是一个整数n(1<=N=10)。接下来的n+1行,每行有n个实数,表示球面上一点的n维坐标。每一个实数精确到小数点
后6位,且其绝对值都不超过20000。
Output
有且只有一行,依次给出球心的n维坐标(n个实数),两个实数之间用一个空格隔开。每个实数精确到小数点
后3位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。
Sample Input
0.0 0.0
-1.0 1.0
1.0 0.0
Sample Output
HINT
提示:给出两个定义: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}+cdots +(a_{n}-b_{n})^{2}}$
n维空间确实不太好想,但我们可以先考虑二维空间,二维空间中给我们圆上的两个点,那么我们可以找到此两点连线的中垂线到这两点的距离一样,所以我们还需要第三个点来确定圆心。那么我们设圆心的坐标为$(x_{1}, x_{2}, …, x_{n})$,其半径设为$r$。那么总共n+1个方程,n+1个未知数。
然后我们来列个方程:$sum_{i=1}^{n}(x_{i}-a_{i})^{2}=r^2$,有n+1个方程,解多元方程的话我们会高斯消元,但这是线性的,这个多元元2次方程不是线性的,所以我们只能想法了。
我们将相邻两个式子相减就会得到n个不含r的方程:$sum_{j=1}^{n}(o_{j}-x_{i,j})^{2}-sum_{j=1}^{n}(o_{j}-x_{i+1,j})^{2}=0$
整理可得n个不含二次项的方程:$sum_{j=1}^{n}(o_{j}-x_{i,j})^{2}-(o_{j}-x_{i+1,j})^{2}=sum_{j=1}^{n}x_{i,j}^{2}-x_{i+1,j}^{2}+2o_{j}x_{i+1,j}-2o_{j}x_{i,j}$其中o表示圆心
移项得$sum_{j=1}^{n} 2o_{j}x_{i+1,j}-2o_{j}x_{i,j}=sum_{i=1}^{n}x_{i+1,j}^{2}-x_{i,j}^{2}$,即$o_{j}$的系数为$2(x_{i+1,j}-x_{i,j})$
用线性方程就可以解答了。那么我们可以用高斯消元来解决这题,实际上高斯消元就是把矩阵化成上三角矩阵而已。然后就是一个个减下来
当然这是一种做法,还有一种比较麻烦的做法就是用模拟退火,这个做法之所以比较麻烦就是调参的时候可能会WA掉或者T掉十来发。。。这题在网上似乎也并没有完美用模拟退火A了的解答。
模拟退火的具体讲解:https://www.cnblogs.com/no-true/p/9737193.html
https://blog.csdn.net/daaikuaichuan/article/details/81381875
随机一个点当圆心,然后和每个点比较,求出平均距离r,如果到这个点的距离大于r,说明离这个点远了,就给圆心施加一个向这个点的力;
若小于r,说明近了,就施加一个远离这个点的力。所有点比较完后,把假设的圆心按合力方向移动一个距离,距离和当前温度(T)有关。时间越久,温度越低
T要从一个比较大的数(为了开始时ans能尽快移动到圆心附近)到小于精度的数(保证答案精度),迭代次数尽量多(就是不TLE的情况下T每次乘的数尽量接近1)
为了是结果更精确,我们的距离不开方。
以下是AC代码:
#include <bits/stdc++.h> using namespace std; const double esp=1e-8; double g[20][20],mat[20][20]; double pw(double a) { return a*a; } int main() { int n; scanf ("%d",&n); for (int i=1; i<=n+1; i++) for (int j=1; j<=n; j++) scanf ("%lf",&g[i][j]); double s; for (int i=1; i<=n; i++){ s=0; for (int j=1; j<=n; j++){ mat[i][j]=2*(g[i+1][j]-g[i][j]); s+=pw(g[i+1][j])-pw(g[i][j]); } mat[i][n+1]=s; } for (int i=1; i<=n; i++){ int now=i; for (int j=i+1; j<=n; j++) if (fabs(mat[j][i])>fabs(mat[now][i])) now=j; if (now!=i) swap(mat[now],mat[i]); for (int k=i+1; k<=n; k++){ double p=mat[k][i]/mat[i][i]; for (int j=i; j<=n+1; j++) mat[k][j]-=mat[i][j]*p; } } for (int i=n; i>=1; i--){ for (int j=i+1; j<=n; j++) mat[i][n+1]-=mat[j][n+1]*mat[i][j]; mat[i][n+1]/=mat[i][i]; } for (int i=1; i<=n; i++) printf("%.3f%c",mat[i][n+1],i==n?' ':' '); return 0; }
模拟退火洛谷AC代码:(BZOJ上会T)
#include <bits/stdc++.h> using namespace std; const double inf=10000; const double esp=1e-5; const double lim=0.9999; double p[20][20],ans[20],f[20],dis[20]; double dist(double *a,double *b,int n) { double sum=0; for (int i=1; i<=n; i++) sum+=(a[i]-b[i])*(a[i]-b[i]); return sum; } int main() { int n; scanf ("%d",&n); for (int i=1; i<=n+1; i++) for (int j=1; j<=n; j++) scanf ("%lf",&p[i][j]); for (int i=1; i<=n; i++) ans[i]=0; for (double T=inf; T>=esp; T*=lim){ double avg=0; for (int i=1; i<=n+1; i++){ dis[i]=dist(ans,p[i],n); avg+=dis[i]; } avg/=n+1; for (int i=1; i<=n; i++) f[i]=0; for (int i=1; i<=n+1; i++) for (int j=1; j<=n; j++) f[j]+=(dis[i]-avg)*(p[i][j]-ans[j])/avg; for (int i=1; i<=n; i++) ans[i]+=f[i]*T; } for (int i=1; i<=n; i++) printf("%.3f ",ans[i]); }