计算几何模板

计算几何

通用

typedef double db;
const db eps = 1e-8;
const db inf = 1e20;
const db pi = acos(-1.0);
const int maxp = 1010;
// Compares a double to zero
int sgn(db x){
    if(fabs(x) < eps) return 0;
    if(x < 0) return -1;
    return 1;
}
//square of a double 
inline db sqr(db x){return x*x;}

一. 二维几何

1. 点

所有的点操作封装在(Point) 结构体里面，通过成员函数的方式调用相应功能

标号	成员函数	含义	备注
1	bool operator == (Point b)const	判断两个点是否相等
2	bool operator < (Point b)const	判断两个点的顺序， 按照((x,y)) 的关键字顺序
3	Point operator - (const Point &b)const	两个点相减，返回向量
4	db operator ^ (const Point &b)const	向量的叉积，返回有向面积
5	db operator * (const Point &b)const	向量的点积
6	db len()	向量长度
7	db len2()	向量长度平方
8	db distance(Point p)	距离
9	Point operator + (const Point &b)const	向量加法
10	Point operator * (const db &k) const	向量乘以标量
11	Point operator /(const db &k)const	向量除以标量
12	db rad(Point a, Point b)	返回 ((pa, pb)) 的夹角	4，5
13	Point trunc(db r)	将向量长度调整为 r	6
14	Point rotleft()	将向量左转 90 度
15	Point rotright()	将向量右转 90 度
16	Point rotate(Point p, db angle)	将当前点绕 p 逆时针转 angle 度

struct Point{
    db x, y;
    Point(){}
    Point(db x, db y):x(x), y(y){}
    void input(){
        scanf("%lf%lf",&x, &y);//如果为longdouble则需要改为Lf
    }
};

基本函数

// 1.比较两点是否相同
bool operator == (Point b)const {
    return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
}
// 2.按照(x,y)的优先级比较，因题目不同要修改
bool operator < (Point b)const{
    return sgn(x-b.x) == 0 ? sgn(y - b.y) < 0 : x < b.x;
}
// 3.常用，求向量
Point operator - (const Point &b)const{
    return Point(x - b.x, y - b.y);
}
// 4. 叉积
db operator ^ (const Point &b)const{
    return x * b.y - y * b.x;
}
// 5. 点积
db operator * (const Point &b)const{
    return x * b.x + y * b.y;
}
// 6. 返回长度
db len(){
    return hypot(x, y);
} 
// 7. 返回长度平方
db len2(){
    return x * x + y * y;
}
// 8. 返回两点距离 
db distance(Point p){
    return hypot(x - p.x, y - p.y);
}
// 9. 向量加
Point operator + (const Point &b)const{
    return Point(x + b.x, y + b.y);
}
// 10. 标量乘
Point operator * (const db &k) const {
    return Point(x * k,  y * k);
}
// 11. 标量除
Point operator /(const db &k)const {
    return Point(x / k, y / k);
}

1.1 计算(vec{pa}) 与 (vec{pb})夹角

// 12.
// 计算 pa 和 pb 的夹角
// 就是求这个点看 a，b所成的夹角
// LightOJ 1203
db rad(Point a, Point b){
    Point p = *this;
    return fabs(atan2( fabs((a-p)^(b-p)), (a-p) * (b-p) ));
}

1.2 转换长度

// 13. 化为长度为 r 的向量
Point trunc(db r){
    db l = len();
    if(!sgn(l)) return *this;
    r /= l;
    return Point(x*r, y*r);
}

1.3 旋转

// 14. 逆时针旋转90度
Point rotleft(){
    return Point(-y, x);
}
// 15. 顺时针转90度
Point rotright(){
    return Point(y, -x);
}
// 16. 绕点 p 逆时针旋转 angle
Point rotate(Point p, db angle){
    Point v = (*this) - p;
    db c = cos(angle), s = sin(angle);
    return Point(p.x + v.x * c - v.y * s, p.y + v.x * s + v.y * c);
}

2. 线

标号	成员函数	含义	备注
17	bool operator == (Line v)	判断两条线段是否相等	1
18	Line (Point p, db angle)	根据一个点和一个倾斜角确定直线
19	Line (db a, db b, db c)	根据(ax+by+c=0) 确定直线
20	void adjust()	调整线段两个端点	2
21	db length()	线段长度	8
22	db angle()	返回直线倾斜角([0,pi))
23	int relation(Point p)	返回点和直线的关系(1.在左侧，2.在右侧,3.在直线上)	3，4
24	bool pointonseg(Point p)	点在线段上的判断，包括端点	3，4，5
25	bool parallel(Line v)	两向量平行（重合）	3，4
26	int segcrosseg(Line v)	两线段相交判断（2. 规范相交，1. 非规范相交，0. 不相交）	3，4，5
27	int linecorssseg(Line v)	直线与线段相交判断（含义与26相同，v代表线段）	3，4
28	int linecrossline(Line v)	两直线的关系（0. 平行，1. 重合，2.相交）	23，25
29	Point crosspoint(Line v)	两直线的交点，要保证两直线不平行或重合	3，4，10
30	db dispointtoline(Point p)	点到直线的距离	3，4，21
31	db dispointtoseg(Point p)	点到线段的距离	3，5，8，30
32	db dissegtoseg(Line v)	线段到线段的距离（前提是两线段不相交，相交距离为0）	31
33	Point lineprog(Point p)	返回 p 在直线上的投影	3，5，7，9，10
34	Point symmetrypoint(Point p)	返回 p 关于直线的对称点	33

struct Line{
    Point s, e;
    Line(){}
    Line(Point s, Point e):s(s),e(e){}
    void input(){
        s.input();
        e.input();
    }
};

基本函数

// 17. 判断直线是否相等
bool operator == (Line v){
    return (s == v.s) && (e == v.e);
}
// 20. 调整直线两点顺序
void adjust(){
    if(e < s) swap(s, e);
}
// 21. 求直线长度
db length(){
    return s.distance(e);
}

2.1 根据一点和倾斜角确定直线

// 18. 根据一个点和倾斜角angle确定直线， 0 <= angle < pi
Line (Point p, db angle){
    s = p;
    if(sgn(angle - pi / 2) == 0){
        e = (s + Point(0, 1));
    }else{
        e = (s + Point(1, tan(angle)));
    }
}

2.2 (ax+by+c=0) 求直线

// 19. ax + by + c = 0
Line (db a, db b, db c){
    if(sgn(a) == 0){ //  y = -c / b
        s = Point(0, -c/b);
        e = Point(1, -c/b);
    } else if (sgn(b) == 0){ // x = -c / a
        s = Point(-c/a, 0);
        e = Point(-c/a, 1);
    }else{//(0, -c/b), (1, (-c-a)/b)
        s = Point(0, -c/b);
        e = Point(1, (-c-a)/b);
    }
}

2.3 求直线倾斜角度

返回值范围 ([0,pi))

// 22. 返回直线倾斜角 0 <= angle < pi
db angle(){
    db k = atan2(e.y - s.y, e.x - s.x);
    if(sgn(k) < 0) k += pi;
    if(sgn(k - pi) == 0) k -= pi;
    return k;
}

2.4 点和直线关系

/*
    	23. 
        点和直线的关系
        1 在左侧
        2 在右侧
        3 在直线上
*/
int relation(Point p){
    int c = sgn((p-s) ^ (e-s));
    if(c < 0) return 1;
    else if(c > 0) return 2;
    else return 3;
}

2.5 点在线段上的判断

// 24. 点在线段上的判断，包括端点 第二个判断改为小于表示不包括端点
bool pointonseg(Point p){
    return sgn((p-s)^(e-s)) == 0 && sgn((p-s) * (p-e)) <= 0;
}

2.6 判断两向量平行

// 25. 两向量平行（对应直线平行或重合）
bool parallel(Line v){
    return sgn((e - s) ^ (v.e - v.s)) == 0;
}

2.7 线段相交

/*
    	26. 
        两线段相交判断
        2 规范相交
        1 非规范相交
        0 不相交
*/
int segcrosseg(Line v){
    int d1 = sgn((e - s) ^ (v.s - s));  //v.s 在 线段的哪一侧
    int d2 = sgn((e - s) ^ (v.e - s));  //v.e 在 线段的哪一侧
    int d3 = sgn((v.e - v.s) ^ (s - v.s));  
    int d4 = sgn((v.e - v.s) ^ (e - v.s));
    if((d1^d2) == -2 && (d3^d4) == -2) return 2; // 跨立实验满足 一个是-1一个是1

    // 1. v.s在线段上 || v.e 在线段上 || s 在另外一条线段上 || e在另外一条线段上
    return (d1 == 0 &&  sgn((v.s - s) * (v.s - e)) <= 0) ||
        (d2 == 0 && sgn((v.e - s) * (v.e - e)) <= 0) ||
        (d3 == 0 && sgn((s - v.s) * (s - v.e)) <= 0) || 
        (d4 == 0 && sgn((e - v.s) * (e - v.e)) <= 0);
}

2.8 直线与线段相交判断

/*
    	27. 
        直线与线段相交判断
        *this line -v seg 
        2 规范相交
        1 非规范相交
        0 不相交
*/
int linecorssseg(Line v){ // v是线段
    int d1 = sgn((e - s) ^ (v.s - s));
    int d2 = sgn((e - s) ^ (v.e - s));
    if((d1 ^ d2) == -2) return 2;
    return (d1 == 0 || d2 == 0);
}

2.9 两直线关系

/*
    	28. 
        两直线关系
        0 平行
        1 重合
        2 相交
*/
int linecrossline(Line v){
    // 如果平行，则看点是否在直线上
    if((*this).parallel(v)) return v.relation(s) == 3;
    return 2;
}

2.10 求两直线交点

/*
    	29. 
        求两直线的交点
        要保证两直线不平行或重合
    */
Point crosspoint(Line v){
    db a1 = (v.e - v.s) ^ (s - v.s);
    db a2 = (v.e - v.s) ^ (e - v.s);
    return Point((s.x * a2 - e.x * a1) / (a2 - a1), (s.y * a2 - e.y * a1) / (a2 - a1));
}

2.11 点到直线距离

// 30. 点到直线的距离
db dispointtoline(Point p){
    return fabs((p-s) ^ (e-s)) / length();
}

2.12 点到线段距离

// 31. 点到线段的距离
db dispointtoseg(Point p){
    if(sgn((p-s)*(e-s)) < 0 || sgn((p-e) * (s-e)) < 0)
        return min(p.distance(s), p.distance(e));
    return dispointtoline(p);
} 

2.13 线段与线段距离

/*
    	32. 
        返回线段到线段的距离
        前提是两线段不相交，相交距离就是 0 了
*/
db dissegtoseg(Line v){
    return min(min(dispointtoseg(v.s), dispointtoseg(v.e)), min(v.dispointtoseg(s), v.dispointtoseg(e)));
}

2.14 点在直线上的投影

/*
        33. 返回 p 在直线上的投影
*/
Point lineprog(Point p){
    return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );
}

2.15 关于直线的对称点

//  34. 返回点 p 关于直线的对称点
Point symmetrypoint(Point p){
    Point q = lineprog(p);
    return Point(2*q.x - p.x, 2 * q.y - p.y);
}

3. 圆

struct circle{
    Point p;
    db r;
    circle(){}
    circle(Point p, db r):p(p), r(r){}
    void input(){
        p.input();
        scanf("%lf", &r);
    }
    // 通过圆心角确定圆上的一个点
    Point point(double a){
        return Point(p.x + cos(a) * r, p.y + sin(a) * r);
    }
}

基本函数

bool operator == (circle v){
    return (p == v.p) && sgn(r - v.r) == 0;
}
bool operator < (circle v)const{
    return ((p<v.p) || (p == v.p) && sgn(r - v.r) < 0);
}
// 面积
db area(){
    return pi * r * r;
}
// 周长
db circumference(){
    return 2 * pi * r;
}

3.1 三角形外接圆

/*
        三角形的外接圆
        需要Point 的 + / rotate() 以及 Line 的crosspoint()
        利用两条边的中垂线得到圆心
        UVA 12304
*/
circle(Point a, Point b, Point c){
    Line u = Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));
    Line v = Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));
    p = u.crosspoint(v);
    r = p.distance(a);
}

3.2 三角形内切圆

/*
        三角形的内切圆
        bool t 没有作用，只是为了和上面外接圆函数区别
        UVA 12304
*/
circle(Point a, Point b, Point c, bool t){
    Line u, v;
    // u 为角 a 的平分线
    db m  = atan2(b.y-a.y, b.x-a.x), n = atan2(c.y - a.y, c.x - a.x);
    u.s = a;
    u.e = u.s + Point(cos((n+m)/2), sin((n+m)/2));
    // v 为角 b 的平分线
    m = atan2(a.y-b.y, a.x-b.x), n = atan2(c.y-b.y, c.x-b.x);
    v.s = b;
    v.e = v.s + Point(cos((n+m)/2), sin((n+m)/2));
    p = u.crosspoint(v);
    r = Line(a,b).dispointtoseg(p);
}

3.3 点和圆的关系

 /*
        点和圆的关系
        0 圆外
        1 圆上
        2 圆内
*/
int relation(Point b){
    db dst = b.distance(p);
    if(sgn(dst - r) < 0) return 2;
    else if(sgn(dst - r) == 0) return 1;
    return 0;
}

3.4 线段和圆, 直线和圆的关系

/*
        线段和圆的关系
        比较的是圆心到线段的距离和半径的关系
        2 交
        1 切
        0 不交
*/
int relationseg(Line v){
    db dst = v.dispointtoseg(p);
    if(sgn(dst - r) < 0) return 2;
    else if(sgn(dst - r) == 0) return 1;
    return 0;
}
int relationline(Line v){
    db dst = v.dispointtoline(p);
    if(sgn(dst - r) < 0) return 2;
    else if(sgn(dst - r) == 0) return 1;
    return 0;
}

3.5 两圆的关系

/*
        两圆的关系
        5 相离
        4 外切
        3 相交
        2 内切
        1 内含
*/
int relationcircle(circle v){
    db d = p.distance(v.p);
    if(sgn(d - r - v.r) > 0) return 5;
    if(sgn(d - r - v.r) == 0) return 4;
    db l = fabs(r - v.r);
    if(sgn(d - r - v.r) < 0 && sgn(d - l) > 0) return 3;
    if(sgn(d - l) == 0) return 2;
    if(sgn(d - l) < 0) return 1;
}

3.6 求两个圆的交点

/*
        求两个圆的交点，返回0表示没有交点，返回1是一个交点，2是两个交点
*/
int pointcrosscircle(circle v, Point &p1, Point &p2){
    int rel = relationcircle(v);
    if(rel == 1 || rel == 5) return 0;
    db d = p.distance(v.p);
    db l = (d * d + r * r - v.r * v.r) / (2 * d);
    db h = sqrt(r * r - l * l);
    Point tmp = p + (v.p - p).trunc(l);
    p1 = tmp + ((v.p - p).rotleft().trunc(h));
    p2 = tmp + ((v.p - p).rotright().trunc(h));
    if(rel == 2 || rel == 4)return 1;
    return 2;
}

3.7 求直线与圆的交点，返回交点个数

// 求直线与圆的交点，返回交点个数
int pointcrossline(Line v, Point &p1, Point &p2){
    if(!(*this).relationline(v)) return 0;
    Point a = v.lineprog(p);
    db d = v.dispointtoline(p);
    d = sqrt(r * r - d * d);
    if(sgn(d) == 0){
        p1 = a;
        p2 = a;
        return 1;
    }
    p1 = a + (v.e - v.s).trunc(d);
    p2 = a - (v.e - v.s).trunc(d);
    return 2;
}

3.8 通过固定两点，半径为(r1)的圆

// 得到 通过a，b两点，半径为r1的两个圆
int getcircle(Point a, Point b, db r1, circle &c1, circle &c2){
    circle x(a, r1), y(b, r1);
    int t = x.pointcrosscircle(y, c1.p, c2.p);
    if(!t) return 0;
    c1.r = c2.r = r;
    return t;
}

3.9 得到与直线 (u) 相切，过点 (q)，半径为 (r1) 的圆

// 得到与直线 u 相切，过点 q， 半径为 r1 的圆
int getcircle(Line u, Point q, db r1, circle &c1, circle &c2){
    db dis = u.dispointtoline(q);
    if(sgn(dis - r1 * 2) > 0) return 0;
    if(sgn(dis) == 0){
        c1.p = q + ((u.e - u.s).rotleft().trunc(r1));
        c2.p = q + ((u.e - u.s).rotright().trunc(r1));
        c1.r = c2.r = r1;
        return 2;
    }
    Line u1 = Line((u.s + (u.e - u.s).rotleft().trunc(r1)), (u.e + (u.e - u.s).rotleft().trunc(r1)));
    Line u2 = Line((u.s + (u.e - u.s).rotright().trunc(r1)), (u.e + (u.e - u.s).rotright().trunc(r1)));
    circle cc = circle(q, r1);
    Point p1, p2;
    if(!cc.pointcrossline(u1, p1, p2)) cc.pointcrossline(u2, p1, p2);
    c1 = circle(p1, r1);
    if(p1 == p2){
        c2 = c1;
        return 1;
    }
    c2 = circle(p2, r1);
    return 2;
}

3.10 同时与直线 (u),(v)相切，半径为(r1)的圆

// 同时与直线u，v相切，半径为r1的圆 , u,v不平行
int getcircle(Line u, Line v, db r1, circle &c1, circle &c2, circle &c3, circle &c4){
    if(u.parallel(v)) return 0;
    Line u1 = Line(u.s + (u.e - u.s).rotleft().trunc(r1), u.e + (u.e - u.s).rotleft().trunc(r1));
    Line u2 = Line(u.s + (u.e - u.s).rotright().trunc(r1), u.e + (u.e - u.s).rotright().trunc(r1));
    Line v1 = Line(v.s + (v.e - v.s).rotleft().trunc(r1), v.e + (v.e - v.s).rotright().trunc(r1));
    Line v2 = Line(v.s + (v.e - v.s).rotright().trunc(r1), v.e + (v.e - v.s).rotright().trunc(r1));

    c1.r = c2.r = c3.r = c4.r = r1;
    c1.p = u1.crosspoint(v1);
    c2.p = u1.crosspoint(v2);
    c3.p = u2.crosspoint(v1);
    c4.p = u2.crosspoint(v2);
    return 4;
}

3.11 同时与不相交圆 (cx),(cy)相切，半径为(r1)的圆

// 同时与不相交圆 cx, cy 相切，半径为r1的圆
int getcircle(circle cx, circle cy, db r1, circle &c1, circle &c2){
    circle x(cx.p, r1+cx.r), y(cy.p, r1+cy.r);
    int t = x.pointcrosscircle(y, c1.p, c2.p);
    if(!t) return 0;
    c1.r = c2.r = r1;
    return t;
}

3.12 过一点作圆的切线

// 过一点作圆的切线 (先判断点和圆的关系)
int tangentline(Point q, Line &u, Line &v){
    int x = relation(q);
    if(x == 2) return 0; //圆内
    if(x == 1){ //圆上
        u = Line(q, q+(q-p).rotleft());
        v = u;
        return 1;
    }
    db d = p.distance(q);
    db l = r * r / d;
    db h = sqrt(r * r - l * l);
    u = Line(q, p + ((q - p).trunc(l) + (q-p).rotleft().trunc(h)));
    v = Line(q, p + (q - p).trunc(l) + (q - p).rotright().trunc(h));
    return 2;
}

3.13 求两圆相交的面积

// 求两圆相交的面积
db areacircle(circle v){
    int rel = relationcircle(v);
    if(rel >= 4) return 0.0;
    if(rel <= 2) return min(area(), v.area()); //内部
    db d = p.distance(v.p);
    db hf = (r + v.r + d) / 2.0;
    db ss = 2 * sqrt(hf*(hf - r) * (hf - v.r) * (hf - d));
    db a1 = acos((r * r + d * d - v.r * v.r) / (2.0 * r * d));
    a1 = a1 * r * r;
    db a2 = acos((v.r * v.r + d * d - r * r) / (2.0 * v.r * d));
    a2 = a2 * v.r * v.r;
    return a1 + a2 - ss;
}

3.14 求圆和三角形( riangle pab) 的相交面积

// 求圆和三角形 pab 的相交面积
// POJ3675 HDU3982 HDU2892
db areatriangle(Point a, Point b){
    if(sgn((p-a)^(p-b)) == 0)return 0.0;
    Point q[5];
    int len = 0;
    q[len++] = a;
    Line l(a, b);
    Point p1, p2;
    if(pointcrossline(l, q[1], q[2]) == 2){
        if(sgn((a-q[1]) * (b-q[1])) < 0) q[len++] = q[1];
        if(sgn((a-q[2]) * (b-q[2])) < 0) q[len++] = q[2];
    }
    q[len++] = b;
    if(len == 4 && sgn((q[0] - q[1]) * (q[2] - q[1])) > 0) swap(q[1], q[2]);
    db res = 0;
    for(int i=0;i<len-1;i++){
        if(relation(q[i]) == 0 || relation(q[i+1]) == 0){
            db arg = p.rad(q[i], q[i+1]);
            res += r * r * arg / 2.0;
        }else{
            res += fabs((q[i] - p) ^ (q[i+1] - p)) / 2.0;
        }
    }
    return res;
}

3.15 求两个圆的公切线

• 非结构体成员函数

// a[i] 存放第 i 条公切线与 圆A 的交点
int getTangents(circle A, circle B, Point*a, Point *b){
    int cnt = 0;
    // 以A为半径更大的那个圆进行计算
    if(A.r < B.r) return getTangents(B, A, b, a);
    db d2 = (A.p-B.p).len2();  // 圆心距平方
    db rdiff = A.r - B.r;		// 半径差
    db rsum = A.r + B.r;		//半径和
    if(d2 < rdiff * rdiff) return 0; 	// 情况1，内含,没有公切线
    Vector AB = B.p - A.p;				// 向量AB，其模对应圆心距
    db base = atan2(AB.y, AB.x);		// 求出向量AB对应的极角
    if(d2 == 0 && A.r == B.r) return -1;// 情况3，两个圆重合，无限多切线
    if(d2 == rdiff * rdiff){ 			// 情况2，内切，有一条公切线
        a[cnt] = A.point(base);			
        b[cnt] = B.point(base);cnt++;
        return 1;
    }
    // 求外公切线
    db ang = acos((A.r - B.r) / sqrt(d2)); //求阿尔法
    // 两条外公切线
    a[cnt] = A.point(base+ang); b[cnt] = B.point(base+ang); cnt++;
    a[cnt] = A.point(base-ang); b[cnt] = B.point(base-ang); cnt++;
    if(d2 == rsum * rsum){  // 情况5，外切，if里面求出内公切线
        a[cnt] = A.point(base); b[cnt] = B.point(pi+base); cnt++;
    }
    else if(d2 > rsum * rsum){	//情况6，相离，再求出内公切线
        db ang = acos((A.r + B.r) / sqrt(d2));
        a[cnt] = A.point(base + ang); b[cnt] = B.point(pi+base+ang);cnt++;
        a[cnt] = A.point(base - ang); b[cnt] = B.point(pi+base-ang);cnt++;
    }
    // 此时，d2 < rsum * rsum 代表情况 4 只有两条外公切线
    return cnt;
}

3.16 求圆的并

struct circles{
    circle c[N];
    double ans[N];
    double pre[N];
    int n;
    circles(){}
    void add(circle cc){
        c[n++] = cc;
    }
    // x 包含在 y 中
    bool inner(circle x, circle y){
        if(x.relationcircle(y) != 1) return 0;
        return sgn(x.r - y.r) <= 0 ? 1 : 0;
    }
    double areaarc(double th, double r){
        return 0.5 * r * r * (th - sin(th));
    }
    // 圆的面积并，去掉内含的圆
    void init_or(){
        int i, j, k = 0;
        bool mark[N] = {0};
        for(i = 0; i < n; i++){
            for(j = 0; j < n;j ++){
                if(i != j && !mark[j]){
                    if((c[i] == c[j]) || inner(c[i], c[j])) break;
                }
            }
            if(j < n) mark[i] = 1;
        }
        for(i = 0; i < n; i++){
            if(!mark[i]){
                c[k++] = c[i];
            }
        }
        n = k;
    }
    // 圆的面积交，去掉内含的圆
    void init_and(){
        int i, j, k = 0;
        bool mark[N] = {0};
        for(i = 0; i < n; i++){
            for(j = 0; j < n;j ++){
                if(i != j && !mark[j]){
                    if((c[i] == c[j]) || inner(c[j], c[i])) break;
                }
            }
            if(j < n) mark[i] = 1;
        }
        for(i = 0; i < n; i++){
            if(!mark[i]){
                c[k++] = c[i];
            }
        }
        n = k;
    }
    void getarea(){
        memset(ans, 0, sizeof ans);
        vector<pair<double,int> > v;
        for(int i=0;i<n;i++){
            v.clear();
            v.push_back(make_pair(-pi, 1));
            v.push_back(make_pair(pi, -1));   
            for(int j=0;j<n;j++){
                if(i != j){
                    Point q = (c[j].p - c[i].p);
                    db ab = q.len(), ac = c[i].r, bc = c[j].r;
                    if(sgn(ab + ac - bc) <= 0){
                        v.push_back(make_pair(-pi, 1));
                        v.push_back(make_pair(pi, -1));
                        continue;
                    }
                    if(sgn(ab + bc - ac) <= 0) continue;
                    if(sgn(ab - ac - bc) > 0) continue;
                    double th = atan2(q.y, q.x), fai = acos((ac*ac + ab * ab - bc * bc) / (2.0*ac*ab));
                    double a0 = th - fai;
                    if(sgn(a0 + pi) < 0) a0 += 2 * pi;
                    db a1 = th + fai;
                    if(sgn(a1 - pi) > 0) a1 -= 2 * pi;
                    if(sgn(a0 - a1) > 0){
                        v.push_back(make_pair(a0, 1));
                        v.push_back(make_pair(pi, -1));
                        v.push_back(make_pair(-pi, 1));
                        v.push_back(make_pair(a1, -1));
                    }
                    else {
                        v.push_back(make_pair(a0, 1));
                        v.push_back(make_pair(a1, -1));
                    }
                }
            }
            sort(v.begin(),v.end());
            int cur = 0;
            for(int j=0;j<v.size();j++){
                if(cur && sgn(v[j].first - pre[cur])){
                    ans[cur] += areaarc(v[j].first - pre[cur], c[i].r);
                    ans[cur] += 0.5 * (c[i].point(pre[cur])^c[i].point(v[j].first));
                }
                cur += v[j].second;
                pre[cur] = v[j].first;
            }
        }
    }
}cs;

4. 多边形

struct polygon{
    int n;
    Point p[maxp];
    Line l[maxp];
    void input(int n){
        this->n=n;
        for(int i=0;i<n;i++) p[i].input();
    }
}

基本函数

// 得到所有边
void getline(){
    for(int i=0;i<n;i++){
        l[i] = Line(p[i],p[(i+1)%n]);
    }
}
// 以 p0 为标准极角排序
struct cmp{
    Point p;
    cmp(const Point &p0){p = p0;}
    bool operator()(const Point &aa, const Point &bb){
        Point a = aa, b = bb;
        int d = sgn((a-p)^(b-p));
        if(d == 0){
            return sgn(a.distance(p) - b.distance(p)) < 0;
        }
        return d > 0;
    }
};
/*
        进行极角排序
        首先找打最左下角的点
        需要重载好Point的 < 操作符
*/
void norm(){
    Point mi = p[0];
    for(int i=1;i<n;i++) mi = min(mi, p[i]);
    sort(p, p+n, cmp(mi));
}
// 得到周长
db getcircumference(){
    db sum = 0;
    for(int i=0;i<n;i++){
        sum += p[i].distance(p[(i+1)%n]);
    }
    return sum;
}
// 得到面积
db getarea(){
    db sum = 0;
    // 以原点为划分点
    for(int i=0;i<n;i++){
        sum += (p[i]^p[(i+1)%n]);
    }
    return fabs(sum)/2;
}
// 得到方向，1 表示逆时针，0表示顺时针
bool getdir(){
    db sum = 0;
    for(int i=0;i<n;i++){
        sum += (p[i] ^ p[(i+1)%n]);
    }
    if(sgn(sum) > 0) return 1;
    return 0;
}
// 得到重心
Point getbarycentre(){
    Point ret(0,0);
    db area = 0;
    for(int i=1;i<n-1;i++){
        db tmp = (p[i]-p[0])^(p[i+1]-p[0]);
        if(sgn(tmp) == 0) continue;
        area += tmp;
        ret.x += (p[0].x + p[i].x + p[i+1].x) / 3 * tmp;
        ret.y += (p[0].y + p[i].y + p[i+1].y) / 3 * tmp;
    }
    if(sgn(area)) ret = ret / area;
    return ret;
}

4.1 凸包

/*
     得到凸包 凸包点编号0 ~ n-1
*/
void getconvex(polygon &convex){
    sort(p, p+n);
    convex.n = n;
    for(int i=0;i<min(n, 2);i++){
        convex.p[i] = p[i];
    }
    if(convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n --;
    if(n <= 2) return;
    int &top = convex.n;
    top = 1;
    for(int i=2;i<n;i++){
        while(top && sgn((convex.p[top] - p[i]) ^ (convex.p[top-1] - p[i])) <= 0)
            top--;
        convex.p[++top] = p[i];
    }
    int temp = top;
    convex.p[++top] = p[n-2];
    for(int i=n-3;i>=0;i--){
        while(top != temp && sgn((convex.p[top] - p[i]) ^ (convex.p[top-1] - p[i])) <= 0)
            top--;
        convex.p[++top] = p[i];
    }
    if(convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n --;
    convex.norm();
}
void Graham(polygon &convex){
    norm();
    int &top = convex.n;
    top = 0;
    if(n == 1){
        top = 1;
        convex.p[0] = p[0];
        return ;
    }
    if(n == 2){
        top = 2;
        convex.p[0] = p[0];
        convex.p[1] = p[1];
        if(convex.p[0] == convex.p[1]) top--;
        return;
    }
    convex.p[0] = p[0];
    convex.p[1] = p[1];
    top = 2;
    for(int i=2;i<n;i++){
        while(top > 1 && sgn((convex.p[top-1] - convex.p[top-2]) ^ (p[i]-convex.p[top-2])) <= 0) 
            top--;
        convex.p[top++] = p[i];
    }
    if(convex.n == 2 && (convex.p[0] == convex.p[1])) convex.n--;
}

4.2 判断是不是凸

// 判断是不是凸的
bool isconvex(){
    bool s[3];
    memset(s, false, sizeof s);
    for(int i=0;i<n;i++){
        int j = (i + 1) % n;
        int k = (j + 1) % n;
        s[sgn((p[j] - p[i]) ^ (p[k] - p[i])) + 1] = true;
        if(s[0] && s[2]) return false;
    }
    return true;
}

4.3 判断点和任意多边形的关系

/*
        判断点和任意多边形的关系
        3 点上
        2 边上
        1 内部
        0 外部
*/
int relationpoint(Point q){
    for(int i=0;i<n;i++){
        if(p[i] == q) return 3; 
    }
    getline();
    for(int i=0;i<n;i++){
        if(l[i].pointonseg(q)) return 2;
    }
    int cnt = 0;
    for(int i=0;i<n;i++){
        int j = (i + 1) % n;
        int k = sgn((q-p[j])^(p[i]-p[j]));
        int u = sgn(p[i].y - q.y);
        int v = sgn(p[j].y - q.y);
        if(k > 0 && u < 0 && v >= 0) cnt ++;
        if(k < 0 && v < 0 && u >= 0) cnt --;
    }
    return cnt != 0;
}

4.4 直线 u 切割 凸多边形左侧

//直线 u 切割凸多边形左侧
//注意直线方向
//HDU3982

void convexcut(Line u, polygon &po){
    int &top = po.n;
    top = 0;
    for(int i=0;i<n;i++){
        int d1 = sgn((u.e-u.s)^(p[i]-u.s));
        int d2 = sgn((u.e-u.s)^(p[(i+1)%n]-u.s));
        if(d1 >= 0) po.p[top++] = p[i];
        if(d1 * d2 < 0) po.p[top++] = u.crosspoint(Line(p[i], p[(i+1)%n]));
    }
}

4.5 多边形与圆交的面积

// 多边形和圆交的面积
db areacircle(circle c){
    db ans = 0;
    for(int i=0;i<n;i++){
        int j = (i + 1) % n;
        if(sgn((p[j]-c.p)^(p[i]-c.p)) >= 0)
            ans += c.areatriangle(p[i], p[j]);
        else ans -= c.areatriangle(p[i], p[j]);
    }
    return fabs(ans);
}

4.6 多边形与圆的关系

/*
        多边形与圆的关系
        2. 圆完全在多边形内
        1. 圆在多边形里面，碰到了多边形边界
        0. 其他
*/
db relationcircle(circle c){
    getline();
    int x = 2;
    if(relationpoint(c.p) != 1) return 0; // 圆心不在内部
    for(int i=0;i<n;i++){
        if(c.relationseg(l[i]) == 2) return 0;
        if(c.relationseg(l[i]) == 1) x = 1; // 相切
    }
    return x;
}

4.7 最小矩形覆盖

db minRectangleCover(){
    if(n < 3) return 0.0;
    p[n] = p[0];
    db ans = -1;
    int up = 1, r = 1, l;
    for(int i=0;i<n;i++){
        Point vec = p[i+1] - p[i];
        while(sgn( (vec^(p[up+1]-p[i])) - (vec^(p[up]-p[i])) ) >= 0)
            up = (up + 1) % n;
        while(sgn( (vec*(p[r+1]-p[i])) - (vec * (p[r]-p[i])) ) >= 0)
            r = (r + 1) % n;
        if(i == 0) l = r;
        while(sgn( (vec*(p[l+1]-p[i])) - (vec * (p[l]-p[i])) ) <= 0)
            l = (l + 1) % n;

        db d = (p[i] - p[i+1]).len2();
        db tmp = (vec^(p[up]-p[i])) * ( (vec*(p[r]-p[i])) - (vec *(p[l]-p[i]))) / d;
        if(ans < 0 || ans > tmp) ans = tmp;
    }
    return ans;
}

4.8 平面最远点对

// 最远的一对点的距离
db diameter(){
    if(n == 2) return p[0].distance(p[1]);
    int i = 0, j = 0;
    for(int k=0;k<n;k++){
        if(p[i] < p[k]) i = k;
        if(!(p[j] < p[k])) j = k;
    }
    int si = i, sj = j;
    db res = 0;
    while(i != sj || j != si){
        res = max(res, p[i].distance(p[j]));
        int ni = (i+1)%n;
        int nj = (j+1)%n;
        if(sgn((p[ni]-p[i])^(p[nj]-p[j])) <= 0){
            i = ni;
        } else j = nj;
    }
    return res;
}