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;}