*模板--计算几何

二维几何相关（未完全测试）

基础+凸包+旋转卡壳+半平面交

/*************************************************************************
    > File Name: board.cpp
    > Author: yijiull
    > Mail: 1147161372@qq.com 
    > Created Time: 2017年09月22日 星期五 08时32分54秒
 ************************************************************************/
#include <iostream>
#include <cstring>
#include <cstdio>
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const int inf = 0x3f3f3f3f;
struct Point {
    double x,y;
    Point (double x = 0, double y = 0) : x(x), y(y) {}
};
typedef Point Vector;
Vector operator + (Vector a, Vector b) {
    return Vector (a.x + b.x, a.y + b.y);
}
Vector operator * (Vector a, double s) {
    return Vector (a.x * s, a.y * s);
}
Vector operator / (Vector a, double p) {
    return Vector (a.x / p, a.y / p);
}
Vector operator - (Point a, Point b) {
    return Vector (a.x - b.x, a.y - b.y);
}
bool operator < (Point a, Point b) {
    return a.x < b.x || (a.x == b.x && a.y < b.y);
}
int dcmp (double x) {
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}
bool operator == (const Point &a, const Point &b) {
    return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;
}
double Angel (Vector a) {
    return atan2(a.y, a.x);
}
double Dot(Vector a, Vector b) {
    return a.x * b.x + a.y * b.y;
}
double Length (Vector a) {
    return sqrt(Dot(a, a));
}
double Angle (Vector a, Vector b) {
    return acos(Dot(a, b) / Length(a) / Length(b));
}
double Cross (Vector a, Vector b) {
    return a.x * b.y - a.y * b.x;
}
double  Area2 (Point a, Point b, Point c) {
    return Cross(b - a, c - a);
}
Vector Rotate (Vector a, double rad) {  //右手逆时针旋转rad(弧度)
    return Vector (a.x * cos(rad) - a.y * sin(rad), a.x * sin(rad) + a.y * cos(rad));
}
Vector Normal (Vector a) {
    double L = Length(a);
    return Vector(-a.y / L, a.x / L);
}
//两直线交点
Point GetLineIntersection (Point p, Vector v, Point q, Vector w) {
    Vector u = p - q;
    double t1 = Cross(w, u) / Cross(v, w);
    double t2 = Cross(v, u) / Cross(v, w);
    return p + v * t1;  //  return q + w * t2;
}
//点到直线的dis
double DistanceToLine(Point p, Point a, Point b) {
    Vector v1 = b - a, v2 = p - a;
    return fabs(Cross(v1, v2)) / Length(v1);
}
//点到线段的dis
double DistanceToSegment(Point p, Point a, Point b) {
    if(a == b) return Length(a - p);
    Vector v1 = b - a, v2 = p - a, v3 = p - b;
    if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
    else if(dcmp(Dot(v1, v3)) > 0) return Length(v3); 
    else return fabs(Cross(v1, v2)) / Length(v1);
}
//点在直线的投影
Point GetLineProjection(Point p, Point a, Point b) {
    Vector v = b - a;
    return a + v * (Dot(v, p-a) / Dot(v, v));
}
//判断两线段是否规范相交
bool SegmentProperIntersection(Point a1, Point b1, Point a2, Point b2) { 
    double c1 = Cross(b1 - a1, a2 - a1), c2 = Cross(b1 - a1, b2 - a1),
           c3 = Cross(b2 - a2, a1 - a2), c4 = Cross(b2 - a2, b1 - a2);
    return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0;
}
//判断点是否在线段上(即∠APB等于π)
bool OnSegment(Point p, Point a, Point b) {
    return dcmp(Cross(a - p, b - p)) == 0 && dcmp(Dot(a - p, b - p)) < 0;
}
//多边形面积
double ConvexPolygonArea(Point *p, int n) {
    double area = 0;
    for(int i = 1; i < n-1; i++) {
        area += Cross(p[i] - p[0], p[i+1] - p[0]);
    }
    return area / 2;
}
//判断点是否在多边形内or外or上
bool isPointInPolygon(Point p,Point *poly, int n) {
    int cnt = 0;
    poly[n] = poly[0];
    for(int i = 0; i < n; i++) {
        if(OnSegment(p, poly[i], poly[i+1])) return 0;  //在边界上
        int k = dcmp(Cross(poly[i+1]-poly[i], p - poly[i]));
        int d1 = dcmp(poly[i].y - p.y);
        int d2 = dcmp(poly[i+1].y - p.y);
        if(k > 0 && d1 <= 0 && d2 > 0) cnt++;
        if(k < 0 && d2 <= 0 && d1 >= 0) cnt--;
    }
    if(cnt) return -1; //内部
    return 1; //外部
}
//凸包
int ConvexHull(Point *p, int n, Point *ch) {
    int m = 0;
    sort(p, p+n);
    for(int i = 0; i < n; i++) {
        while(m > 1 && Cross(ch[m-1] - ch[m-2], p[i] - ch[m-2]) <= 0) m--;  //去掉在边上的点
        ch[m++] = p[i];
    }
    int k = m;
    for(int i = n-2; i >= 0; i--) {
        while(m > k && Cross(ch[m-1] - ch[m-2], p[i] - ch[m-2]) <= 0) m--;  //
        ch[m++] = p[i];
    }
    if(m > 1) m--;
    return m;
}
//判断是否是稳定凸包
bool IsStableConvexHull(Point *ch, int n) {
    if(n == 0) return 0;
    ch[n] = ch[0];
    ch[n+1] = ch[1];
    for(int i = 0; i < n; i++) {
        if(Cross(ch[i] - ch[(i-1+n)%n], ch[i+1] - ch[i]) != 0 && Cross(ch[i+1] - ch[i] , ch[i+2] - ch[i]) != 0) 
            return 0;
    }
    return 1;
}
//旋转卡壳求凸包直径
double RotateCalipers(Point *ch, int n) {
    double ans = - inf;
    ch[n] = ch[0];
    int q = 1;
    for(int i = 0; i < n; i++) {
        while(Cross(ch[i+1] - ch[i], ch[q] - ch[i]) < Cross(ch[i+1] - ch[i], ch[q+1] -ch[i])) q = (q+1)%n;
        ans = max(ans, max(Length(ch[i] - ch[q]), Length(ch[i+1] -ch[q+1])));
    }
    return ans;
}

//半平面交
struct Line{
    Point p;
    Vector v;
    double rad;
    Line () {}
    Line (Point p, Vector v) : p(p), v(v) {
        rad = atan2(v.y,v.x);
    }
    bool operator < (const Line &L) const {
        return rad < L.rad;
    }
};
bool OnLeft(Line L, Point p) {
    return Cross(L.v, p - L.p) > 0;
}
Point GetLineIntersection (Line a, Line b) {
    Vector u = a.p - b.p;
    double t = Cross(b.v, u) / Cross(a.v, b.v);
    return a.p + a.v*t;
}

int HalfplaneIntersection(Line *L, int n,Point *poly) {
    sort(L, L+n);
    int first,last;
    Point *p = new Point[n];
    Line *q = new Line[n];  //双端队列
    q[first = last = 0] = L[0];
    for(int i = 1; i < n; i++) {
        while(first < last && !OnLeft(L[i], p[last-1])) last--;   //去尾
        while(first < last && !OnLeft(L[i], p[first])) first++; 
        q[++last] = L[i];
        if(dcmp(Cross(q[last].v, q[last-1].v)) == 0) {
            last--;
            if(OnLeft(q[last], L[i].p)) q[last] = L[i];
        }
        if(first < last) p[last-1] = GetLineIntersection (q[last-1],q[last]);
    }
    while(first < last && !OnLeft(q[first], p[last-1])) last--;  //删除无用平面
    if(last - first <= 1) return 0;  //空集
    p[last] = GetLineIntersection (q[last], q[first]);

    int m = 0;
    for(int i = first; i <= last; i++) poly[m++] = p[i];
    return m;
}

int main(){

}

//线段AB与圆的交点
void LineIntersectionCircle(Point A, Point B, Circle cr, Point *p, int &num){
    double x0 = cr.o.x, y0 = cr.o.y;
    double x1 = A.x, y1 = A.y, x2 = B.x, y2 = B.y;
    double dx = x2 - x1, dy = y2 - y1;
    double a = dx * dx + dy * dy;
    double b = 2 * dx *(x1 - x0) + 2 * dy * (y1 -y0);
    double c = (x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0) - cr.r *cr.r;
    double delta = b*b - 4*a*c;
    num = 0;
    if(dcmp(delta) >= 0){
        double t1 = (-b - sqrt(delta)) / (2*a);
        double t2 = (-b + sqrt(delta)) / (2*a);
        if(dcmp(t1-1) <= 0 && dcmp(t1) >= 0) {
            p[num++] = Point(x1 + t1*dx, y1 + t1*dy);
        }
        if(dcmp(t2-1) <= 0 && dcmp(t2) >= 0) {
            p[num++] = Point(x1 + t2*dx, y1 + t2*dy);
        }
    }
}
//扇形面积
double SectorArea(Point a, Point b, Circle cr){
    double theta = Angle(a - cr.o) - Angle(b - cr.o);
    while(theta <= 0) theta += 2*pi;
    while(theta > 2*pi) theta -= 2*pi;
    theta = min(theta, 2*pi - theta);
    return cr.r * cr.r * theta / 2;
}

double cal(Point a, Point b, Circle cr){
    Point tp[3];
    int num = 0;
    Vector ta = a - cr.o;
    Vector tb = b - cr.o;
    bool ina = (Length(ta) - cr.r) < 0;
    bool inb = (Length(tb) - cr.r) < 0;
    if(ina){
        if(inb){
            return fabs(Cross(ta, tb))/2;
        } else{
            LineIntersectionCircle(a, b, cr, tp, num); 
            return SectorArea(b, tp[0], cr) + fabs(Cross(ta, tp[0] - cr.o))/2;
        }
    } else{
        if(inb){
            LineIntersectionCircle(a, b, cr, tp, num); 
            return SectorArea(a, tp[0], cr) + fabs(Cross(tb, tp[0] - cr.o))/2;
        } else {
            LineIntersectionCircle(a, b, cr, tp, num); 
            if(num == 2) {
                return SectorArea(a, tp[0], cr) + SectorArea(b, tp[1], cr) + fabs(Cross(tp[0]-cr.o, tp[1]-cr.o))/2;
            } else {
                return SectorArea(a, b, cr);
            }
        }
    }
}
//圆与多边形的面积并
double CirclePolyArea(Point *p, int n, Circle cr){
   double res = 0;
   p[n] = p[0];
   for(int i = 0; i < n; i++){
       int sgn = dcmp(Cross(p[i] - cr.o, p[i+1] - cr.o));
       if(sgn){
           res += sgn * cal(p[i], p[i+1], cr);
       }
   }
   return res;
}

自适应辛普森积分

double F(double x) {
    //Simpson公式用到的函数
}
double simpson(double a, double b) {  //三点Simpson法，这里要求F是一个全局函数
    double c = a + (b - a) / 2;
    return (F(a) + 4 * F(c) + F(b))*(b - a) / 6;
}
double asr(double a, double b, double eps, double A) { //自适应Simpson公式（递归过程）。已知整个区间[a,b]上的三点Simpson值A
    double c = a + (b - a) / 2;
    double L = simpson(a, c), R = simpson(c, b);
    if (fabs(L + R - A) <= 15 * eps)return L + R + (L + R - A) / 15.0;
    return asr(a, c, eps / 2, L) + asr(c, b, eps / 2, R);
}
double asr(double a, double b, double eps) { //自适应Simpson公式（主过程）
    return asr(a, b, eps, simpson(a, b));
}

多边形重心

/*
 *  求多边形重心
 *  INIT: pnt[]已按顺时针(或逆时针)排好序; | CALL: res = bcenter(pnt, n);
 */
Point bcenter(Point pnt[], int n) {
    Point p, s;
    double tp, area = 0, tpx = 0, tpy = 0;
    p.x = pnt[0].x;
    p.y = pnt[0].y;
    for (int i = 1; i <= n; ++i) {
        //  Point:0 ~ n - 1
        s.x = pnt[(i == n) ? 0 : i].x;
        s.y = pnt[(i == n) ? 0 : i].y;
        tp = (p.x * s.y - s.x * p.y);
        area += tp / 2;
        tpx += (p.x + s.x) * tp;
        tpy += (p.y + s.y) * tp;
        p.x = s.x;
        p.y = s.y;
    }
    s.x = tpx / (6 * area);
    s.y = tpy / (6 * area);
    return s;
}

相关公式（转自f_zyj）

已知圆锥表面积S求最大体积V

    V = S * sqrt(S / (72 * Pi))

三角形重点

设三角形的三条边为a, b, c, 且不妨假设a <= b <= c.
面积

三角形面积可以根据海伦公式求得：

    s = sqrt(p * (p - a) * (p - b) * (p - c));
    p = (a + b + c) / 2;

关键点与A, B, C三顶点距离之和
费马点

该点到三角形三个顶点的距离之和最小。
有个有趣的结论:
若三角形的三个内角均小于120度,那么该点连接三个顶点形成的三个角均为120度;若三角形存在一个内角大于120度,则该顶点就是费马点。
计算公式如下:
若有一个内角大于120度(这里假设为角C),则距离为a + b;若三个内角均小于120度,则距离为sqrt((a * a + b * b + c * c + 4 * sqrt(3.0) * s) / 2)。
内心

角平分线的交点。
令x = (a + b - c) / 2, y = (a - b + c) / 2, z = (-a + b + c) / 2, h = s / p.
计算公式为sqrt(x * x + h * h) + sqrt(y * y + h * h) + sqrt(z * z + h * h)。
重心

中线的交点。
计算公式如下:
2.0 / 3 * (sqrt((2 * (a * a + b * b) - c * c) / 4)
+ sqrt((2 * (a * a + c * c) - b * b) / 4) + sqrt((2 * (b * b + c * c) - a * a) / 4))。
垂心

垂线的交点。
计算公式如下:
3 * (c / 2 / sqrt(1 - cosC * cosC))。
外心

三点求圆心坐标。

Point waixin(Point a, Point b, Point c)
{
    double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1 * a1 + b1 * b1) / 2;
    double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2 * a2 + b2 * b2) / 2;
    double d = a1 * b2 - a2 * b1;
    return Point(a.x + (c1 * b2 - c2 * b1) / d, a.y + (a1 * c2 -a2 * c1) / d);
}

圆的面积交

double area[maxn];
#define sqr(x) (x)*(x)
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[maxn], tp[maxn * 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]);
        }
    }
}

三维几何相关

/*************************************************************************
    > File Name: board.cpp
    > Author: yijiull
    > Mail: 1147161372@qq.com 
    > Created Time: 2017年09月22日 星期五 16时51分54秒
 ************************************************************************/
#include <iostream>
#include <cstring>
#include <cstdio>
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const int inf = 0x3f3f3f3f;
struct Point3{
    double x, y, z;
    Point3(double x=0, double y=0, double z=0) : x(x), y(y), z(z) {} 
};
int dcmp(double x) {
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1; 
}
typedef Point3 Vector3;
Vector3 operator + (Vector3 a, Vector3 b) {
    return Vector3(a.x + b.x, a.y + b.y, a.z + b.z);
}
Vector3 operator - (Point3 a, Point3 b) {
    return Vector3(a.x - b.x, a.y - b.y, a.z - b.z); 
}
Vector3 operator * (Vector3 a, double p) {
    return Vector3(a.x * 3, a.y * 3, a.z * 3);
}
Vector3 operator / (Vector3 a, double p) {
    return Vector3(a.x / p, a.y / p, a.z / p);
}
double Dot(Vector3 a, Vector3 b) {
    return a.x * b.x + a.y * b.y + a.z * b.z; 
}
bool operator == (Point3 a, Point3 b) {
    return a.x==b.x && a.y==b.y && a.z==b.z;
}
double Length(Vector3 a) {
    return sqrt(Dot(a,a));
}
double Angle(Vector3 a, Vector3 b) {
    return acos(Dot(a,b) / Length(a) / Length(b)) ; 
}
Vector3 Normal(Vector3 a) {
    return a / Length(a);
}
//点p到平面p0-n的距离,n是单位向量
double DistanceToPlane(const Point3 & p, const Point3 &p0, const Vector3 &n) {
    return fabs(Dot(p-p0, n));  //如果不取绝对值，得到的是有向距离
}

//点P在平面p0-n上的投影，ｎ必须为单位向量
Point3 GetPlaneProjection(const Point3 p, const Point3 &p0, const Vector3 &n) {
    return p - n * Dot(p - p0, n); 
}
//直线p1-p2到平面p0-n的交点，假定交点唯一存在
Point3 LinePlaneIntersection(Point3 p1, Point3 p2, Point3 p0, Vector3 n) {
    Vector3 v = p2 - p1;
    double t = (Dot(n, p0 -p1) / Dot(n, p2 -p1)); //判断分母是否为０
    return p1 + v * t; //如果是线段，判断ｔ是不是在０和１之间
}

Vector3 Cross(Vector3 a, Vector3 b) {
    return Vector3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x); 
}
double Area2(Point3 a, Point3 b, Point3 c) {
    return Length(Cross(b-a, c-a));
}
//点p在三角形p0,p1,p2中
bool PointInTri(Point3 p, Point3 p0, Point3 p1, Point3 p2) {
    double area1 = Area2(p, p0, p1);
    double area2 = Area2(p, p1, p2);
    double area3 = Area2(p, p2, p0);
    return dcmp(area1 + area2 + area3 - Area2(p0, p1, p2)) == 0;
}
//三角形p0,p1,p2是否和线段ab相交
bool TriSegIntersection(Point3 p0, Point3 p1, Point3 p2, Point3 a, Point3 b, Point3 &p) {
    Vector3 n = Cross(p1 - p0, p2 - p0);
    if(dcmp(Dot(n, b - a)) == 0) return false;  //无交点
    double t = Dot(n, p0 - a) / Dot(n, b - a); 
    if(dcmp(t) < 0 || dcmp(t-1) > 0) return false;  //交点不在直线上
    p = a + (b - a) * t;
    return PointInTri(p, p0, p1, p2);  //是否在三角形内
}
//点p到直线ab的距离
double DistanceToLine(Point3 p, Point3 a, Point3 b) {
    Vector3 v1 = b - a, v2 = p - a;
    return Length(Cross(v1, v2)) / Length(v1);
}
//点p到线段ab的距离
double DisstanceToSegment(Point3 p, Point3 a, Point3 b) {
    if(a == b) return Length(p - a);
    Vector3 v1 = b - a, v2 = p - a, v3 = p - b;
    if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
    else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
    else return Length(Cross(v1, v2)) / Length(v1);
}
//四面体abcd的体积
double Volume6(Point3 a, Point3 b, Point3 c, Point3 d) {
    return Dot(d-a, Cross(b-a, c-a));
}

//三维凸包
int vis[200][200];
struct Face{
    int v[3];
    Vector3 normal(Point3 *p) const {
        return Cross(p[v[1]] - p[v[0]], p[v[2]] - p[v[0]]);
    }
    int cansee(Point3 *p, int i) const {
        return Dot(p[i] - p[v[0]], normal(p)) > 0 ? 1 : 0;
    }
};
//倍增法求三维凸包
//没有考虑特殊情况(如四点共面),实践中,请在调用之前对输入点进行微小扰动
vector<Face> CH3D(Point3 *p, int n) {
    vector<Face> cur;
    memset(vis, 0, sizeof(vis));
    //由于已经扰动,前三个点不共线
    cur.push_back((Face){ {0, 1, 2} } );
    cur.push_back((Face){ {2, 1, 0} } );
    for(int i = 3; i < n; i++) {
        vector<Face> next;
        for(int j = 0; j < cur.size(); j++) {
            Face &f = cur[j];
            int res = f.cansee(p, i);
            if(!res) next.push_back(f);
            for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;
        }
        for(int j = 0; j < cur.size(); j++) {
            for(int k = 0; k < 3; k++) {
                int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
                if(vis[a][b] != vis[b][a] && vis[a][b])  //(a, b) 是分界线, 左边对(a, b) 可见
                    next.push_back((Face) { {a, b, i} } );
            }
        }
        cur = next;
    }
    return cur;
}
//扰动!!!
double rand01() {
    return rand() / (double)RAND_MAX; 
}
double randeps() {
    return (rand01() - 0.5) * eps;
}
Point3 add_noise(Point3 p) {
    return Point3(p.x + randeps(), p.y + randeps(), p.z + randeps());
}

int main(){

}