【BZOJ】1013 球形空间生成器

最新推荐文章于 2020-04-15 19:00:42 发布

y20070316

最新推荐文章于 2020-04-15 19:00:42 发布

阅读量331

点赞数

分类专栏：计算几何高斯消元文章标签：高斯消元 OI BZOJ 计算几何

本文链接：https://blog.csdn.net/y20070316/article/details/50614262

版权

计算几何同时被 2 个专栏收录

3 篇文章 0 订阅

订阅专栏

高斯消元

1 篇文章 0 订阅

订阅专栏

该博客介绍了如何使用高斯消元法解决BZOJ问题1013，即球形空间生成器。通过建立数学模型，分析了如何根据圆心坐标和圆上两点的坐标，形成等式并进行求解。博主指出，尽管时间复杂度为O(n^3)，但当n=100时，通过避免浮点误差，该方法仍然可行。

摘要由CSDN通过智能技术生成

分析

设圆心的坐标为 $(r_1,r_2,...,r_n)$ ，
点 $A(a_1,a_2,...,a_n)$ 和点 $B(b_1,b_2,...,b_n)$ 在圆上。

则有
$dist=\sqrt{(r_1-a_1)^2+...+(r_n-a_n)^2}=\sqrt{(r_1-b_1)^2+...+(r_n-b_n)^2}$
$\therefore (r_1-a_1)^2+...+(r_n-a_n)^2=(r_1-b_1)^2+...+(r_n-b_n)^2$
$\therefore \sum {r_i}^2-\sum 2a_ir_i+\sum {a_i}^2=\sum {r_i}^2-\sum 2b_ir_i+\sum {b_i}^2$
$\therefore \sum (2b_i-2a_i)r_i=\sum ({b_i}^2)-\sum ({a_i}^2)$

一共有 $n+1$ 个点，可以列出 $n$ 个等式。
可以预处理出 $\sum ({a_i}^2)$ ，每个等式用 $O(n)$ 可求出。
用高斯消元即可解决。

时间复杂度： $O(n^3)$
空间复杂度： $O(n^2)$

所以数据出到 $n=100$ 其实也是可以的，只要能避免浮点误差带来的太大影响。

代码

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
using namespace std;

const int N=15;
const double eps=1e-5;

int n; double x[N][N];
double a[N][N],sum[N],res[N];

void init(void)
{
    scanf("%d",&n);
    for (int i=1;i<=n+1;i++)
        for (int j=1;j<=n;j++)
            scanf("%lf",&x[i][j]);

    for (int i=1;i<=n+1;i++)
        for (int j=1;j<=n;j++)
            sum[i]+=x[i][j]*x[i][j];
    for (int i=1;i<=n;i++)
    {
        for (int j=1;j<=n;j++)
            a[i][j]=x[i+1][j]-x[i][j];
        a[i][n+1]=(sum[i+1]-sum[i])/2;
    }
}

inline int cmp(double i,double j)
{
    if (fabs(i-j)<eps) return 0;
    return i<j?-1:1;
}

inline void swap(int i,int j)
{
    for (int k=1;k<=n+1;k++) a[i][k]+=a[j][k];
    for (int k=1;k<=n+1;k++) a[j][k]=a[i][k]-a[j][k];
    for (int k=1;k<=n+1;k++) a[i][k]-=a[j][k];
}

void gauss(void)
{
    double r;
    for (int i=1;i<=n;i++)
    {
        for (int j=i+1;j<=n;j++)
            if (!cmp(a[i][i],0)||cmp(abs(a[i][i]),abs(a[j][i]))>0) swap(i,j);
        for (int j=i+1;j<=n;j++)
            if (cmp(a[i][j],0))
            {
                r=a[j][i]/a[i][i];
                for (int k=i;k<=n+1;k++) a[j][k]-=a[i][k]*r;
            }
    }
    for (int i=n;i;i--)
    {
        for (int j=i+1;j<=n;j++) a[i][n+1]-=a[i][j]*res[j];
        res[i]=a[i][n+1]/a[i][i];
    }

    for (int i=1;i<n;i++) printf("%0.3lf ",res[i]);
    printf("%0.3lf\n",res[n]);
}

int main(void)
{   
    init();
    gauss();    

    return 0;
}