The area of the union of circles

原来的求圆交函数有问题，交上去WA，用了个别人的

#include<stdio.h>
#include<math.h>
#include<iostream>
#include<algorithm>
using namespace std;
const double eps=1e-8;
const double PI=acos(-1.0);
struct Circle{
    double x,y,r,angle;
    int d;
    Circle(){}
    Circle(double xx,double yy){x=xx;y=yy;}
};
struct Node{//入点+1，可以表示覆盖了多少层
    double angle;
    int flag;
    Node(){}
    Node(double a,int b){angle=a;flag=b;}
};
Circle c[1005];
Node p[2010];
int n;
double dist(Circle a,Circle b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int cmp(Circle a,Circle b){
    return a.r<b.r;
}
void init(){
    sort(c,c+n,cmp);
    int i,j;
    for(i=0;i<n;i++){//计算每个圆被覆盖的次数 
    	//c[i].d=1;
    	for(j=i+1;j<n;j++){
    		//if(c[i].x==c[j].x && c[i].y==c[j].y ||c[i].r<eps) continue;
    		if(c[j].r-c[i].r-dist(c[i],c[j])>-eps) c[i].d++;
    	}   
    }
    
}
int getarc(Circle p,Circle q,double &a,double &b){//a是小角，b是大角
    double d=dist(p,q);
	if(p.r + q.r - d < 0) return 0;//外离 
  	if(fabs(p.r-q.r) - d >-eps) return 0;//内含及内切 
    double tmp1=acos((d*d+p.r*p.r-q.r*q.r)/(2*d*p.r));
    double tmp2=atan2(q.y-p.y,q.x-p.x);
    //printf("%f %f %f %f\n",p.x,p.y,q.x,q.y);
    a=tmp2-tmp1;    b=tmp1+tmp2;
    if(b>PI) b-=2*PI;
    if(a<-PI) a+=2*PI;
    return 2;
}
int CirCrossCir(Circle p1,Circle p2,Circle &cp1,Circle &cp2) {//两圆求交点数 ,并用cp1和cp2存储其交点
    double mx = p2.x - p1.x, sx = p2.x + p1.x, mx2 = mx * mx;
    double my = p2.y - p1.y, sy = p2.y + p1.y, my2 = my * my;
    double sq = mx2 + my2, d = -(sq - (p1.r - p2.r)*(p1.r - p2.r)) * (sq - (p1.r + p2.r)*(p1.r + p2.r));//sq为两圆距离的平方和
    if (d + eps < 0) return 0; if (d < eps) d = 0; else d = sqrt(d);
    double x = mx * ((p1.r + p2.r) * (p1.r - p2.r) + mx * sx) + sx * my2;
    double y = my * ((p1.r + p2.r) * (p1.r - p2.r) + my * sy) + sy * mx2;
    double dx = mx * d, dy = my * d; sq *= 2;
    cp1.x = (x - dy) / sq; cp1.y = (y + dx) / sq;
    cp2.x = (x + dy) / sq; cp2.y = (y - dx) / sq;
    if (d > eps) return 2; else return 1;
}
int cmp1(Node a,Node b){
    if(fabs(a.angle-b.angle)>eps) return a.angle<b.angle;
    else return a.flag>b.flag;
}
double cross(Circle a,Circle b){
    return a.x*b.y-a.y*b.x;
}
double getarea(double a,double r){
    return r*r*(a/2.0);
}
void solve(){
    double area=0;
    
    for(int i=0;i<n;i++){
        int cnt=0,flag=0;
        for(int j=0;j<n;j++){
            if(i==j) continue;
            Circle cp1,cp2;
            double a,b;
            //if(getarc(c[i],c[j],a,b)<2) continue;
            if(CirCrossCir(c[i],c[j],cp2,cp1)<2) continue;
            a = atan2(cp1.y - c[i].y, cp1.x - c[i].x);//求圆心指向交点的向量的极角，注意这里atan2(y,x)函数要先y后x
            b = atan2(cp2.y - c[i].y, cp2.x - c[i].x);
            //
            //printf("%f %f\n",a,b);
            p[cnt++]=Node(a,1); p[cnt++]=Node(b,-1);//a是小角，b是大角
            if(a-b>eps) flag++;//记录圆的最左点被覆盖了多少次
        }
        p[cnt++]=Node(-PI,flag);    p[cnt++]=Node(PI,-flag);
        sort(p,p+cnt,cmp1);
        int s=flag+c[i].d;//+1是因为每段圆弧至少被c[i]这个圆覆盖
        for(int j=1;j<cnt;j++){
            if(s==1) {//被覆盖了1次及以上的面积
                Circle a,b;
                a.x=c[i].x+c[i].r*cos(p[j-1].angle); a.y=c[i].y+c[i].r*sin(p[j-1].angle);
                b.x=c[i].x+c[i].r*cos(p[j].angle); b.y=c[i].y+c[i].r*sin(p[j].angle);
                double k=p[j].angle-p[j-1].angle;
                
                //printf("%f %f %f %f  :%f\n",a.x,a.y,b.x,b.y,cross(a,b)*0.5);
                area+=(k-sin(k))*c[i].r*c[i].r*0.5;//弓形面积
                area+=cross(a,b)*0.5;//有向面积
            }
            s+=p[j].flag;
        }
    }
    
    printf("%.5f\n",area);
}
int main(){
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
    #endif // ONLINE_JUDGE

    scanf("%d",&n);
    for(int i=0;i<n;i++) {
    	scanf("%lf%lf%lf",&c[i].x,&c[i].y,&c[i].r);
    	c[i].d=1;//至少被自己覆盖 
    }
    init();
    solve();
    return 0;
}