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;
}

posted @ 2022-04-14 17:32 zyy2001 阅读(68) 评论(0) 编辑收藏举报

刷新页面返回顶部

acwing_zyy

别人躺着我站着，别人走着我跑着。

207. 球形空间产生器

题目链接