LA 2572 - Viva Confetti

這道題是訓練指南上面的例題，判斷小圓弧的時候，求出一個圓跟其它所有圓的交點，然後判斷兩個相鄰交點之間的中點。如果這個中點沒有在這個圓上面的所有圓裡面，那麼這個交點就沒有被擋住，也就是說這個圓可見。然後這個交點下面的第一個圓也是可見的。

做的時候看了一下别人的部落格，然後自作聰明的把判斷兩個圓相交的函數改了。本來的函數裡面是求圓C1跟C2的交點，求完焦點以後隻把容器C1裡面加入交點的極角，這樣循環的時候就需要  

for(int i=0;i<n;++i){

      for(int t=0;t<n;++t){

           CircleCircleIntersection(i,t);

      }

    }

我想把它改成

for(int i=0;i<n;++i){

     for(int t=i+1;t<n;++t){

          CircleCircleIntersection(i,t);

      }

    }

這樣在圓相交的函數裡面一次就需要修改C1跟C2的容器。這樣就一直錯。開始沒看出來什麼函數錯了，一直在找，後來發現是求C2的極角的時候出問題了，修改過來就對了。但是還是不能一次把兩個圓的交點一起求。

找了好長時間錯誤，終于找到錯那了。。我加的那些計算每錯，錯的地方是精确度。題目要求精确度是加減五乘十的負十三次方，我把eps設定成了十的負十三次放。把精确度改過來以後，怎麼寫怎麼對。。

代碼如下：

#include<iostream>
#include<cmath>
#include<cstring>
#include<vector>
#include<algorithm>
using namespace std;
const int maxn=105;
const double pi=cos(-1);

struct Point{
    double x,y;
    Point(double x,double y):x(x),y(y){}
};

typedef Point Vector;

Vector operator + (Vector A,Vector B){return Vector(A.x+B.x,A.y+B.y);}

Vector operator - (Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}

Vector operator * (Vector a,double p){return Vector(a.x*p,a.y*p);}

Vector operator / (Vector a,double p){return Vector(a.x/p,a.y/p);}

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

const double eps=1e-14;
int dcmp(double x){
    if(fabs(x)<eps)return 0;
    return x<0?-1:1;
}

bool operator == (const Point& a ,Point& b){
    return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}

struct Circle{
    Point c;
    double r;
    Circle(Point c=Point(0,0),double r=0):c(c),r(r){}
    Point point(double a){
        return Point(c.x+cos(a)*r, c.y+sin(a)*r);
    }
};

double Dot(Vector a,Vector b){
    return a.x*b.x+a.y*b.y;
}

double Length(Vector A){
    return sqrt(Dot(A,A));
}

double Angle(Vector a,Vector b){//ÏòÁ¿¼Ð½Ç 
    return acos(Dot(a,b)/Length(a)/Length(b));
}

double angle(Vector a){//ÏòÁ¿µÄ¼«½Ç 
    return atan2(a.y,a.x);
}
int n;
vector<double> pointAng[maxn];
Circle C[maxn];
int vis[maxn];

bool PointInCircle(Point p,Circle C){
    if(dcmp(Length(p-C.c)-C.r)<0)return 1;
    return 0;
}

bool CircleInCircle(int c1,int c2){
    Circle C1=C[c1];
    Circle C2=C[c2];
    double d=Length(C1.c-C2.c);
    if(dcmp(d)==0){
        if(dcmp(C1.r-C2.r)<=0)return 1;
        return 0;
    }
    if(dcmp(C2.r-C1.r-d)>0)return 1;
    return 0;
}

int CircleCircleIntersection(int c1,int c2){
//    cout<<"c1="<<c1<<' '<<"c2="<<c2<<endl;
    Circle C1=C[c1];
    Circle C2=C[c2];
    double d=Length(C1.c-C2.c);
    if(dcmp(d)==0){
        if(dcmp(C1.r-C2.r)==0)return -1;
        return 0;
    }
    if(dcmp(C1.r+C2.r-d)<0)return 0;
    if(dcmp(fabs(C1.r-C2.r)-d)>0)return 0;
    double a=angle(C2.c-C1.c);
    double da=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));
//    cout<<"///"<<C1.r<<' '<<C2.r<<' '<<d<<' '<<(C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d)<<endl;
    Point p1 =C1.point(a-da),p2=C1.point(a+da);
//    cout<<c1<<' '<<c2<<' '<<"C1 "<<C1.c.x<<' '<<C1.c.y<<' '<<C1.r<<' '<<C2.c.x<<' '<<C2.c.y<<' '<<C2.r<<" zz "<<p1.x<<' '<<p1.y<<' '<<p2.x<<' '<<p2.y<<"a="<<a<<' '<<"da="<<da<<endl;
    if(p1==p2)return 1;
    pointAng[c1].push_back(a-da);
    pointAng[c1].push_back(a+da);
//    double a2=pi+a;
//    if(a2<-pi)a2+=(2*pi);
//    if(a2>pi)a2-=(2*pi);
//    da=Angle(p1-C2.c,C1.c-C2.c);
    double pa1=angle(p1-C2.c);
    double pa2=angle(p2-C2.c);
    pointAng[c2].push_back(pa1);
    pointAng[c2].push_back(pa2);
//    a=angle(C1.c-C2.c);
//    da=acos((C2.r*C2.r+d*d-C1.r*C1.r)/(2*C2.r*d));
//    pointAng[c2].push_back(a-da);
//    pointAng[c2].push_back(a+da);
    return 2;
} 

void solve(){
    memset(vis,0,sizeof(vis));
    for(int i=0;i<n;++i){
        for(int t=i+1;t<n;++t){
//            cout<<"i="<<i<<' '<<"t="<<t<<endl;
            CircleCircleIntersection(t,i);
        }
    }
    for(int i=0;i<n;++i){
        if(!pointAng[i].size())continue;
        sort(pointAng[i].begin(),pointAng[i].end());
        vector<double>::iterator iter=unique(pointAng[i].begin(),pointAng[i].end());
        pointAng[i].resize(iter-pointAng[i].begin());
    }
    for(int i=0;i<n;++i){
        int ez=pointAng[i].size();
        if(!ez){
            bool can=1;
            for(int t=i+1;t<n;++t){
                if(CircleInCircle(i,t)){can=0;break;}
            }
            if(can)vis[i]=1;
        }
        else{
            pointAng[i].push_back(pointAng[i][0]);
            for(int j=0;j<pointAng[i].size();++j){
                Point midpoint=C[i].point((pointAng[i][j]+pointAng[i][j+1])/2);
                bool can=1;
                for(int k=i+1;k<n;++k){
                    if(PointInCircle(midpoint,C[k])){can=0;break;}
                }
    //            cout<<i<<' '<<j<<' '<<pointAng[i].size()<<' '<<midpoint.x<<' '<<midpoint.y<<' '<<(midpoint.x-C[i].c.x)*(midpoint.x-C[i].c.x)+(midpoint.y-C[i].c.y)*(midpoint.y-C[i].c.y)<<' '<<can<<endl;
                if(can){
                    vis[i]=1;
                    for(int k=i-1;k>=0;--k){
                        if(PointInCircle(midpoint,C[k])){
//                            cout<<"tttttttttt"<<endl;
                            vis[k]=1;break;}
                    }
                }
            }        
        }
    }
    int ans=0;
    for(int i=0;i<n;++i){
        if(vis[i])ans++;
    }
    cout<<ans<<endl;
}

int main(){
//    freopen("data.txt","r",stdin);
    ios::sync_with_stdio(false);
    while(cin>>n){
        if(n==0)break;
        for(int i=0;i<n;++i){
            cin>>C[i].c.x>>C[i].c.y>>C[i].r;
        }
        solve();
    }
    return 0;
}