第20次CSP认证第四题星际旅行(计算几何分类讨论超详解)

本文链接：https://blog.csdn.net/tongjingqi_/article/details/121115789

官网链接：http://118.190.20.162/view.page?gpid=T110

在这里插入图片描述
题意：题面很简练，就是计算不经过一个圆，求任意两点间最短曲线距离。
思路：分为两种情况，存在与圆周重合的轨迹，不存在与圆周重合的轨迹。
我们可以这样思考，如果存在与圆周重合的轨迹，那么两点的连线一定会进入圆的内部，问题就转化为判断一个线段是否穿过圆。这时我们自然会想到，求圆心和线段所在直线的距离，如果距离大于半径，那一定不会和穿过圆。不过注意，这个是不会穿过圆的充分条件，不是必要条件！！！我们要找到充要条件，才能保证不重不漏。
不是必要条件是因为我们的线段不等于直线，而仅仅是直线的一部分，如果整条直线都没有穿过圆，那么线段也一定不会穿过圆，反之不一定成立。
所以我们还需要想到直线穿过圆，但是线段不穿过圆的情况，首先构造一条穿过圆的直线，然后我们可以观察哪些线段是不会穿过圆的，显然，直线与圆交点处向外做射线的部分，是这样的线段。
在这里插入图片描述

我们再次思考这样的线段满足怎样的条件，我们可以发现，该射线可以绕发出点旋转而不进入圆内，活动范围是180°！
于是我们可以想到，如果要满足这个条件，那么设线段两端点为A,B，圆心为O点，∠ABO和∠BAO必然有一个大于等于90°。
为了验证我们找到的是否是充要条件，这里可以联想到计算机图形学中的线段窗口裁剪算法，我们可以把中心为圆的区域划分成10个区域，然后讨论两个点落在不同区域时所连成线段与圆的关系是否满足刚才的判断条件
在这里插入图片描述
讨论清楚两种情况后，接下来就是用数学公式实现，求点到直线距离可以先求面积然后除以底，求曲线长度可以用余弦定理求出大角之后减去两个小角，详细可见代码。

#include <iostream>
#include<math.h>
#include<string.h>
typedef long long ll;
const double pi=acos(-1);
const int maxn=100+5;
const int maxm=2000+100;
using namespace std;
ll N,M,R;
ll O[maxn],p[maxm][maxn];
double g[maxm][maxm];
ll squdis(ll *a,ll *b){
    ll sum=0;
    for(ll i=0;i<N;i++){
        sum+=(a[i]-b[i])*(a[i]-b[i]);
    }
    return (sum);
}
int main()
{

    cin>>N>>M>>R;
    for(int i=0;i<N;i++){//黑洞不一定在原点！
        scanf("%lld",&O[i]);
    }
    for(int i=0;i<M;i++){
        for(int j=0;j<N;j++){
            scanf("%lld",&p[i][j]);
            p[i][j]-=O[j];
        }
    }
    memset(O,0,sizeof O);
    for(ll i=0;i<M;i++){
  //      double sum=0;
        for(ll j=i+1;j<M;j++){
            if(i==j)continue;
            ll
            doi2=squdis(O,p[i]),
            doj2=squdis(O,p[j]),
            dij2=squdis(p[i],p[j]);
            double
            doi=sqrt((double)(doi2)),
            doj=sqrt((double)doj2),
            dij=sqrt((double)dij2);
            double tp=(double)(doi+doj+dij)/2.0;
            double S=sqrt(tp*(tp-doi)*(tp-doj)*(tp-dij));
            double H=S*2/dij;
            if(H>=R||dij2+doi2<=doj2||dij2+doj2<=doi2){
                g[i][j]=g[j][i]=dij;
                continue;
            }
            double disi2=doi2-R*R,disi=sqrt((double)disi2);
            double disj2=doj2-R*R,disj=sqrt((double)disj2);
            double theta;
            if(fabs(doi+doj-dij)<1e-9)theta=pi;
            else theta=acos((double)(doi2+doj2-dij2)/(2.0*doi*doj));
            //cout<<endl<<(double)(doi2+doj2-dij2)/(2.0*doi*doj)<<"&&"<<theta<<"*"<<pi<<endl;
            double theta1=acos((double)R/doi),theta2=acos((double)R/doj);
            double alpha=theta-theta1-theta2;
            double disc=R*alpha;
            g[i][j]=g[j][i]=disc+disi+disj;
        }
    }
    for(int i=0;i<M;i++){
        double sum=0;
        for(int j=0;j<M;j++){
            if(i==j)continue;
            sum+=g[i][j];
        }
        printf("%.14f\n",double(sum));
    }
    return 0;
}