geo - poligon cut

typedef vector<point> Polygon;
struct Line {
    point A,B;
};
// Cut a polygon with a line. Returns one half.
// To return the other half, reverse the direction of Line l (by negating l.a, l.b)
// The line must be formed using 2 points
int ccw(point a,point b,point c)
{
    return dcmp(cross(b-a,c-a));
}
Polygon polygon_cut(const Polygon& P, Line l) {
    Polygon Q;
    for(int i = 0; i < P.size(); ++i) {
        point A = P[i], B = (i == P.size()-1) ? P[0] : P[i+1];
        if (ccw(l.A, l.B, A) != -1) Q.push_back(A);
        if (ccw(l.A, l.B, A)*ccw(l.A, l.B, B) < 0) {
            point p = LLIN(A,B,l.A, l.B); //areIntersect(Line(A, B), l, p);
            Q.push_back(p);
        }
    }
    return Q;
}
point Reflect(Line l, point p)
{
    point q=  PPL(l.A, l.B, p);
    return p + (q - p) * 2;
}

point MakeCircle(point a,point b, point c)
{
    if(dcmp(cross(b-a,c-a)) == 0) {
        point p[3];
        p[0] = a, p[1] = b, p[2] = c;
        sort(p ,p + 3);
        return (p[0] + p[2])/2;
    }
    return CCC(a,b,c);
}
void SmallEnclosingCircle( Polygon p, point &cen, double &r)
{
    random_shuffle(p.begin(), p.end());
    cen = p[0], r = 0;
    int n = p.size();
    for(int i= 1;i < n;i ++){
        if(dist(cen,p[i]) <= r + EPS) continue;
        cen = p[i], r = 0;
        for(int j = 0; j < i; j ++) {
            if(dist(cen, p[j]) <= r + EPS) continue;
            cen = (p[i] + p[j])/2;
            r = dist(cen, p[j]);
            for(int k = 0; k < j ;k ++) {
                if(dist(cen, p[k]) <= r + EPS) continue;
                cen = MakeCircle(p[i],p[j],p[k]);
                r = dist(cen, p[k]);
            }
        }
    }
}

struct Circle{
    point p;
    double r;
};
// helper functions for commonCircleArea
double cir_area_solve(double a, double b, double c) {
    return acos((a*a + b*b - c*c) / 2 / a / b);
}
double cir_area_cut(double a, double r) {
    double s1 = a * r * r / 2;
    double s2 = sin(a) * r * r / 2;
    return s1 - s2;
}
// Tested: http://codeforces.com/contest/600/problem/D
double commonCircleArea(Circle c1, Circle c2) { //return the common area of two circle
    if (c1.r < c2.r) swap(c1, c2);
    double d = mag(c1.p - c2.p);
    if (d + c2.r <= c1.r + EPS) return c2.r*c2.r*PI;
    if (d >= c1.r + c2.r - EPS) return 0.0;
    double a1 = cir_area_solve(d, c1.r, c2.r);
    double a2 = cir_area_solve(d, c2.r, c1.r);
    return cir_area_cut(a1*2, c1.r) + cir_area_cut(a2*2, c2.r);
}