算法模板の计算几何

1、凸包

inline bool cmp(const POINT &a, const POINT &b) {
    if(a.y == b.y) return a.x < b.x;
    return a.y < b.y;
}
//turn left
inline bool Cross(POINT &sp, POINT &ep, POINT &op) {
    return (sp.x - op.x) * (ep.y - op.y) - (ep.x - op.x) * (sp.y - op.y) >= 0;
}
  
inline double dist(POINT &a, POINT &b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
  
void Graham_scan() {
    std::sort(p, p + n, cmp);
    top = 1;
    stk[0] = 0; stk[1] = 1;
    for(int i = 2; i < n; ++i) {
        while(top && Cross(p[i], p[stk[top]], p[stk[top - 1]])) --top;
        stk[++top] = i;
    }
    int len = top;
    stk[++top] = n - 2;
    for(int i = n - 3; i >= 0; --i) {
        while(top != len && Cross(p[i], p[stk[top]], p[stk[top - 1]])) --top;
        stk[++top] = i;
    }
}

2、旋转卡壳（POJ 3384）

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
 
#define EPS 1e-8
#define MAXN 1000
 
inline int sgn(double x) {
    if(fabs(x) < EPS) return 0;
    return x > 0 ? 1 : -1;
}
 
struct Point {
    double x, y;
    Point(double xx = 0, double yy = 0): x(xx), y(yy) {}
    bool operator == (const Point &b) const {
        return sgn(x - b.x) == 0 && sgn(y - b.y) == 0;
    }
};
//cross
inline double operator ^ (const Point &a, const Point &b) {
    return a.x * b.y - a.y * b.x;
}
 
inline Point operator - (const Point &a, const Point &b) {
    return Point(a.x - b.x, a.y - b.y);
}
 
struct Line {
    Point s, e;
    double ag;
};
 
struct polygon {
    Point v[MAXN];
    int n;
} pg, res;
 
inline double dist(Point &a, Point &b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
 
inline double Cross(Point o, Point s, Point e) {
    return (s - o) ^ (e - o);
}
//cross_point
Point operator * (const Line &a, const Line &b) {
    Point res;
    double u = Cross(a.s, a.e, b.s), v = Cross(a.e, a.s, b.e);
    res.x = (b.s.x * v + b.e.x * u)/(u + v);
    res.y = (b.s.y * v + b.e.y * u)/(u + v);
    return res;
}
 
int parallel(Line a, Line b) {
    double u = (a.e.x - a.s.x) * (b.e.y - b.s.y) - (a.e.y - a.s.y) * (b.e.x - b.s.x);
    return sgn(u) == 0;
}
 
inline void set_vector(double x1, double y1, double x2, double y2, Line &v) {
    v.s.x = x1; v.s.y = y1;
    v.e.x = x2; v.e.y = y2;
    v.ag = atan2(y2 - y1, x2 - x1);
}
 
Line vct[MAXN], deq[MAXN];
 
bool cmp(const Line &a, const Line &b) {
    if(sgn(a.ag - b.ag) == 0)
        return sgn(Cross(b.s, b.e, a.s)) < 0;
    return a.ag < b.ag;
}
 
int half_planes_cross(Line *v, int vn) {
    int i, n;
    //sort(v, v + vn, cmp);
    for(i = n = 1; i < vn; ++i) {
        if(sgn(v[i].ag - v[i-1].ag) == 0) continue;
        v[n++] = v[i];
    }
    int head = 0, tail = 1;
    deq[0] = v[0], deq[1] = v[1];
    for(i = 2; i < n; ++i) {
        if(parallel(deq[tail - 1], deq[tail]) || parallel(deq[head], deq[head + 1]))
            return false;
        while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[tail - 1] * deq[tail])) > 0)
            --tail;
        while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[head] * deq[head + 1])) > 0)
            ++head;
        deq[++tail] = v[i];
    }
    while(head < tail && sgn(Cross(deq[head].s, deq[head].e, deq[tail - 1] * deq[tail])) > 0)
        --tail;
    while(head < tail && sgn(Cross(deq[tail].s, deq[tail].e, deq[head] * deq[head + 1])) > 0)
        ++head;
    if(tail <= head + 1) return false;
    res.n = 0;
    for(i = head; i < tail; ++i)
        res.v[res.n++] = deq[i] * deq[i + 1];
    res.v[res.n++] = deq[head] * deq[tail];
    res.n = unique(res.v, res.v + res.n) - res.v;
    res.v[res.n] = res.v[0];
    return true;
}
 
void moving(Line v[], int vn, double r) {
    for(int i = 0; i < vn; ++i) {
        double dx = v[i].e.x - v[i].s.x, dy = v[i].e.y - v[i].s.y;
        dx = dx / dist(v[i].s, v[i].e) * r;
        dy = dy / dist(v[i].s, v[i].e) * r;
        v[i].s.x += dy; v[i].e.x += dy;
        v[i].s.y -= dx; v[i].e.y -= dx;
    }
}
 
int ix, jx;
 
double dia_roataing_calipers() {
    double dia = 0;
    ix = jx = 0;
    int q = 1;
    for(int i = 0; i < res.n; ++i) {
        while(sgn(Cross(res.v[i+1], res.v[i], res.v[q+1]) - Cross(res.v[i+1], res.v[i], res.v[q])) > 0)
            q = (q + 1) % res.n;
        if(sgn(dist(res.v[i], res.v[q]) - dia) > 0) {
            dia = dist(res.v[i], res.v[q]);
            ix = i; jx = q;
        }
        if(sgn(dist(res.v[i+1], res.v[q]) - dia) > 0) {
            dia = dist(res.v[i+1], res.v[q]);
            ix = i+1; jx = q;
        }
    }
    return dia;
}
 
int main() {
    int n;
    double r;
    while(scanf("%d%lf", &n, &r) != EOF) {
        for(int i = 0; i < n; ++i) scanf("%lf%lf", &pg.v[i].x, &pg.v[i].y);
        pg.v[n] = pg.v[0];
        for(int i = 0; i < n; ++i)
            set_vector(pg.v[i].x, pg.v[i].y, pg.v[i+1].x, pg.v[i+1].y, vct[i]);
        moving(vct, n, r);
        half_planes_cross(vct, n);
        dia_roataing_calipers();
        printf("%.4f %.4f %.4f %.4f
", res.v[ix].x, res.v[ix].y, res.v[jx].x, res.v[jx].y);
    }
}

3、最小圆覆盖

#include <cstdio>
#include <cmath>

const int MAXN = 1010;

struct Point {
    double x, y;
    Point(double xx = 0, double yy = 0): x(xx), y(yy) {}
};

double operator ^ (const Point &a, const Point &b) {
    return a.x * b.y - a.y * b.x;
}

Point operator - (const Point &a, const Point &b) {
    return Point(a.x - b.x, a.y - b.y);
}

double dist(const Point &a, const Point &b) {
    double x = a.x - b.x, y = a.y - b.y;
    return sqrt(x * x + y * y);
}

struct Circle {
    double r;
    Point centre;
};

struct Triangle {
    Point t[3];
    double Area() {
        return fabs((t[1] - t[0]) ^ (t[2] - t[0]))/2;
    }
};

Circle circumcircleOfTriangle(Triangle t) {
    Circle tmp;
    double a, b, c, c1, c2;
    Point A(t.t[0].x, t.t[0].y);
    Point B(t.t[1].x, t.t[1].y);
    Point C(t.t[2].x, t.t[2].y);
    a = dist(B, C);
    b = dist(A, C);
    c = dist(A, B);
    tmp.r = a * b * c / t.Area() / 4;
    c1 = (A.x * A.x + A.y * A.y - B.x * B.x - B.y * B.y) / 2;
    c2 = (A.x * A.x + A.y * A.y - C.x * C.x - C.y * C.y) / 2;
    tmp.centre.x = (c1 * (A.y - C.y) - c2 * (A.y - B.y)) / ((A.x - B.x) * (A.y - C.y) - (A.x - C.x) * (A.y - B.y));
    tmp.centre.y = (c1 * (A.x - C.x) - c2 * (A.x - B.x)) / ((A.y - B.y) * (A.x - C.x) - (A.y - C.y) * (A.x - B.x));
    return tmp;
}

Circle c;
Point a[MAXN];

Circle MinCircle2(int tce, Triangle ce) {
    Circle tmp;
    if(tce == 0) tmp.r = -2;
    else if(tce == 1) {
        tmp.centre = ce.t[0];
        tmp.r = 0;
    }
    else if(tce == 2) {
        tmp.r = dist(ce.t[0], ce.t[1]) / 2;
        tmp.centre.x = (ce.t[0].x + ce.t[1].x) / 2;
        tmp.centre.y = (ce.t[0].y + ce.t[1].y) / 2;
    }
    else if(tce == 3) tmp = circumcircleOfTriangle(ce);
    return tmp;
}

void MinCircle(int t, int tce, Triangle ce) {
    Point tmp;
    c = MinCircle2(tce, ce);
    if(tce == 3) return;
    for(int i = 1; i <= t; ++i) {
        if(dist(a[i], c.centre) > c.r) {
            ce.t[tce] = a[i];
            MinCircle(i - 1, tce + 1, ce);
            tmp = a[i];
            for(int j = i; j >= 2; --j) a[j] = a[j - 1];
            a[1] = tmp;
        }
    }
}

void run(int n) {
    Triangle ce;
    MinCircle(n, 0, ce);
    printf("%.2f %.2f %.2f
", c.centre.x, c.centre.y, c.r);
}

int main() {
    int n;
    while(scanf("%d", &n) != EOF) {
        if(n == 0) break;
        for(int i = 1; i <= n; ++i)
            scanf("%lf%lf", &a[i].x, &a[i].y);
        run(n);
    }
}

4、半平面交+最远点对

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;

#define EPS 1e-8
#define MAXN 1000

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

struct Point {
    double x, y;
    Point(double xx = 0, double yy = 0): x(xx), y(yy) {}
    bool operator == (const Point &b) const {
        return sgn(x - b.x) == 0 && sgn(y - b.y) == 0;
    }
};
//cross
inline double operator ^ (const Point &a, const Point &b) {
    return a.x * b.y - a.y * b.x;
}

inline Point operator - (const Point &a, const Point &b) {
    return Point(a.x - b.x, a.y - b.y);
}

struct Line {
    Point s, e;
    double ag;
};

struct polygon {
    Point v[MAXN];
    int n;
} pg, res;

inline double dist(Point &a, Point &b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

inline double Cross(Point o, Point s, Point e) {
    return (s - o) ^ (e - o);
}
//cross_point
Point operator * (const Line &a, const Line &b) {
    Point res;
    double u = Cross(a.s, a.e, b.s), v = Cross(a.e, a.s, b.e);
    res.x = (b.s.x * v + b.e.x * u)/(u + v);
    res.y = (b.s.y * v + b.e.y * u)/(u + v);
    return res;
}

int parallel(Line a, Line b) {
    double u = (a.e.x - a.s.x) * (b.e.y - b.s.y) - (a.e.y - a.s.y) * (b.e.x - b.s.x);
    return sgn(u) == 0;
}

inline void set_vector(double x1, double y1, double x2, double y2, Line &v) {
    v.s.x = x1; v.s.y = y1;
    v.e.x = x2; v.e.y = y2;
    v.ag = atan2(y2 - y1, x2 - x1);
}

Line vct[MAXN], deq[MAXN];

bool cmp(const Line &a, const Line &b) {
    if(sgn(a.ag - b.ag) == 0)
        return sgn(Cross(b.s, b.e, a.s)) < 0;
    return a.ag < b.ag;
}

int half_planes_cross(Line *v, int vn) {
    int i, n;
    //sort(v, v + vn, cmp);
    for(i = n = 1; i < vn; ++i) {
        if(sgn(v[i].ag - v[i-1].ag) == 0) continue;
        v[n++] = v[i];
    }
    int head = 0, tail = 1;
    deq[0] = v[0], deq[1] = v[1];
    for(i = 2; i < n; ++i) {
        if(parallel(deq[tail - 1], deq[tail]) || parallel(deq[head], deq[head + 1]))
            return false;
        while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[tail - 1] * deq[tail])) > 0)
            --tail;
        while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[head] * deq[head + 1])) > 0)
            ++head;
        deq[++tail] = v[i];
    }
    while(head < tail && sgn(Cross(deq[head].s, deq[head].e, deq[tail - 1] * deq[tail])) > 0)
        --tail;
    while(head < tail && sgn(Cross(deq[tail].s, deq[tail].e, deq[head] * deq[head + 1])) > 0)
        ++head;
    if(tail <= head + 1) return false;
    res.n = 0;
    for(i = head; i < tail; ++i)
        res.v[res.n++] = deq[i] * deq[i + 1];
    res.v[res.n++] = deq[head] * deq[tail];
    res.n = unique(res.v, res.v + res.n) - res.v;
    res.v[res.n] = res.v[0];
    return true;
}

void moving(Line v[], int vn, double r) {
    for(int i = 0; i < vn; ++i) {
        double dx = v[i].e.x - v[i].s.x, dy = v[i].e.y - v[i].s.y;
        dx = dx / dist(v[i].s, v[i].e) * r;
        dy = dy / dist(v[i].s, v[i].e) * r;
        v[i].s.x += dy; v[i].e.x += dy;
        v[i].s.y -= dx; v[i].e.y -= dx;
    }
}

int ix, jx;

double dia_roataing_calipers() {
    double dia = 0;
    ix = jx = 0;
    int q = 1;
    for(int i = 0; i < res.n; ++i) {
        while(sgn(Cross(res.v[i+1], res.v[i], res.v[q+1]) - Cross(res.v[i+1], res.v[i], res.v[q])) > 0)
            q = (q + 1) % res.n;
        if(sgn(dist(res.v[i], res.v[q]) - dia) > 0) {
            dia = dist(res.v[i], res.v[q]);
            ix = i; jx = q;
        }
        if(sgn(dist(res.v[i+1], res.v[q]) - dia) > 0) {
            dia = dist(res.v[i+1], res.v[q]);
            ix = i+1; jx = q;
        }
    }
    return dia;
}

int main() {
    int n;
    double r;
    while(scanf("%d%lf", &n, &r) != EOF) {
        for(int i = 0; i < n; ++i) scanf("%lf%lf", &pg.v[i].x, &pg.v[i].y);
        pg.v[n] = pg.v[0];
        for(int i = 0; i < n; ++i)
            set_vector(pg.v[i].x, pg.v[i].y, pg.v[i+1].x, pg.v[i+1].y, vct[i]);
        moving(vct, n, r);
        half_planes_cross(vct, n);
        dia_roataing_calipers();
        printf("%.4f %.4f %.4f %.4f
", res.v[ix].x, res.v[ix].y, res.v[jx].x, res.v[jx].y);
    }
}

其他1：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef double TYPE;
#define Abs(x) (((x)>0)?(x):(-(x)))
#define Sgn(x) (((x)<0)?(-1):(1))
#define Max(a,b) (((a)>(b))?(a):(b))
#define Min(a,b) (((a)<(b))?(a):(b))
#define EPS 1e-8
#define Infinity 1e+10
#define PI acos(-1.0)//3.14159265358979323846
TYPE Deg2Rad(TYPE deg){return (deg * PI / 180.0);}
TYPE Rad2Deg(TYPE rad){return (rad * 180.0 / PI);}
TYPE Sin(TYPE deg){return sin(Deg2Rad(deg));}
TYPE Cos(TYPE deg){return cos(Deg2Rad(deg));}
TYPE ArcSin(TYPE val){return Rad2Deg(asin(val));}
TYPE ArcCos(TYPE val){return Rad2Deg(acos(val));}
TYPE Sqrt(TYPE val){return sqrt(val);}

//点
struct POINT {
    TYPE x, y;
    POINT(TYPE xx = 0, TYPE yy = 0): x(xx), y(yy) {}
};
// 两个点的距离
TYPE Distance(const POINT &a, const POINT &b) {
    return Sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
//线段
struct SEG {
    POINT a ,b;
    SEG() {};
    SEG(POINT aa, POINT bb): a(aa), b(bb) {}
};
//直线(两点式)
struct LINE {
    POINT a, b;
    LINE() {};
    LINE(POINT aa, POINT bb): a(aa), b(bb) {}
};
//直线(一般式)
struct LINE2 {
    TYPE A,B,C;
    LINE2() {};
    LINE2(TYPE AA, TYPE BB, TYPE CC): A(AA), B(BB), C(CC) {}
};
//两点式化一般式
LINE2 Line2line(const LINE &L) {
    LINE2 L2;
    L2.A = L.b.y - L.a.y;
    L2.B = L.a.x - L.b.x;
    L2.C = L.b.x * L.a.y - L.a.x * L.b.y;
    return L2;
}
// 引用返回直线 Ax + By + C =0 的系数
void Coefficient(const LINE &L, TYPE &A, TYPE &B, TYPE &C) {
    A = L.b.y - L.a.y;
    B = L.a.x - L.b.x;
    C = L.b.x * L.a.y - L.a.x * L.b.y;
}
void Coefficient(const POINT &p,const TYPE &a,TYPE &A,TYPE &B,TYPE &C) {
    A = Cos(a);
    B = Sin(a);
    C = - (p.y * B + p.x * A);
}
//直线相交的交点
POINT Intersection(const LINE &A, const LINE &B) {
    TYPE A1, B1, C1;
    TYPE A2, B2, C2;
    Coefficient(A, A1, B1, C1);
    Coefficient(B, A2, B2, C2);
    POINT I;
    I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);
    I.y =   (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1);
    return I;
}
//点到直线的距离
TYPE Ptol(const POINT &p,const LINE2 &L) {
    return fabs(L.A * p.x + L.B * p.y + L.C)/Sqrt(L.A * L.A + L.B * L.B);
}
//两线段叉乘
TYPE Cross2(const SEG &p, const SEG &q) {
    return (p.b.x - p.a.x) * (q.b.y - q.a.y) - (q.b.x - q.a.x) * (p.b.y - p.a.y);
}
//有公共端点O叉乘,并判拐,若正p0->p1->p2拐向左
TYPE Cross(const POINT &a, const POINT &b, const POINT &o) {
    return (a.x - o.x) * (b.y - o.y) - (b.x - o.x) * (a.y - o.y);
}
//判等(值，点，直线)
bool IsEqual(const TYPE &a, const TYPE &b) {
    return (Abs(a - b) < EPS);
}
bool IsEqual(const POINT & a, const POINT & b) {
    return IsEqual(a.x, b.x) && IsEqual(a.y, b.y);
}
bool IsEqual(const LINE & A, const LINE & B) {
    TYPE A1, B1, C1;
    TYPE A2, B2, C2;
    Coefficient(A, A1, B1, C1);
    Coefficient(B, A2, B2, C2);
    return IsEqual(A1 * B2, A2 * B1) && IsEqual(A1 * C2, A2 * C1) && IsEqual(B1 * C2, B2 * C1);
}
// 判断点是否在线段上
bool IsOnSeg(const SEG &seg, const POINT &p) {
    return (IsEqual(p, seg.a) || IsEqual(p, seg.b)) ||
        (((p.x - seg.a.x) * (p.x - seg.b.x) < 0 ||
        (p.y - seg.a.y) * (p.y - seg.b.y) < 0) &&
        (IsEqual(Cross(seg.b, p, seg.a), 0)));
}

//判断两条线断是否相交，端点重合算相交
bool IsIntersect(const SEG &u, const SEG &v) {
    return (Cross(v.a, u.b, u.a) * Cross(u.b, v.b, u.a) >= 0) &&
        (Cross(u.a, v.b, v.a) * Cross(v.b, u.b, v.a) >= 0) &&
        (Max(u.a.x, u.b.x) >= Min(v.a.x, v.b.x)) &&
        (Max(v.a.x, v.b.x) >= Min(u.a.x, u.b.x)) &&
        (Max(u.a.y, u.b.y) >= Min(v.a.y, v.b.y)) &&
        (Max(v.a.y, v.b.y) >= Min(u.a.y, u.b.y));
}
//判断线段P和直线Q是否相交，端点重合算相交
bool Isonline(const SEG &p, const LINE &q) {
    SEG p1,p2,q0;
    p1.a = q.a; p1.b = p.a;
    p2.a = q.a; p2.b = p.b;
    q0.a = q.a; q0.b = q.b;
    return (Cross2(p1,q0) * Cross2(q0,p2) >= 0);
}
//判断两条线断是否平行
bool IsParallel(const LINE &A, const LINE &B) {
    TYPE A1, B1, C1;
    TYPE A2, B2, C2;
    Coefficient(A, A1, B1, C1);
    Coefficient(B, A2, B2, C2);
    return (A1 * B2 == A2 * B1);
}
//判断两条直线是否相交
bool IsIntersect(const LINE & A, const LINE & B) {
  return !IsParallel(A, B);
}
// 矩形
struct RECT {
    POINT a; // 左下点
    POINT b; // 右上点
    RECT() {}
    RECT(const POINT &aa, const POINT &bb): a(aa), b(bb) {}
};
//矩形化标准
RECT Stdrect(const RECT &q) {
    RECT p = q;
    if(p.a.x > p.b.x) swap(p.a.x , p.b.x);
    if(p.a.y > p.b.y) swap(p.a.y , p.b.y);
    return p;
}
//根据下标返回矩形的边
SEG Edge(const RECT &rect, int idx) {
    SEG edge;
    while (idx < 0) idx += 4;
    switch (idx % 4) {
        case 0: //下边
            edge.a = rect.a;
            edge.b = POINT(rect.b.x, rect.a.y);
            break;
        case 1: //右边
            edge.a = POINT(rect.b.x, rect.a.y);
            edge.b = rect.b;
            break;
        case 2: //上边
            edge.a = rect.b;
            edge.b = POINT(rect.a.x, rect.b.y);
            break;
        case 3: //左边
            edge.a = POINT(rect.a.x, rect.b.y);
            edge.b = rect.a;
            break;
        default:
            break;
    }
    return edge;
}
//矩形的面积
TYPE Area(const RECT &rect) {
    return (rect.b.x - rect.a.x) * (rect.b.y - rect.a.y);
}
//两个矩形的公共面积
TYPE CommonArea(const RECT &A, const RECT &B) {
    TYPE area = 0.0;
    POINT LL(Max(A.a.x, B.a.x), Max(A.a.y, B.a.y));
    POINT UR(Min(A.b.x, B.b.x), Min(A.b.y, B.b.y));
    if((LL.x <= UR.x) && (LL.y <= UR.y)) {
        area = Area(RECT(LL, UR));
    }
    return area;
}
// 圆
struct CIRCLE {
    TYPE x, y, r;
    CIRCLE() {}
    CIRCLE(TYPE xx, TYPE yy, TYPE zz): x(xx), y(yy), r(zz) {}
};
//圆心
POINT Center(const CIRCLE &circle) {
    return POINT(circle.x, circle.y);
}
//圆的面积
TYPE Area(const CIRCLE &circle) {
    return PI * circle.r * circle.r;
}
//两个圆的公共面积
TYPE CommonArea(const CIRCLE &A, const CIRCLE &B)
{
    TYPE area = 0.0;
    const CIRCLE & M = (A.r > B.r) ? A : B;
    const CIRCLE & N = (A.r > B.r) ? B : A;
    TYPE D = Distance(Center(M), Center(N));
    if((D < M.r + N.r) && (D > M.r - N.r)) {
        TYPE cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D);
        TYPE cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D);
        TYPE alpha = 2.0 * ArcCos(cosM);
        TYPE beta = 2.0 * ArcCos(cosN);
        TYPE TM = 0.5 * M.r * M.r * Sin(alpha);
        TYPE TN = 0.5 * N.r * N.r * Sin(beta);
        TYPE FM = (alpha / 360.0) * Area(M);
        TYPE FN = (beta / 360.0) * Area(N);
        area = FM + FN - TM - TN;
    }
    else if(D <= M.r - N.r) {
        area = Area(N);
    }
    return area;
}
//判断圆是否在矩形内(不允许相切)
bool IsInCircle(const CIRCLE &circle, const RECT &rect) {
    return (circle.x - circle.r > rect.a.x) &&
        (circle.x + circle.r < rect.b.x) &&
        (circle.y - circle.r > rect.a.y) &&
        (circle.y + circle.r < rect.b.y);
}
//判断矩形是否在圆内(不允许相切)
bool IsInRect(const CIRCLE & circle, const RECT & rect) {
    POINT c,d;
    c.x = rect.a.x; c.y = rect.b.y;
    d.x = rect.b.x; d.y = rect.a.y;
    return (Distance( Center(circle) , rect.a ) < circle.r) &&
        (Distance( Center(circle) , rect.b ) < circle.r) &&
        (Distance( Center(circle) , c ) < circle.r) &&
        (Distance( Center(circle) , d ) < circle.r);
}
//判断矩形是否与圆相离(不允许相切)
bool Isoutside(const CIRCLE & circle, const RECT & rect) {
    POINT c,d;
    c.x = rect.a.x; c.y=rect.b.y;
    d.x = rect.b.x; d.y=rect.a.y;
    return ((Distance( Center(circle) , rect.a ) > circle.r) &&
        (Distance( Center(circle) , rect.b ) > circle.r) &&
        (Distance( Center(circle) , c ) > circle.r) &&
        (Distance( Center(circle) , d ) > circle.r) &&
        (rect.a.x > circle.x || circle.x > rect.b.x || rect.a.y > circle.y || circle.y > rect.b.y)) ||
        ((circle.x - circle.r > rect.b.x) ||
        (circle.x + circle.r < rect.a.x) ||
        (circle.y - circle.r > rect.b.y) ||
        (circle.y + circle.r < rect.a.y));
}
//三角形
struct TRIANGLE {
    POINT a, b, c;
    TRIANGLE() {};
    TRIANGLE(const POINT & aa, const POINT & bb, const POINT & cc):
        a(aa), b(bb), c(cc) {}
};
//三角形标准化
TRIANGLE StdTRIANGLE(TYPE len1, TYPE len2, TYPE len3) //把3边长转化成三角形
{
    POINT a, b, c;     //已知边长后将其转化成坐标三角形
    a.x = a.y = 0.0;
    b.x = len1;
    b.y = 0.0;
    c.x = (len2 * len2 - len3 * len3 + len1 * len1)/len1/2.0;
    c.y = Sqrt(len2 * len2 - c.x * c.x);
    TRIANGLE temp(a,b,c);
    return temp;
};
//三角形费马点(即三角形到三顶点之和最小的点)
POINT fermentPONT(TRIANGLE &t) {
    POINT u, v;
    double step = fabs(t.a.x) + fabs(t.a.y) + fabs(t.b.x) + fabs(t.b.y) + fabs(t.c.x) + fabs(t.c.y);
    int i, j, k;
    u.x = (t.a.x + t.b.x + t.c.x)/3;
    u.y = (t.a.y + t.b.y + t.c.y)/3;
    while (step > EPS)
    for(k = 0; k < 10; step /= 2, ++k)
    for(i = -1; i <= 1; ++i)
    for(j = -1; j <= 1; ++j) {
        v.x = u.x + step * i;
        v.y = u.y + step * j;
        if(Distance(u,t.a) + Distance(u,t.b) + Distance(u,t.c) > Distance(v,t.a) + Distance(v,t.b) + Distance(v,t.c))
            u = v;
    }
    return u;
};
//三角形重心(中线交点)
POINT InCenter(const TRIANGLE &t) {
    return POINT((t.a.x + t.b.x + t.c.x)/3, (t.a.y + t.b.y + t.c.y)/3);
}
//三角形外心(中垂线交点)
POINT CcCenter(const TRIANGLE &t) {
    POINT u, v;
    LINE A, B;
    A.a = t.a; A.b = t.b; B.a = t.b; B.b = t.c;
    TYPE A1, B1, C1;
    TYPE A2, B2, C2;
    Coefficient(A, B1, A1, C1);
    Coefficient(B, B2, A2, C2);
    B1 = -B1; B2 = -B2;
    C1 = -((A.a.x + A.b.x) * A1 + (A.a.y + B.a.y) * B1)/2;
    C2 = -((B.a.x + B.b.x) * A2 + (B.a.y + B.b.y) * B2)/2;
    POINT I(0, 0);
    I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);
    I.y =   (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1);
    return I;
}
//三角形内心(角平分线交点)
POINT incenter(const TRIANGLE &t) {
    LINE u, v;
    TYPE m, n;
    u.a = t.a;
    m = atan2(t.b.y - t.a.y, t.b.x - t.a.x);
    n = atan2(t.c.y - t.a.y, t.c.x - t.a.x);
    u.b.x = u.a.x + cos((m + n)/2);
    u.b.y = u.a.y + sin((m + n)/2);
    v.a = t.b;
    m = atan2(t.a.y - t.b.y, t.a.x - t.b.x);
    n = atan2(t.c.y - t.b.y, t.c.x - t.b.x);
    v.b.x = v.a.x + cos((m + n)/2);
    v.b.y = v.a.y + sin((m + n)/2);
    return Intersection(u, v);
};
//三角形内切圆面积
TYPE InArea(const TRIANGLE &t) {
    TYPE p, a, b, c, s;
    a = Distance(t.a, t.b);
    c = Distance(t.a, t.c);
    b = Distance(t.c, t.b);
    s = (a + b + c)/2;
    p = s * (s - a) * (s - b) * (s - c);
    s = p * PI/s/s;
    return s;
}
//三角形外接圆面积
TYPE OutArea(const TRIANGLE & t) {
    TYPE a,b,c;
    a = (t.a.x - t.b.x) * (t.a.x - t.b.x) + (t.a.y - t.b.y) * (t.a.y - t.b.y);
    b = (t.a.x - t.c.x) * (t.a.x - t.c.x) + (t.a.y - t.c.y) * (t.a.y - t.c.y);
    c = (t.c.x - t.b.x) * (t.c.x - t.b.x) + (t.c.y - t.b.y) * (t.c.y - t.b.y);
    a = PI * sqrt(c/(1-(a+b-c)*(a+b-c)/a/b/4));
    return a;
}
// 多边形 ,逆时针或顺时针给出x,y
struct POLY {
    int n;     //n个点
    TYPE *x, *y;   //x,y为点的指针，首尾必须重合
    POLY(): n(0), x(NULL), y(NULL) {};
    POLY(int nn, const TYPE *xx, const TYPE *yy) {
        n = nn;
        x = new TYPE[n + 1];
        memcpy(x, xx, n*sizeof(TYPE));
        x[n] = xx[0];
        y = new TYPE[n + 1];
        memcpy(y, yy, n*sizeof(TYPE));
        y[n] = yy[0];
    }
};
//多边形顶点
POINT Vertex(const POLY &poly, int idx) {
  // idx %= poly.n;
  return POINT(poly.x[idx], poly.y[idx]);
}
//多边形的边
SEG Edge(const POLY &poly, int idx) {
  // idx %= poly.n;
  return SEG(POINT(poly.x[idx], poly.y[idx]), POINT(poly.x[idx + 1], poly.y[idx + 1]));
}
//多边形的周长
TYPE Perimeter(const POLY & poly) {
    TYPE p = 0.0;
    for (int i = 0; i < poly.n; ++i)
        p += Distance(Vertex(poly, i), Vertex(poly, i + 1));
    return p;
}
//求多边形面积
TYPE Area(const POLY & poly) {
    if (poly.n < 3) return TYPE(0);
    double s = poly.y[0] * (poly.x[poly.n - 1] - poly.x[1]);
    for (int i = 1; i < poly.n; ++i) {
        s += poly.y[i] * (poly.x[i - 1] - poly.x[(i + 1) % poly.n]);
    }
    return s/2;
}
//多边形的重心
POINT InCenter(const POLY & poly) {
    TYPE S,Si,Ax,Ay;
    POINT p;
    Si = (poly.x[poly.n - 1] * poly.y[0] - poly.x[0] * poly.y[poly.n - 1]);
    S = Si;
    Ax = Si * (poly.x[0] + poly.x[poly.n - 1]);
    Ay = Si * (poly.y[0] + poly.y[poly.n - 1]);
    for(int i = 1; i < poly.n; ++i){
        Si = (poly.x[i-1] * poly.y[i] - poly.x[i] * poly.y[i-1]);
        Ax += Si * (poly.x[i-1] + poly.x[i]);
        Ay += Si * (poly.y[i-1] + poly.y[i]);
        S += Si;
    }
    S = S * 3;
    return POINT(Ax/S,Ay/S);
}
//判断点是否在多边形上
bool IsOnPoly(const POLY &poly, const POINT &p) {
    for(int i = 0; i < poly.n; ++i){
        if( IsOnSeg(Edge(poly, i), p) ) return true;
    }
    return false;
}
//判断点是否在多边形内部
bool IsInPoly(const POLY &poly, const POINT &p) {
    SEG L(p, POINT(Infinity, p.y));
    int count = 0;
    for (int i = 0; i < poly.n; i++) {
    SEG S = Edge(poly, i);
    if (IsOnSeg(S, p)) {
        return true; //如果想让 在poly上则返回 true,
    }
    if(!IsEqual(S.a.y, S.b.y)) {
        POINT & q = (S.a.y > S.b.y)?(S.a):(S.b);
        if(IsOnSeg(L, q)) ++count;
        else if(!IsOnSeg(L, S.a) && !IsOnSeg(L, S.b) && IsIntersect(S, L)) {
          ++count;
        }
    }
  }
  return (count % 2 != 0);
}
//判断是否简单多边形
bool IsSimple(const POLY & poly) {
    if(poly.n < 3) return false;
    SEG L1, L2;
    for(int i = 0; i < poly.n - 1; ++i) {
        L1 = Edge(poly, i);
        for(int j = i + 1; j < poly.n; ++j) {
            L2 = Edge(poly, j);
            if(j == i + 1) {
                if(IsOnSeg(L1, L2.b) || IsOnSeg(L2, L1.a)) return false;
            }
            else if(j == poly.n - i - 1) {
              if (IsOnSeg(L1, L2.a) || IsOnSeg(L2, L1.b)) return false;
            }
            else if(IsIntersect(L1, L2)) return false;
        } // for j
    } // for i
    return true;
}
// 点阵的凸包，返回一个多边形(求法一，已测试)，不适用于点少于三个的情况
POLY ConvexHull(const POINT * set, int n) {
    POINT * points = new POINT[n];
    memcpy(points, set, n * sizeof(POINT));
    TYPE * X = new TYPE[n];
    TYPE * Y = new TYPE[n];
    int i, j, k = 0, top = 2;
    for(i = 1; i < n; ++i) {
        if ((points[i].y < points[k].y) || ((points[i].y == points[k].y) && (points[i].x < points[k].x)))
            k = i;
    }
    swap(points[0], points[k]);
    for(i = 1; i < n - 1; ++i) {
        k = i;
        for(j = i + 1; j < n; ++j) {
            if((Cross(points[j], points[k], points[0]) > 0) ||
              ((Cross(points[j], points[k], points[0]) == 0) &&
              (Distance(points[0], points[j]) < Distance(points[0], points[k]))))
                k = j;
        }
        swap(points[i], points[k]);
    }
    X[0] = points[0].x; Y[0] = points[0].y;
    X[1] = points[1].x; Y[1] = points[1].y;
    X[2] = points[2].x; Y[2] = points[2].y;
    for(i = 3; i < n; i++) {
        while(Cross(points[i], POINT(X[top], Y[top]),POINT(X[top - 1], Y[top - 1])) >= 0)
            --top;
        ++top;
        X[top] = points[i].x;
        Y[top] = points[i].y;
    }
    delete [] points;
    POLY poly(++top, X, Y);
    delete [] X;
    delete [] Y;
    return poly;
}
// 点阵的凸包，返回一个多边形 (求法二，未测试)
bool cmp(const POINT& aa, const POINT& bb) {
    if(aa.y != bb.y) return aa.y < bb.y;
    else return aa.x < bb.x;
}
POLY ConvexHull2(const POINT * set, int n) {// 不适用于点少于三个的情况
    POINT *a = new POINT[n];
    memcpy(a, set, n * sizeof(POINT));
    sort(a, a + n, cmp);
    TYPE * X = new TYPE[n];
    TYPE * Y = new TYPE[n];
    X[0] = a[0].x;
    Y[0] = a[0].y;
    X[1] = a[1].x;
    Y[1] = a[1].y;
    int i,j;
    for(i = 2, j = 2; i < n; ++i) {
        X[j] = a[i].x;
        Y[j] = a[i].y;
        while(((X[j]-X[j-1]) * (Y[j-1]-Y[j-2]) - (Y[j]-Y[j-1]) * (X[j-1]-X[j-2])) >= 0 && j > 1) {
            X[j-1] = X[j];
            Y[j-1] = Y[j];
            --j;
        }
        ++j;
    }
    X[j] = a[n-2].x;
    Y[j] = a[n-2].y;
    ++j;
    for(i = n - 3; i >= 0; --i) {
        X[j] = a[i].x;
        Y[j] = a[i].y;
        while(((X[j]-X[j-1]) * (Y[j-1]-Y[j-2]) - (Y[j]-Y[j-1]) * (X[j-1]-X[j-2])) >= 0) {
            X[j-1] = X[j];
            Y[j-1] = Y[j];
            --j;
        }
        ++j;
    }
    delete [] a;
    POLY poly(j, X, Y);
    delete [] X;
    delete [] Y;
    return poly;
}

int main() {
    double a,b,c,l1,l2,l3,l4;
    POINT d;
    int n;
    scanf("%d",&n);
    while(n--)
    {
        scanf("%lf%lf%lf",&a,&b,&c);
        TRIANGLE t=StdTRIANGLE(a,b,c);

        d=fermentPONT(t);//费马点
        l1=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);

        d=incenter(t);//内心
        l2=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);

        d=InCenter(t);//重心
        l3=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);

        d=CcCenter(t);//外心
        l4=Distance(d,t.a)+Distance(d,t.b)+Distance(d,t.c);

        printf("%.3lf %.3lf %.3lf %.3lf
",l1,l2,l3,l4);
    }
    return 0;
}

其他2：

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;

struct Tpoint{
    double x, y;
};

inline double length(Tpoint A, Tpoint B){
    double x = A.x - B.x, y = A.y - B.y;
    return sqrt(x * x + y * y);
}

inline double across(Tpoint A, Tpoint B, Tpoint C){
    double x1 = C.x - A.x, x2 = C.x - B.x;
    double y1 = C.y - A.y, y2 = C.y - B.y;
    return x1 * y2 - x2 * y1;
}

double distance(Tpoint A, Tpoint B, Tpoint C){
    return abs(across(A,B,C))/length(B,C);
}
//http://iask.sina.com.cn/b/14697945.html
int main(){
    Tpoint A, B, C, cen;
    while(true){
        cin>>A.x>>A.y>>B.x>>B.y>>C.x>>C.y;
        /*
        srand(time(0));
        A.x = (double)rand()/rand(); A.y = (double)rand()/rand();
        B.x = (double)rand()/rand(); B.y = (double)rand()/rand();
        C.x = (double)rand()/rand(); C.y = (double)rand()/rand();
        cout<<A.x<<' '<<A.y<<' '<<B.x<<' '<<B.y<<' '<<C.x<<' '<<C.y<<endl;*/
        double a = length(B,C), b = length(C, A), c = length(A, B);
        if(abs(across(A, B, C)) < 0.01) {cout<<"Wrong Input!"<<endl; continue;}
        cen.x = (a * A.x + b * B.x + c * C.x)/(a + b + c);
        cen.y = (a * A.y + b * B.y + c * C.y)/(a + b + c);
        cout<<cen.x<<' '<<cen.y<<endl;
        cout<<distance(cen, A, B)<<' '<<distance(cen, B, C)<<' '<<distance(cen, C, A)<<endl;
        //cout<<abs(across(A, B, C))/(a + b + c)<<endl;
        system("pause");
    }
}



#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;

struct Tpoint{
    double x, y, z;
};

inline double length(Tpoint A, Tpoint B){
    double x = A.x - B.x, y = A.y - B.y, z = A.z - B.z;
    return sqrt(x * x + y * y + z * z);
}

inline double across(Tpoint A, Tpoint B, Tpoint C){
    Tpoint p1, p2;
    p1.x = C.x - A.x; p1.y = C.y - A.y; p1.z = C.z - A.z;
    p2.x = C.x - B.x; p2.y = C.y - B.y; p2.z = C.z - B.z;
    double x = p1.y * p2.z - p1.z * p2.y;
    double y = p1.z * p2.x - p1.x * p2.z;
    double z = p1.x * p2.y - p1.y * p2.x;
    //cout<<x<<' '<<y<<' '<<z<<endl;
    return sqrt(x * x + y * y + z * z);
}

double distance(Tpoint A, Tpoint B, Tpoint C){
    return abs(across(A,B,C))/length(B,C);
}

void neijieyuan(Tpoint A, Tpoint B, Tpoint C){
    Tpoint cen;

    srand(time(0));
    A.x = (double)rand()/rand(); A.y = -(double)rand()/rand();
    B.x = -(double)rand()/rand(); B.y = (double)rand()/rand();
    C.x = -(double)rand()/rand(); C.y = (double)rand()/rand();
    cout<<A.x<<' '<<A.y<<' '<<B.x<<' '<<B.y<<' '<<C.x<<' '<<C.y<<endl;
    double a = length(B,C), b = length(C, A), c = length(A, B);
    //cout<<a<<' '<<b<<' '<<c<<endl;
    if(abs(across(A, B, C)) < 0.01) {cout<<"Wrong Input!"<<endl;return;}
    cen.x = (a * A.x + b * B.x + c * C.x)/(a + b + c);
    cen.y = (a * A.y + b * B.y + c * C.y)/(a + b + c);
    cen.z = (a * A.z + b * B.z + c * C.z)/(a + b + c);
    cout<<cen.x<<' '<<cen.y<<' '<<cen.z<<endl;
    cout<<distance(cen, A, B)<<' '<<distance(cen, B, C)<<' '<<distance(cen, C, A)<<endl;
}

int main(){
    Tpoint A, B, C;
    while(true){
        //cin>>A.x>>A.y>>A.z>>B.x>>B.y>>B.z>>C.x>>C.y>>C.z;
        //cout<<across(A, B, C)<<endl;
        neijieyuan(A, B, C);
        //cout<<abs(across(A, B, C))/(a + b + c)<<endl;
        system("pause");
    }
}



＃include <stdio.h>
＃include <math.h>

typedef struct TPoint

{

    //平面点

    double x;

    double y;

}TPoint;


typedef struct TTriangle

{

    TPoint t[2];

}TTriangle;

 

typedef struct TCircle

{

    //圆

    double r;

    TPoint centre;

}TCircle;

 

double distance(const TPoint p1, const TPoint p2)

{

    //计算平面上两个点之间的距离

   return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));   

}

 

double triangleArea(const TTriangle t)

{

    //已知三角形三个顶点的坐标，求三角形的面积

    return fabs(t.t[0].x * t.t[1].y + t.t[1].x * t.t[2].y + t.t[2].x * t.t[0].y

      - t.t[1].x * t.t[0].y - t.t[2].x * t.t[1].y - t.t[0].x * t.t[2].y) / 2;  

}

 

TCircle circumcircleOfTriangle(const TTriangle t)

{

    //三角形的外接圆

    TCircle tmp;

    double a, b, c, c1, c2;

    double xA, yA, xB, yB, xC, yC;

    a = distance(t.t[0], t.t[1]);

    b = distance(t.t[1], t.t[2]);

    c = distance(t.t[2], t.t[0]);

    printf("a = %lf b = %lf, c = %lf
", a, b, c);

    //根据S = a * b * c / R / 4;求半径R

    tmp.r = a * b * c / triangleArea(t) / 4;

   

   xA = t.t[0].x; yA = t.t[0].y;

    xB = t.t[1].x; yB = t.t[1].y;

    xC = t.t[2].x; yC = t.t[2].y;

    c1 = (xA * xA + yA * yA - xB * xB - yB * yB) / 2;

    c2 = (xA * xA + yA * yA - xC * xC - yC * yC) / 2;

   

    tmp.centre.x = (c1 * (yA - yC) - c2 * (yA - yB)) /

         (xA - xB) * (yA - yC) - (xA - xC) * (yA - yB);

    tmp.centre.y = (c1 * (xA - xC) - c2 * (xA - xB)) /

         (yA - yB) * (xA - xC) - (yA - yC) * (xA - xB);

        

    return tmp;    

}

 

TCircle incircleOfTriangle(const TTriangle t)

{

    //三角形的内接圆

    TCircle tmp;

    double a, b, c, angleA, angleB, angleC, p, p2, p3, f1, f2;

    double xA, yA, xB, yB, xC, yC;

    a = distance(t.t[0], t.t[1]);

    b = distance(t.t[1], t.t[2]);

    c = distance(t.t[2], t.t[0]);

   

    printf("a = %lf b = %lf, c = %lf
", a, b, c);

    tmp.r = 2 * triangleArea(t) / (a + b +c);

    angleA = acos((b * b + c * c - a * a) / (2 * b * c));

    angleB = acos((a * a + c * c - b * b) / (2 * a * c));

    angleC = acos((a * a + b * b - c * c) / (2 * a * b));

    p = sin(angleA / 2);

    p2 = sin(angleB / 2);

    p3 = sin(angleC / 2);

   

    xA = t.t[0].x; yA = t.t[0].y;

    xB = t.t[1].x; yB = t.t[1].y;

    xC = t.t[2].x; yC = t.t[2].y;

   

    f1 = ((tmp.r / p2) * (tmp.r / p2) - (tmp.r / p) * (tmp.r / p) +

         xA * xA - xB * xB + yA * yA - yB * yB) / 2;

    f2 = ((tmp.r / p3) * (tmp.r / p3) - (tmp.r / p) * (tmp.r / p) +

         xA * xA - xC * xC + yA * yA - yC * yC) / 2;

    tmp.centre.x = (f1 * (yA - yC) - f2 * (yA - yB)) /

                   ((xA - xB) * (yA - yC) - (xA - xC) * (yA - yB));

    tmp.centre.y = (f1 * (xA - xC) - f2 * (xA - xB)) /

                   ((yA - yB) * (xA - xC) - (yA - yC) * (xA - xB));

    return tmp;  

}

 

int main()

{

    TTriangle t;

    TCircle circumcircle, incircle;

    while(scanf("%lf%lf%lf%lf%lf%lf",

                    &t.t[0].x, &t.t[0].y, &t.t[1].x, &t.t[1].y,

                             &t.t[2].x, &t.t[2].y) != EOF)

    {

          incircle = incircleOfTriangle(t);

          circumcircle = circumcircleOfTriangle(t);

          printf("%lf %lf 
", incircle.r, circumcircle.r);

          printf("%.3lf
", incircle.r / circumcircle.r);

    }      

    return 0;

}