[ICPC 2017 ECNA B]Craters

Link

Solution

求出所有可能在凸包上的切线,然后半平面交。。。
二维量的eps比较微妙

Code

const int XN=1e5+11;
const double eps=1e-4,inf=1e100,pi=acos(-1);

int sgn(double const &x) {
    return (x>-eps)-(x<eps);
}

struct Point {
    double x,y;

    double Length() const {
        return sqrt(x*x+y*y);
    }

    double Length2() const {
        return x*x+y*y;
    }

    Point Unit() const {
        double len=sqrt(x*x+y*y);
        return {x/len,y/len};
    }

    Point TurnRight() const {
        return {y,-x};
    }

    friend Point operator +(Point const &a,Point const &b) {
        return {a.x+b.x,a.y+b.y};
    }

    friend Point operator -(Point const &a,Point const &b) {
        return {a.x-b.x,a.y-b.y};
    }

    friend Point operator *(Point const &p,double const &k) {
        return {p.x*k,p.y*k};
    }

    friend double Inner(Point const &a,Point const &b) {
        return a.x*b.x+a.y*b.y;
    }

    friend double Outer(Point const &a,Point const &b) {
        return a.x*b.y-a.y*b.x;
    }

    friend double Dist(Point const &a,Point const &b) {
        return sqrt(Inner(a-b,a-b));
    }

    friend bool operator <(Point const &a,Point const &b) {
        return a.x<b.x || (a.x==b.x && a.y<b.y);
    }

    friend bool operator ==(Point const &a,Point const &b) {
        return a.x==b.x && a.y==b.y;
    }
}p[XN];

struct Circle {
    Point o;
    double r;

    bool Include(Circle const &a) {
        return sgn(a.r+Dist(o,a.o)-r)<=0;
    }
}o[XN];

struct Line {
    Point p,v;

    friend bool operator <(Line const &a,Line const &b) {
        return a.p<b.p || (a.p==b.p && a.v<b.v);
    }
};

double Dist(Point const &a,Line const &l) {
    return fabs(Outer(a-l.p,l.v))/l.v.Length();
}

bool OnLeft(Point const &p,Line const &ln) {
    return sgn(Outer(ln.v,p-ln.p))>0;
}

bool Paral(Line const &ln1,Line const &ln2) {
    return sgn(Outer(ln1.v,ln2.v))==0;
}

Point Cross(Line const &ln1,Line const &ln2) {
    if(Paral(ln1,ln2)) {
        return ln1.p+ln1.v*inf;
    } else {
        double t=Outer(ln2.v,ln1.p-ln2.p)/Outer(ln1.v,ln2.v);
        return ln1.p+ln1.v*t;
    }
}

Point CutPoint(Point const &p,Circle const &a) {
    if(sgn(Dist(p,a.o)-a.r)==0)
        return p;
    else {
        double x=Dist(p,a.o),cosa=a.r/x,sina=sqrt(1-cosa*cosa);
        Point v=(p-a.o).Unit(),
              t={v.x*cosa-v.y*sina,v.y*cosa+v.x*sina};
        return a.o+t*a.r;
    }
}

Point Mirror(Point const &p,Line const &l) {
    Point foot=l.p+l.v.Unit()*(Inner(p-l.p,l.v)/l.v.Length());
    return foot+(foot-p);
}

Line Cut(Circle const &a,Circle const &b) {
    /*a--b   b|
      ----   a|
      */
    if(sgn(a.r-b.r)==0) {
        return {a.o+(b.o-a.o).Unit().TurnRight()*a.r,b.o-a.o};
    } else {
        double d=Dist(a.o,b.o),
               x=a.r*d/(b.r-a.r);
        Point c=a.o+(a.o-b.o).Unit()*x,
              p=CutPoint(c,a);
        if(a.r>b.r) {
            p=Mirror(p,{a.o,b.o-a.o});
            return {p,c-p};
        } else
            return {c,p-c};
    }
}

std::vector<Line> Intersect(std::vector<Line> &ln) {
    std::sort(ln.begin(),ln.end(),[](Line const &a,Line const &b)->bool {
                auto Quad=[](Point const &a)->int {
                    return sgn(a.x)>0 && sgn(a.y)>=0?1:
                            (sgn(a.x)<=0 && sgn(a.y)>0?2:
                            (sgn(a.x)<0 && sgn(a.y)<=0?3:4));
                };
                int qa=Quad(a.v),qb=Quad(b.v);
                return qa==qb?sgn(Outer(b.v,a.v))<0:qa<qb;
            });
    static Point Qp[XN];static Line Qln[XN];
    int n=ln.size();
    int head,tail;Qln[head=tail=1]=ln[0];
    for(int i=1;i<n;++i) {
        while(tail-head>=1 && !OnLeft(Qp[tail-1],ln[i]))
            tail--;
        while(tail-head>=1 && !OnLeft(Qp[head],ln[i]))
            head++;
        Qln[++tail]=ln[i];
        if(Paral(Qln[tail-1],Qln[tail])){
            int d1=OnLeft(Qln[tail].p,Qln[tail-1]),d2=OnLeft(Qln[tail-1].p,Qln[tail]);
            if(!d1 || !d2) {
                tail--;
                if(d1)
                    Qln[tail]=ln[i];
            }
        }
        if(tail-head>=1)
            Qp[tail-1]=Cross(Qln[tail-1],Qln[tail]);
    }
    while(tail-head>=1 && !OnLeft(Qp[tail-1],Qln[head]))
        tail--;
    if(tail-head>=1)
        return std::vector<Line>(Qln+head,Qln+tail+1);
    else
        return {};
}

double ArcLen(double const &len,double const &r) {
    return r*asin(len/2/r)*2;
}

int main() {
    freopen("input","r",stdin);
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);
    srand(0x1f1f1f1f);
    int n;std::cin>>n;
    for(int i=1;i<=n;++i) {
        std::cin>>o[i].o.x>>o[i].o.y>>o[i].r;
        o[i].r+=10;
    }
    static bool ban[XN];
    for(int i=1;i<=n;++i)
        if(!ban[i]) {
            for(int j=1;j<=n;++j)
                if(i!=j && o[i].Include(o[j]))
                    ban[j]=1;
        }
    int cc=0;
    for(int i=1;i<=n;++i)
        if(!ban[i])
            o[++cc]=o[i];
    n=cc;
    double Ans=0;
    if(n==1)
        Ans=2*pi*o[1].r;
    else {
        std::map<Line,std::pair<int,int>> id;
        std::vector<Line> v;
        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j) 
                if(i!=j) {
                    Line cut=Cut(o[i],o[j]);
                    /*
                    for(auto x : {o[i],o[j]})
                        std::cerr<<Dist(x.o,cut)<<' '<<x.r<<'\n';
                    */
                    bool tak=1;
                    for(int k=1;k<=n;++k)
                        tak&=OnLeft(o[k].o,cut) && sgn(Dist(o[k].o,cut)-o[k].r)>=0;
                    if(tak)
                        v.push_back(cut);
                }
        v=Intersect(v);
        int s=v.size();
        std::vector<std::array<Point,2>> cp(s);
        std::vector<int> ids(s);
        for(int i=0;i<s;++i) {
            double mx=-inf,mn=inf;
            for(int j=1;j<=n;++j) {
                Point pts[2];
                pts[1]=Mirror(pts[0]=CutPoint(v[i].p,o[j]),{v[i].p,o[j].o-v[i].p});
                for(int d=0;d<2;++d)
                    if(sgn(Outer(pts[d]-v[i].p,v[i].v))==0) {
                        double t=Inner(pts[d]-v[i].p,v[i].v)/v[i].v.Length2();
                        if(Enlarge(mx,t)) {
                            cp[(i+1)%s][0]=pts[d];
                            ids[(i+1)%s]=j;
                        }
                        if(Reduce(mn,t)) {
                            cp[i][1]=pts[d];
                        }
                        break;
                    }
            }
        }
        for(int i=0;i<s;++i)
            Ans+=Dist(cp[i][1],cp[(i+1)%s][0])+ArcLen(Dist(cp[i][0],cp[i][1]),o[ids[i]].r);
    }
    std::cout<<std::fixed<<std::setprecision(10)<<Ans<<'\n';
    return 0;
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注