SPOJ - CIRU The area of the union of circles （圆形面积并）

题意：

给n个圆的圆心与半径，求这些圆的面积并。

思路：

两个办法，
一个是很优美的几何做法：http://www.cnblogs.com/oyking/p/3424999.html，

一个是用自适应Simpson积分来来求圆并：http://blog.csdn.net/qpswwww/article/details/44201333 。

代码：

//几何方法
#include<bits/stdc++.h>
using namespace std;

const int inf = 0x3f3f3f3f;
const double eps = 1e-8 ;
const double pi = acos(-1.0);

inline double sqr(double x) {
    return x * x ;
}
inline int sgn(double x) {
    if (fabs(x)  < eps) return 0 ;
    return x > 0? 1 : -1 ;
}

struct Point{
    double x,y ;
    Point(){}
    Point(double _x, double _y): x(_x), y(_y){}
    double norm() { return sqrt(sqr(x) + sqr(y)) ; }
    void input() { scanf("%lf%lf", &x, &y) ; }

    friend Point operator + (Point a, Point b) { return Point(a.x + b.x, a.y + b.y) ; }
    friend Point operator - (Point a, Point b) { return Point(a.x - b.x, a.y - b.y) ; }
    friend Point operator * (Point a, double b) { return Point(a.x * b, a.y * b) ; }
    friend Point operator / (Point a, double b) { return Point(a.x / b, a.y / b) ; }
    friend bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0 ; }
    friend bool operator < (Point a, Point b) { return sgn(a.x - b.x) < 0 || (sgn(a.x - b.x) == 0 && sgn(a.y - b.y) < 0) ; }

    Point rotate(Point p, double ang) {
        Point v = (*this) - p ;
        Point ret ;
        ret.x = v.x * cos(ang) - v.y * sin(ang) ;
        ret.y = v.y * cos(ang) + v.x * sin(ang) ;
        return ret + p ;
    }
};
double det(Point a, Point b) { return a.x * b.y - b.x * a.y; }
int CircleInterCircle(Point c1, double r1, Point c2, double r2, Point *res)
{
    double d = (c1 - c2).norm() ;
    if (sgn(d) == 0) {
        if (sgn(r1 - r2) == 0) return -1 ;
        return 0 ;
    }
    if (sgn(r1 + r2 - d) < 0) return 0 ;
    if (sgn(fabs(r1 - r2) - d) > 0) return 0 ;

    double ang = atan2((c2 - c1).y, (c2 - c1).x) ;
    double vang = acos( (sqr(r1) + sqr(d) - sqr(r2)) / (2*r1*d) ) ;
    res[0] = Point(c1.x+r1, c1.y).rotate(c1, ang + vang) ;
    res[1] = Point(c1.x+r1, c1.y).rotate(c1, ang - vang) ;
    if (res[0] == res[1]) return 1 ;
    return 2 ;
}
struct Region {
    double st, ed;
    Region(){}
    Region(double _st, double _ed): st(_st), ed(_ed) {}
    bool operator < (const Region &a)const {
        return sgn(st - a.st) < 0 || (sgn(st - a.st) == 0 && sgn(ed - a.ed) < 0) ;
    }
};
struct Circle {
    Point c ;
    double r ;
    vector<Region> reg ;
    Circle(){}
    Circle(Point _c, double _r): c(_c), r(_r) {}

    void add(const Region &r) { reg.push_back(r); }
    double area(double ang = pi) { return ang * sqr(r); }
    Point makepoint(double ang) { return Point(c.x + r*cos(ang) , c.y + r*sin(ang)); }

    bool operator < (const Circle &a)const {
        return sgn(r - a.r) < 0 || (sgn(r - a.r) == 0 && c < a.c) ;
    }
    bool operator == (const Circle &a)const {
        return sgn(r - a.r) == 0 && c == a.c ;
    }

};

double Angle(Point p) { return atan2(p.y, p.x); }
double area_CircleUnion(Circle *cir, int n)
{
    bool ok[n+5] ;
    memset(ok, true, sizeof ok) ;
    double ans = 0 ;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) if (ok[j]) {
            if (i == j) continue ;
            double d = (cir[i].c - cir[j].c).norm() ;
            if (sgn(d + cir[i].r - cir[j].r) <= 0) { //!!!!
                ok[i] = false ;
                break ;
            }
        }
    }
    for (int i = 0; i < n; i++) if (ok[i]) {
        Point p[2] ;
        bool flag = false ;
        for (int j = 0; j < n; j++) if(ok[j]) {
            if (i == j) continue ;
            int k = CircleInterCircle(cir[i].c, cir[i].r, cir[j].c, cir[j].r, p) ;
            if (k != 2) continue ;
            flag = true ;

            double ang1 = Angle(p[1] - cir[i].c), ang2 = Angle(p[0] - cir[i].c) ;
            if (sgn(ang1) < 0) ang1 += 2*pi;
            if (sgn(ang2) < 0) ang2 += 2*pi;
            if (sgn(ang1 - ang2) > 0) cir[i].add(Region(ang1, 2*pi)), cir[i].add(Region(0, ang2)) ;
            else cir[i].add(Region(ang1, ang2)) ;
        }
        if (!flag) {
            ans += cir[i].area() ;
            continue ;
        }
        sort(cir[i].reg.begin(), cir[i].reg.end()) ;
        int cnt = 1 ;
        for (int j = 1; j < int(cir[i].reg.size()); j++) {
            if (sgn(cir[i].reg[cnt - 1].ed - cir[i].reg[j].st) >= 0) {
                cir[i].reg[cnt - 1].ed = max(cir[i].reg[cnt - 1].ed, cir[i].reg[j].ed) ;
            }
            else {
                cir[i].reg[cnt++] = cir[i].reg[j] ;
            }
        }
        cir[i].add(Region()) ;
        cir[i].reg[cnt] = cir[i].reg[0] ;
        for (int j = 0; j < cnt; j++) {
            p[0] = cir[i].makepoint(cir[i].reg[j].ed) ;
            p[1] = cir[i].makepoint(cir[i].reg[j+1].st) ;
            ans += det(p[0], p[1]) / 2.0 ;
            double ang = cir[i].reg[j + 1].st - cir[i].reg[j].ed ;
            if (sgn(ang) < 0) ang += 2*pi ;
            ans += 0.5 * sqr(cir[i].r) * (ang - sin(ang)) ;
        }
    }
    return ans ;
}
int n ;
Circle cir[1010] ;
int main()
{
    while (cin >> n) {
        for (int i = 0; i < n; i++) {
            cir[i].c.input() ; scanf("%lf", &cir[i].r) ;
        }
        printf("%.3f\n", area_CircleUnion(cir, n) + eps) ;
    }
    return 0;
}