【THUSC2017】【LOJ2982】宇宙廣播 計算幾何 高斯消元

題目大意

　　有 \(n\) 個 \(n\) 維空間中的球，求這些球的全部公切面。

　　保證不會無解或有無窮多組解。

　　\(n\leq 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_2\cdots +a_nx_n-d\rvert}{\sqrt{a_1^2+a_2^2+\cdots+a_n^2}} \]
　　對於第 \(i\) 個球，有 \(\lvert\sum_{j}a_jx_{i,j}-d\rvert=r\)。

　　枚舉全部 \(\lvert\sum_{j}a_jx_{i,j}-d\rvert\) 的符號，把絕對值符號去掉。

　　而後解方程解出 \(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\n");
        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\n",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("\n");
        }
}
int main()
{
//  freopen("0.in","r",stdin);
    int t;
    scanf("%d",&t);
    while(t--)
        solve();
    return 0;
}