zoukankan      html  css  js  c++  java
  • 【BZOJ1013】[JSOI2008] 球形空间产生器(高斯消元)

    点此看题面

    大致题意: 给定一个(n)维球体上的(n+1)个点,请你求出这个球体的圆心的位置。

    列出方程

    这一看就是一道解方程题

    我们可以设这个球体的圆心的位置为((x_1,x_2,..x_n)),并设每个点到圆心的距离为(dis)

    借助题目中给出的公式,我们可以得到以下方程:

    (egin{cases}sqrt{(x_1-a_{1,1})^2+(x_2-a_{1,2})^2+...+(x_n-a_{1,n})^2}=dis\sqrt{(x_1-a_{2,1})^2+(x_2-a_{2,2})^2+...+(x_n-a_{2,n})^2}=dis\......\sqrt{(x_1-a_{n+1,1})^2+(x_2-a_{n+1,2})^2+...+(x_n-a_{n+1,n})^2}=disend{cases})

    方程的转化

    原方程看起来十分麻烦,又有平方,又有开方,很难解,因此我们要将它转化一下。

    将方程两边同时平方,可得:

    (egin{cases}(x_1-a_{1,1})^2+(x_2-a_{1,2})^2+...+(x_n-a_{1,n})^2=dis^2\(x_1-a_{2,1})^2+(x_2-a_{2,2})^2+...+(x_n-a_{2,n})^2=dis^2\......\(x_1-a_{n+1,1})^2+(x_2-a_{n+1,2})^2+...+(x_n-a_{n+1,n})^2=dis^2end{cases})

    但是,这些方程全部都是二次方程,好像非常难做。

    因此,我们考虑将每个方程展开:

    (egin{cases}x_1^2-2a_{1,1}x_1+a_{1,1}^2+x_2^2-2a_{1,2}x_2+a_{1,2}^2+...+x_n^2-2a_{1,n}x_n+a_{1,n}^2=dis^2\x_1^2-2a_{2,1}x_1+a_{2,1}^2+x_2^2-2a_{2,2}x_2+a_{2,2}^2+...+x_n^2-2a_{2,n}x_n+a_{2,n}^2=dis^2\......\x_1^2-2a_{n+1,1}x_1+a_{n+1,1}^2+x_2^2-2a_{n+1,2}x_2+a_{n+1,2}^2+...+x_n^2-2a_{n+1,n}x_n+a_{n+1,n}^2=dis^2\end{cases})

    这时候,我们显然可以看出每个方程中左边都有(x_1^2+x_2^2+...+x_n^2),右边都有(dis^2),不难想到,将第(1sim n)个方程分别减去第(n+1)个方程,便可以得到一个新的方程组,而且是一次的:

    (egin{cases}2(a_{n+1,1}-a_{1,1})·x_1+...+2(a_{n+1,n}-a_{1,n})·x_n=a_{n+1,1}^2-a_{1,1}^2+...+a_{n+1,n}^2-a_{1,n}^2\2(a_{n+1,1}-a_{2,1})·x_1+...+2(a_{n+1,n}-a_{2,n})·x_n=a_{n+1,1}^2-a_{2,1}^2+...+a_{n+1,n}^2-a_{2,n}^2\...\2(a_{n+1,1}-a_{n,1})·x_1+...+2(a_{n+1,n}-a_{n,n})·x_n=a_{n+1,1}^2-a_{n,1}^2+...+a_{n+1,n}^2-a_{n,n}^2end{cases})

    由于(a)数组是题目中给出的,我们就可以直接高斯消元了。

    代码

    #include<bits/stdc++.h>
    #define max(x,y) ((x)>(y)?(x):(y))
    #define min(x,y) ((x)<(y)?(x):(y))
    #define abs(x) ((x)<0?-(x):(x))
    #define LL long long
    #define ull unsigned long long
    #define N 10
    using namespace std;
    int n;
    namespace Gauss
    {
    	const double eps=1e-10;//eps是一个极小值,防止精度误差
    	double a[N+5][N+5],s[N+5];
    	inline void swap(double &x,double &y)
    	{
    		double t=x;x=y,y=t;
    	}
    	inline void GetData()//将读入的数据转化为方程的系数
    	{
    		register int i,j;
    		for(i=1;i<=n+1;++i) for(j=1;j<=n;++j) scanf("%lf",&a[i][j]);
    		for(i=1;i<=n;++i) for(s[i]=0,j=1;j<=n;++j) 
    			s[i]+=a[n+1][j]*a[n+1][j]-a[i][j]*a[i][j],a[i][j]=2*(a[n+1][j]-a[i][j]);
    	}
    	inline void Find_line(int x)//找到一个行数大于等于x且第x个元素系数不为0的方程,将其移至第x行
    	{
    		register int i=x,j;
    		while(i<=n&&fabs(a[i][x])<eps) ++i;
    		for(j=1;j<=n;++j) swap(a[x][j],a[i][j]);
    	}
    	inline void PrintAns()
    	{
    		register int i,j,k;
    		for(i=1;i<=n;++i)
    		{
    			for(Find_line(i),j=i+1;j<=n;++j)//消去[i+1~n]中每一行第i个元素
    			{
    				register double delta=-a[j][i]/a[i][i];
    				for(s[j]+=s[i]*delta,k=i;k<=n;++k) a[j][k]+=a[i][k]*delta;
    			}
    		}
    		for(i=n;i;--i) for(s[i]/=a[i][i],j=i-1;j;--j) s[j]-=a[j][i]*s[i];//计算出第i个未知数的值,并将第i个元素的值代入第1~i-1行的式子中消去第i个未知数
    		for(i=1;i<=n;++i) printf("%.3lf ",s[i]);//输出每一个未知数的值
    	}
    }
    int main()
    {
    	return scanf("%d",&n),Gauss::GetData(),Gauss::PrintAns(),0;
    }
    
  • 相关阅读:
    实例
    LR接口测试---webservices
    LR常用函数整理
    Codeforces Round #639 (Div. 2) A. Puzzle Pieces
    Codeforces Round #640 (Div. 4)全部七题
    POJ3177 Redundant Paths(e-DCC+缩点)
    洛谷P3469 [POI2008]BLO-Blockade(割点)
    洛谷P3275 [SCOI2011]糖果(缩点+拓扑序DP)
    POJ1236 Network of Schools(强连通分量)
    P3387 【模板】缩点(Tarjan求强连通分量)
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ1013.html
Copyright © 2011-2022 走看看