AOJ 1323: Encycling Circles

この記事は、AOJ-ICPC Advent Calendarの記事です。


円全体を考えるより、中心の動ける範囲を考えたほうがやりやすい。
中心が動ける範囲は、たくさんの円の積集合。これの境界線上で中心を動かして、長さをもとめていく。この集合上でとがった場所に中心があるときの計算が面倒だったりする。
3つの円が同じ点で交わるケースに注意。

#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
using namespace std;
const double EPS = 1e-10;
const double INF = 1e+10;
const double PI = acos(-1);

Pt c[110];
double r[110];
vector<Pt> p[110];
vector<pair<int,int> > tp;
int v[11100];
double ad[110][110];
int main(){
    int a;
    double b;
    while(scanf("%d%lf",&a,&b),a){
        for(int i=0;i<a;i++){
            double x,y;
            scanf("%lf%lf%lf",&x,&y,r+i);
            c[i]=Pt(x,y);
        }
        for(int i=0;i<11100;i++)v[i]=0;
        for(int i=0;i<a;i++)p[i].clear();
        tp.clear();
        bool dame=false;
        for(int i=0;i<a;i++)if(r[i]>b+EPS)dame=true;
        for(int i=0;i<a;i++){
            for(int j=i+1;j<a;j++){
                if(iCC(c[i],b-r[i],c[j],b-r[j])==0)dame=true;
            }
        }
        if(dame){printf("0.0\n");continue;}
        for(int i=0;i<a;i++)for(int j=0;j<a;j++){
            if(iCC(c[i],b-r[i],c[j],b-r[j])==1){
                pair<Pt,Pt> t=pCC(c[i],b-r[i],c[j],b-r[j]);
                p[i].push_back(t.first);
                p[i].push_back(t.second);
            }
        }
        for(int i=0;i<a;i++){
            int n=p[i].size();
            if(n==0){
                p[i].push_back(c[i]+Pt(cos(1),sin(1)));
                p[i].push_back(c[i]+Pt(cos(1),sin(1)));
            }else{
                for(int j=0;j<n;j++){
                    for(int k=0;k<n-1;k++){
                        if((p[i][k]-c[i]).arg()>(p[i][k+1]-c[i]).arg()){
                            swap(p[i][k],p[i][k+1]);
                        }
                    }
                }
                p[i].push_back(p[i][0]);
            }
        }
        double ret=0;
        for(int i=0;i<a;i++)for(int j=0;j<p[i].size()-1;j++){
            if(p[i].size()>2&&(p[i][j]-p[i][j+1]).ABS()<EPS)continue;
            double th=((p[i][j]-c[i]).arg()+(p[i][j+1]-c[i]).arg()+((j==p[i].size()-2)?PI*2:0))/2;
            Pt t=c[i]+Pt(cos(th),sin(th))*(b-r[i]);
            bool ok=true;
            for(int k=0;k<a;k++){
                if((t-c[k]).ABS()>b-r[k]+EPS){ok=false;break;}
            }
            if(ok){
            //  printf("%f %f %f %f\n",p[i][j].x,p[i][j].y,p[i][j+1].x,p[i][j+1].y);
                ret+=(b*2-r[i])*(((j==p[i].size()-2)?PI*2:0)+(p[i][j+1]-c[i]).arg()-(p[i][j]-c[i]).arg());
                tp.push_back(make_pair(i,j));
            }
        }
        int sz=tp.size();
        int at=0;
        int prev=-1;
        for(int i=0;i<a;i++)for(int j=0;j<a;j++){
            if(iCC(c[i],b-r[i],c[j],b-r[j])==1){
                pair<Pt,Pt> t=pCC(c[i],b-r[i],c[j],b-r[j]);
                double th=(t.first-c[j]).arg()-(t.first-c[i]).arg();
                if(th<0)th+=2*PI;
                if(th>PI)th=2*PI-th;
                ad[i][j]=th*b;
//              printf("%d %d: %f\n",i,j,ad[i][j]);
            }
        }
        for(int i=0;i<sz;i++){
            for(int j=0;j<sz;j++){
                if(at==j||(prev==j&&sz>2))continue;
                bool ok=false;
                if((p[tp[at].first][tp[at].second]-p[tp[j].first][tp[j].second]).ABS()<EPS)ok=true;
                if((p[tp[at].first][tp[at].second+1]-p[tp[j].first][tp[j].second]).ABS()<EPS)ok=true;
                if((p[tp[at].first][tp[at].second]-p[tp[j].first][tp[j].second+1]).ABS()<EPS)ok=true;
                if((p[tp[at].first][tp[at].second+1]-p[tp[j].first][tp[j].second+1]).ABS()<EPS)ok=true;
                if(ok){
                    ret+=ad[tp[at].first][tp[j].first];
                //  printf("%d %d\n",tp[at].first,tp[j].first);
            //      printf("%f\n",ret);
                    prev=at;at=j;break;
                }
            }
        }
        printf("%.12f\n",ret);
    }
}