多圆面积交、并

看到个板子，觉得不错。

 typedef long long LL;  
typedef unsigned long long ULL;  
typedef vector <int> VI;  
const int INF = 0x3f3f3f3f;  
const double eps = 1e-10;  
const int MOD = 100000007;  
const int MAXN = 1000010;  
const double PI = acos(-1.0);  
#define sqr(x) ((x)*(x))  
const int N = 1010;  
double area[N];  
int n;  
  
int dcmp(double x)  
{  
    if (x < -eps) return -1;  
    else return x > eps;  
}  
  
struct cp  
{  
    double x, y, r, angle;  
    int d;  
    cp() {}  
    cp(double xx, double yy, double ang = 0, int t = 0)  
    {  
        x = xx;  
        y = yy;  
        angle = ang;  
        d = t;  
    }  
    void get()  
    {  
        scanf("%lf%lf%lf", &x, &y, &r);  
        d = 1;  
    }  
} cir[N], tp[N * 2];  
  
double dis(cp a, cp b)  
{  
    return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));  
}  
  
double cross(cp p0, cp p1, cp p2)  
{  
    return (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x);  
}  
  
int CirCrossCir(cp p1, double r1, cp p2, double r2, cp &cp1, cp &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 - sqr(r1 - r2)) * (sq - sqr(r1 + r2));  
    if (d + eps < 0) return 0;  
    if (d < eps) d = 0;  
    else d = sqrt(d);  
    double x = mx * ((r1 + r2) * (r1 - r2) + mx * sx) + sx * my2;  
    double y = my * ((r1 + r2) * (r1 - r2) + 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;  
}  
  
bool circmp(const cp& u, const cp& v)  
{  
    return dcmp(u.r - v.r) < 0;  
}  
  
bool cmp(const cp& u, const cp& v)  
{  
    if (dcmp(u.angle - v.angle)) return u.angle < v.angle;  
    return u.d > v.d;  
}  
  
double calc(cp cir, cp cp1, cp cp2)  
{  
    double ans = (cp2.angle - cp1.angle) * sqr(cir.r)  
                 - cross(cir, cp1, cp2) + cross(cp(0, 0), cp1, cp2);  
    return ans / 2;  
}  
  
void CirUnion(cp cir[], int n)  
{  
    cp cp1, cp2;  
    sort(cir, cir + n, circmp);  
    for (int i = 0; i < n; ++i)  
        for (int j = i + 1; j < n; ++j)  
            if (dcmp(dis(cir[i], cir[j]) + cir[i].r - cir[j].r) <= 0)  
                cir[i].d++;  
    for (int i = 0; i < n; ++i)  
    {  
        int tn = 0, cnt = 0;  
        for (int j = 0; j < n; ++j)  
        {  
            if (i == j) continue;  
            if (CirCrossCir(cir[i], cir[i].r, cir[j], cir[j].r,  
                            cp2, cp1) < 2) continue;  
            cp1.angle = atan2(cp1.y - cir[i].y, cp1.x - cir[i].x);  
            cp2.angle = atan2(cp2.y - cir[i].y, cp2.x - cir[i].x);  
            cp1.d = 1;  
            tp[tn++] = cp1;  
            cp2.d = -1;  
            tp[tn++] = cp2;  
            if (dcmp(cp1.angle - cp2.angle) > 0) cnt++;  
        }  
        tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, PI, -cnt);  
        tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, -PI, cnt);  
        sort(tp, tp + tn, cmp);  
        int p, s = cir[i].d + tp[0].d;  
        for (int j = 1; j < tn; ++j)  
        {  
            p = s;  
            s += tp[j].d;  
            area[p] += calc(cir[i], tp[j - 1], tp[j]);  
        }  
    }  
}  
  
void solve()  
{  
    scanf("%d", &n);  
    for (int i = 0; i < n; ++i)  
        cir[i].get();  
    memset(area, 0, sizeof(area));  
    CirUnion(cir, n);  
    //去掉重复计算的  
    for (int i = 1; i <= n; ++i)  
    {  
        area[i] -= area[i + 1];  
    }  
    //area[i]为重叠了i次的面积  
    //tot 为总面积  
    double tot = 0;  
    for(int i=1; i<=n; i++) tot += area[i];  
    printf("%f\n", tot);  
}  
  
int main()  
{  
    //freopen("input.txt", "r", stdin);  
    return 0;  
}  

 多圆面积并

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#define ld double
#define eps 1e-13
using namespace std;
int n,top,st,ed;
ld xl[1001],xr[1001];
ld ans;
bool del[1001];
struct data{ld x,y,r;}t[1001],sk[1001];
struct line{ld l,r;}p[1001];
ld dis(data a,data b)
{return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
bool cmp1(data a,data b){return a.r<b.r;}
bool cmp2(data a,data b){return a.x-a.r<b.x-b.r;}
bool cmp3(line a,line b){return a.l<b.l;}
void ini()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
       {scanf("%lf%lf%lf",&t[i].x,&t[i].y,&t[i].r);}
   sort(t+1,t+n+1,cmp1);
   for(int i=1;i<=n;i++)
      for(int j=i+1;j<=n;j++)
         if(dis(t[i],t[j])<=t[j].r-t[i].r)
               {del[i]=1;break;}
   for(int i=1;i<=n;i++)if(!del[i])sk[++top]=t[i];n=top;
   sort(sk+1,sk+n+1,cmp2);
}
ld getf(ld x)
{
    int sz=0,i,j;ld r,len=0,dis;
    for(i=st;i<=ed;i++)
   {
       if(x<=xl[i]||x>=xr[i])continue;
       dis=sqrt(sk[i].r-(x-sk[i].x)*(x-sk[i].x));
       p[++sz].l=sk[i].y-dis;p[sz].r=sk[i].y+dis;
   }
    sort(p+1,p+sz+1,cmp3);
    for(i=1;i<=sz;i++)
    {
        r=p[i].r;
        for(j=i+1;j<=sz;j++)
        {
            if(p[j].l>r)break;
            if(r<p[j].r)r=p[j].r;
        }
        len+=r-p[i].l;i=j-1;
    }
    return len;
}
ld cal(ld l,ld fl,ld fmid,ld fr)
{return (fl+fmid*4+fr)*l/6;}
ld simpson(ld l,ld mid,ld r,ld fl,ld fmid,ld fr,ld s)
{
    ld m1=(l+mid)/2,m2=(r+mid)/2;
    ld f1=getf(m1),f2=getf(m2);
    ld g1=cal(mid-l,fl,f1,fmid),g2=cal(r-mid,fmid,f2,fr);
    if(fabs(g1+g2-s)<eps)return g1+g2;
   return simpson(l,m1,mid,fl,f1,fmid,g1)+simpson(mid,m2,r,fmid,f2,fr,g2);
}
void work()
{
    int i,j;ld l,r,mid,fl,fr,fmid;
    for(i=1;i<=n;i++){xl[i]=sk[i].x-sk[i].r;xr[i]=sk[i].x+sk[i].r;sk[i].r*=sk[i].r;}
    for(i=1;i<=n;i++)
    {
        l=xl[i];r=xr[i];
        for(j=i+1;j<=n;j++)
        {
            if(xl[j]>r)break;
         if(xr[j]>r)r=xr[j];
        }
        st=i;ed=j-1;i=j-1;
      mid=(l+r)/2;
      fl=getf(l);fr=getf(r);fmid=getf(mid);
      ans+=simpson(l,mid,r,fl,fmid,fr,cal(r-l,fl,fmid,fr));
    }
}
int main()
{
    ini();
    work();
    printf("%.3lf",ans);
    return 0;
}