【CSP】202009-4 星际旅行 90%

题目大意：

n维空间内有一半径为r的球体，空间中球体之外有m个点，在不穿过球体的条件下求这m个点两两间的最短曲线距离。

分析：

显然有两种情况：1.两点连线不经过球体；2.两点连线穿过球体。

第一种情况显然，考虑第二种情况：将球心、两点作为研究平面，可以发现最短曲线一定包括两条线段和一段圆弧。由两点间直线最短可知，如果圆弧左右端点如果不是切点的话一定会形成三条（曲）边，和一定大于切线段长度。

之后，主要使用余弦定理和反三角函数就可以处理问题，这部分主要考察代码而不是思维能力（所以当我知道第一次编译并且编译过的时候只剩1个bug的时候我感觉我很强大（逃））

需要注意的是，cmath库中的acos对于浮点数运算自带的误差不能很好的处理，当数值接近1或-1的时候要手动处理一下，防止微小越界出现NaN的情况。

最后，由于是计算两两间距离，可以用数组存储结果进行一点小优化（虽然也只是85->90）。

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <vector>
#include <set>
#include <string.h>
#include <cmath>

#define up(l,r,i) for(int i=l;i<=r;i++)
#define dn(l,r,i) for(int i=r;i>=l;i--)

typedef long long ll;

using namespace std;

inline int _max(const int& a,const int& b){return a>b?a:b;}
inline int _min(const int& a,const int& b){return a<b?a:b;}

#define PI 3.14159265

const int MAXN = 120;
const int MAXM = 2020;

inline bool equal(const double& a, const double &b){
    if(a >= b)
        return ((a-b)<0.00000000001)?1:0;
    return ((b-a)<0.00000000001)?1:0;
}

int n,m;
double r;

struct Vec{
    double w[MAXN];
    double mod,dmod;
    Vec():mod(0),dmod(0){memset(w,0,sizeof(w));}
    void print(){
        printf("[");
        up(1,n,i){
            printf("%.4f, ",w[i]);
        }
        printf("]\n");
    }
}v[MAXM];

inline double mod(Vec& a){
    double ret = 0;
    up(1,n,i)
        ret += a.w[i]*a.w[i];
    ret = sqrt(ret);
    a.mod = ret;
    if(a.mod) a.dmod = 1/a.mod;
    return a.mod;
}

Vec operator - (const Vec& x){
    Vec ret;
    up(1,n,i)
        ret.w[i] = -x.w[i];
    ret.mod = x.mod;
    ret.dmod = x.dmod;
    return ret;
}

Vec operator - (const Vec& x,const Vec& y){
    Vec ret;
    up(1,n,i)
        ret.w[i] = x.w[i] - y.w[i];
    mod(ret);
    return ret;
}

double dot(const Vec& a, const Vec& b){
    double ret = 0.f;
    up(1,n,i)
        ret += a.w[i]*b.w[i];
    return ret;
}

double getangle(const Vec& a, const Vec& b){
    double Cos = dot(a,b)*a.dmod*b.dmod;
    double ret = acos(Cos);
    if(equal(Cos,1.f) || equal(Cos,-1.f)){
        ret = (PI-Cos*PI)/2;
    }
    return ret;
}

inline double r2a(const double& a){
    return a*180/PI;
}

inline double pow2(const double& a){
    return a*a;
}

bool ifcross(const Vec& a, const Vec& b){//判断两点连线是否与圆相交
    Vec l = b-a;//如果某点与圆心连线与两点连线的夹角大于90度那么不可能相交
    double a1 = getangle(-a,l),a2 = getangle(-b,-l);
    //printf("某点与圆心连线与两点连线的夹角分别为%.4f,%.4f\n",r2a(a1),r2a(a2));
    if(a1 > PI/2 || a2 > PI/2) return false;
    double d = a.mod*sin(a1);
    if(d <= r) return true;
    return false;
}

double solve(const Vec& a, const Vec& b){
    if(!ifcross(a,b))
        return (a-b).mod;
    double aa,a1,a2,rem,d1,d2,ans = 0;
    aa = getangle(a,b);
    a1 = acos(r/a.mod);
    a2 = acos(r/b.mod);
    rem = aa-a1-a2;
    rem *= r;
    d1 = sqrt(pow2(a.mod)-pow2(r));
    d2 = sqrt(pow2(b.mod)-pow2(r));
    ans = rem+d1+d2;
    //printf("与圆心连线夹角为%.4f,切线夹角分别为%.4f,%.4f\n",aa,a1,a2);
    return ans;
}

double ans[MAXM][MAXM];

int main()
{
    //freopen("y.in","r",stdin);

    ios::sync_with_stdio(false);

    cin>>n>>m;
    cin>>r;
    Vec o;
    up(1,n,i)
        cin>>o.w[i];
    
    up(1,m,i){
        up(1,n,j){
            cin>>v[i].w[j];
        }

        v[i] = v[i] - o;
        mod(v[i]);
    }

    up(1,m,i){
        double sum = 0;
        up(1,m,j){
            if(i == j) continue;
            //printf("正在计算P%d,P%d之间的最短距离\n",i,j);
            if(j < i) sum += ans[j][i];
            else {ans[i][j] = solve(v[i],v[j]);
                sum += ans[i][j];
            }
        }
        printf("%.12f\n",sum);
    }




    return 0;
}