题解：P17213 [ICPC 2017 Nanning R] The Ball

· · 题解

题意简述

给定若干三维半空间及非负卦限， 求它们交集内能完整放置的最大球半径。 若半径没有上界，输出 Infinity。

解题思路

设球心为 c=(x,y,z)，半径为 r。 考虑一个半空间：

AX+BY+CZ\leq D

记法向量为 v=(A,B,C)。 球中任意点可以写成 c+u，其中 \lVert u\rVert\leq r。 由柯西不等式，线性函数在球上的最大值为：

v\cdot c+\max_{\lVert u\rVert\leq r}v\cdot u =v\cdot c+r\lVert v\rVert

最大值在 u 与 v 同向时取得。 所以整颗球位于该半空间内，当且仅当：

Ax+By+Cz+r\sqrt{A^2+B^2+C^2}\leq D

球还要位于非负卦限内。 球到三个坐标面的距离分别为 x,y,z，因此需要：

\begin{aligned} r & \leq x \\ r & \leq y \\ r & \leq z \end{aligned}

于是原问题等价于以下线性规划：

\begin{aligned} \max & r \\ \text{s.t.} & A_ix+B_iy+C_iz+\sqrt{A_i^2+B_i^2+C_i^2}r\leq D_i \\ & -x+r\leq0 \\ & -y+r\leq0 \\ & -z+r\leq0 \\ & x,y,z,r\geq0 \end{aligned}

线性规划只有 4 个原变量， 但约束右端项可能为负，不能直接从松弛变量基开始。 因此使用两阶段单纯形法。

第一阶段加入一个人工变量。 先选右端项最小的约束，让人工变量入基， 再以消除人工变量为目标进行换基。 若第一阶段最优值仍小于 0，说明约束组不可行。 题目保证交集具有正体积，实际输入不会得到这一结果。

得到可行基后，把仍在基中的人工变量换出， 再进入以 r 为目标的第二阶段。 每次选择检验数为负的非基变量入基， 并在所有该列系数为正的约束中执行最小比值检验。

入基变量和离基变量都按变量编号打破平局， 即使用 Bland 规则，避免退化情况下循环。 若某个能改善目标的入基列没有正系数， 该变量可以无限增大，线性规划无界，输出 Infinity。

上面的半空间约束来自球上线性函数的精确最大值， 三个坐标面约束也恰好等价于球不越过坐标面。 因此线性规划的可行解与题目中的可行球一一对应， 目标值正是球半径。

一次主元变换只是从一条等式中解出入基变量， 再代入其余等式，不会改变可行域。 最小比值检验保证新的基变量仍然非负。 当所有检验数均非负时， 任何非基变量从零开始增大都不能继续改善目标， 当前基可行解就是最优解。

表中只有 6 列，一次主元变换的时间复杂度为 O(n)。 单纯形法的理论最坏复杂度为指数级， 但本题只有 4 个原变量和至多 103 条约束。

参考代码

#include <bits/stdc++.h>
using namespace std;

using ld=long double;
const int N=105;
const int M=6;
const ld eps=1e-12L;
const ld inf=1e100L;
ld a[N][M];
int row[N],col[M],m;
void pivot(int x,int y)
{
    ld v=a[x][y];
    for(int i=0;i<=m+1;i++)
    {
        if(i==x)continue;
        for(int j=0;j<M;j++)
        {
            if(j==y)continue;
            a[i][j]-=a[x][j]*a[i][y]/v;
        }
    }
    for(int j=0;j<M;j++)if(j!=y)a[x][j]/=v;
    for(int i=0;i<=m+1;i++)if(i!=x)a[i][y]/=-v;
    a[x][y]=1/v;
    swap(row[x],col[y]);
}
bool simplex(int phase)
{
    int obj=phase==1?m+1:m;
    while(1)
    {
        int y=-1;
        for(int j=0;j<M-1;j++)
        {
            if(phase==2&&col[j]==-1)continue;
            if(a[obj][j]<-eps&&(y==-1||col[j]<col[y]))y=j;
        }
        if(y==-1)return 1;
        int x=-1;
        for(int i=0;i<m;i++)
        {
            if(a[i][y]<=eps)continue;
            if(x==-1)
            {
                x=i;
                continue;
            }
            ld p=a[i][M-1]/a[i][y];
            ld q=a[x][M-1]/a[x][y];
            if(p<q-eps||(abs(p-q)<=eps&&row[i]<row[x]))x=i;
        }
        if(x==-1)return 0;
        pivot(x,y);
    }
}
ld linear()
{
    int x=0;
    for(int i=1;i<m;i++)if(a[i][M-1]<a[x][M-1])x=i;
    if(a[x][M-1]<-eps)
    {
        pivot(x,4);
        if(!simplex(1)||a[m+1][M-1]<-eps)return -inf;
        for(int i=0;i<m;i++)
        {
            if(row[i]!=-1)continue;
            int y=0;
            for(int j=1;j<M-1;j++)
            {
                if(abs(a[i][j])>abs(a[i][y]))y=j;
            }
            if(abs(a[i][y])>eps)pivot(i,y);
        }
    }
    if(!simplex(2))return inf;
    return a[m][M-1];
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int T;
    cin>>T;
    while(T--)
    {
        memset(a,0,sizeof a);
        int n;
        cin>>n;
        m=n+3;
        for(int i=0;i<n;i++)
        {
            int x,y,z,d;
            cin>>x>>y>>z>>d;
            a[i][0]=x;
            a[i][1]=y;
            a[i][2]=z;
            a[i][3]=sqrtl((ld)x*x+(ld)y*y+(ld)z*z);
            a[i][M-1]=d;
        }
        for(int i=0;i<3;i++)
        {
            a[n+i][i]=-1;
            a[n+i][3]=1;
        }
        for(int i=0;i<m;i++)
        {
            row[i]=4+i;
            a[i][4]=-1;
        }
        for(int i=0;i<4;i++)col[i]=i;
        col[4]=-1;
        a[m][3]=-1;
        a[m+1][4]=1;
        ld ans=linear();
        if(ans>=inf/2)cout<<"Infinity";
        else cout<<fixed<<setprecision(4)<<max(ans,0.0L);
        cout<<'\n';
    }
    return 0;
}