• 【THUSC2017】【LOJ2982】宇宙广播 计算几何 高斯消元


    题目大意

      有 (n)(n) 维空间中的球,求这些球的所有公切面。

      保证不会无解或有无穷多组解。

      (nleq 10)

    题解

      你可以认为这是一道传统题。

      记公切面为 (a_1x_1+a_2x_2+cdots+ a_nx_n=d),满足 (sum_ia_i^2=1)

      一个点 (x_1,x_2,ldots,x_n) 到这个面的距离为

    [frac{lvert a_1x_1+a_2x_2cdots +a_nx_n-d vert}{sqrt{a_1^2+a_2^2+cdots+a_n^2}} ]

      对于第 (i) 个球,有 (lvertsum_{j}a_jx_{i,j}-d vert=r)

      枚举所有 (lvertsum_{j}a_jx_{i,j}-d vert) 的符号,把绝对值符号去掉。

      然后解方程解出 (a_i=D_1d+D_2),带进 (sum_ia_i^2=1) 中解方程就可以解出 (d)

      时间复杂度:(O(n^32^n))

      直接解方程可能会有奇怪情况,所以要把所有坐标偏移一个值。

    代码

    const ldb eps=1e-6;
    ldb DELTA[]={0,0.235,0.3746,0.1246,0.2153,0.364,0.1253,0.6124,0.164,0.7543,0.35476};
    int n;
    ldb a[20][20];
    ldb r[20];
    ldb ans[2100][20][20];
    ldb b[20][20];
    ldb c[20][20];
    ldb e[20];
    int cnt;
    void gao(ldb d)
    {
    	for(int i=1;i<=n;i++)
    		e[i]=c[i][n+1]*d+c[i][n+2];
    	cnt++;
    	for(int i=1;i<=n;i++)
    	{
    		ldb s=0;
    		for(int j=1;j<=n;j++)
    			s+=a[i][j]*e[j];
    		s-=d;
    		for(int j=1;j<=n;j++)
    			ans[cnt][i][j]=a[i][j]-s*e[j];
    	}
    }
    void calc()
    {
    	memcpy(c,b,sizeof b);
    	int flag=0;
    	for(int i=1;i<=n;i++)
    	{
    		int x=i;
    		for(int j=i;j<=n;j++)
    			if(fabs(c[j][i])>fabs(c[x][i]))
    				x=j;
    		if(fabs(c[x][i])<eps)
    		{
    			flag=1;
    			continue;
    		}
    		if(x!=i)
    			swap(c[x],c[i]);
    		ldb v=1/c[i][i];
    		for(int j=1;j<=n+2;j++)
    			c[i][j]*=v;
    		for(int j=1;j<=n;j++)
    			if(j!=i)
    			{
    				ldb v=c[j][i];
    				for(int k=1;k<=n+2;k++)
    					c[j][k]-=c[i][k]*v;
    			}
    	}
    	if(flag)
    	{
    //		ldb d;
    //		for(int i=1;i<=n;i++)
    //			if(!c[i][i])
    //			{
    //				if(flag)
    //				{
    //				}
    //				else
    //				{
    //				}
    		printf("error
    ");
    		return;
    	}
    	ldb A=0,B=0,C=0;
    	for(int i=1;i<=n;i++)
    	{
    		A+=c[i][n+1]*c[i][n+1];
    		B+=2*c[i][n+1]*c[i][n+2];
    		C+=c[i][n+2]*c[i][n+2];
    	}
    	C-=1;
    	if(fabs(A)<eps)
    	{
    		if(fabs(B)<eps)
    		{
    			return;
    		}
    		ldb d=-C/B;
    		gao(d);
    	}
    	else
    	{
    		ldb delta=B*B-4*A*C;
    		if(delta<-eps)
    			return;
    		else
    			delta=max(delta,(ldb)0);
    		if(delta<eps)
    		{
    			ldb d=-B/(2*A);
    			gao(d);
    		}
    		else
    		{
    			ldb d1=(-B+sqrt(delta))/(2*A);
    			gao(d1);
    			ldb d2=(-B+sqrt(delta))/(2*A);
    			gao(d2);
    		}
    	}
    }
    void dfs(int x)
    {
    	if(x>n)
    	{
    		calc();
    		return;
    	}
    	for(int i=1;i<=n;i++)
    		b[x][i]=a[x][i];
    	b[x][n+1]=1;
    	b[x][n+2]=r[x];
    	dfs(x+1);
    	if(fabs(r[x])<eps)
    		return;
    	b[x][n+2]=-r[x];
    	dfs(x+1);
    }
    void solve()
    {
    	scanf("%d",&n);
    	for(int i=1;i<=n;i++)
    	{
    		for(int j=1;j<=n;j++)
    		{
    			scanf("%Lf",&a[i][j]);
    			a[i][j]-=DELTA[j];
    		}
    		scanf("%Lf",&r[i]);
    	}
    	cnt=0;
    	dfs(1);
    	printf("%d
    ",cnt);
    	for(int i=1;i<=cnt;i++)
    		for(int j=1;j<=n;j++)
    		{
    			for(int k=1;k<=n;k++)
    				printf("%.20Lf ",ans[i][j][k]+DELTA[k]);
    			printf("
    ");
    		}
    }
    int main()
    {
    //	freopen("0.in","r",stdin);
    	int t;
    	scanf("%d",&t);
    	while(t--)
    		solve();
    	return 0;
    }
    
  • 相关阅读:
    spring boot☞Swagger2文档构建及单元测试
    Spring Boot☞ 配置文件详解:自定义属性、随机数、多环境配置等
    Spring Boot☞HelloWorld开篇
    vmware workstation 12 密钥
    Spring Boot之JdbcTemplate多数据源配置与使用
    xStream完美转换XML、JSON(转)
    JAXB和XStream比较
    @RestController注解下返回到jsp视图页面
    C++中继承 声明基类析构函数为虚函数作用,单继承和多继承关系的内存分布
    mov offset和lea的区别
  • 原文地址:https://www.cnblogs.com/ywwyww/p/10294712.html
Copyright © 2020-2023  润新知