- /***********三分法求函数极值*************/
- void solve()
- {
- double L, R, m, mm, mv, mmv;
- while (L + eps < R)
- {
- m = (L + R) / 2;
- mm = (m + R) / 2;
- mv = calc(m);
- mmv = calc(mm);
- if (mv <= mmv) R = mm; //三分法求最大值时改为mv>=mmv
- else L = m;
- }
- }
- /*************基础***********/
- int dcmp(double x) {
- if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
- }
- 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 - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); }
- Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
- Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }
- bool operator < (const Point& a, const Point& b) {
- return a.x < b.x || (a.x == b.x && a.y < b.y);
- }
- bool operator == (const Point& a, const Point &b) {
- return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
- }
- 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 angle(Vector v) { return atan2(v.y, v.x); }
- double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; }
- Vector vecunit(Vector x){ return x / Length(x);} //单位向量
- Vector Normal(Vector x) { return Point(-x.y, x.x) / Length(x);} //垂直法向量
- Vector Rotate(Vector A, double rad) {
- return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
- }
- double Area2(const Point A, const Point B, const Point C) { return Length(Cross(B-A, C-A)); }
- /****************直线与线段**************/
- //求直线p+tv和q+tw的交点 Cross(v, w) == 0无交点
- Point GetLineIntersection(Point p, Vector v, Point q, Vector w)
- {
- Vector u = p-q;
- double t = Cross(w, u) / Cross(v, w);
- return p + v*t;
- }
- //点p在直线ab的投影
- Point GetLineProjection(Point P, Point A, Point B) {
- Vector v = B-A;
- return A+v*(Dot(v, P-A) / Dot(v, v));
- }
- //点到直线距离
- double DistanceToLine(Point P, Point A, Point B) {
- Vector v1 = B - A, v2 = P - A;
- return fabs(Cross(v1, v2)) / Length(v1); // 如果不取绝对值,得到的是有向距离
- }
- //点在p线段上
- bool OnSegment(Point p, Point a1, Point a2) {
- return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;
- }
- // 过两点p1, p2的直线一般方程ax+by+c=0
- // (x2-x1)(y-y1) = (y2-y1)(x-x1)
- void getLineGeneralEquation(const Point& p1, const Point& p2, double& a, double& b, double &c) {
- a = p2.y-p1.y;
- b = p1.x-p2.x;
- c = -a*p1.x - b*p1.y;
- }
- //点到线段距离
- double DistanceToSegment(Point p, Point a, Point b)
- {
- if(a == b) return Length(p-a);
- 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);
- }
- //两线段最近距离
- double dis_pair_seg(Point p1, Point p2, Point p3, Point p4)
- {
- return min(min(DistanceToSegment(p1, p3, p4), DistanceToSegment(p2, p3, p4)),
- min(DistanceToSegment(p3, p1, p2), DistanceToSegment(p4, p1, p2)));
- }
- //线段相交判定
- bool SegmentItersection(Point a1, Point a2, Point b1, Point b2)
- {
- double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1),
- c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);
- return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0;
- }
- // 有向直线。它的左边就是对应的半平面
- struct Line {
- Point P; // 直线上任意一点
- Vector v; // 方向向量
- double ang; // 极角,即从x正半轴旋转到向量v所需要的角(弧度)
- Line() {}
- Line(Point P, Vector v):P(P),v(v){ ang = atan2(v.y, v.x); }
- bool operator < (const Line& L) const {
- return ang < L.ang;
- }
- };
- //两直线交点
- Point GetLineIntersection(Line a, Line b) {
- return GetLineIntersection(a.p, a.v, b.p, b.v);
- }
- // 点p在有向直线L的左边(线上不算)
- bool OnLeft(const Line& L, const Point& p) {
- return Cross(L.v, p-L.P) > 0;
- }
- // 二直线交点,假定交点惟一存在
- Point GetLineIntersection(const Line& a, const 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;
- }
- // 半平面交主过程
- vector<Point> HalfplaneIntersection(vector<Line> L) {
- int n = L.size();
- sort(L.begin(), L.end()); // 按极角排序
- int first, last; // 双端队列的第一个元素和最后一个元素的下标
- vector<Point> p(n); // p[i]为q[i]和q[i+1]的交点
- vector<Line> q(n); // 双端队列
- vector<Point> ans; // 结果
- q[first=last=0] = L[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(fabs(Cross(q[last].v, q[last-1].v)) < eps) { // 两向量平行且同向,取内侧的一个
- 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 ans; // 空集
- p[last] = GetLineIntersection(q[last], q[first]); // 计算首尾两个半平面的交点
- // 从deque复制到输出中
- for(int i = first; i <= last; i++) ans.push_back(p[i]);
- return ans;
- }
- /***********多边形**************/
- //多边形有向面积
- double PolygonArea(vector<Point> p) {
- int n = p.size();
- 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;
- }
- // 点集凸包
- // 如果不希望在凸包的边上有输入点,把两个 <= 改成 <
- // 注意:输入点集会被修改
- vector<Point> ConvexHull(vector<Point>& p) {
- // 预处理,删除重复点
- sort(p.begin(), p.end());
- p.erase(unique(p.begin(), p.end()), p.end());
- int n = p.size();
- int m = 0;
- vector<Point> ch(n+1);
- 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(n > 1) m--;
- ch.resize(m);
- return ch;
- }
- // 凸包直径返回 点集直径的平方
- int diameter2(vector<Point>& points) {
- vector<Point> p = ConvexHull(points);
- int n = p.size();
- if(n == 1) return 0;
- if(n == 2) return Dist2(p[0], p[1]);
- p.push_back(p[0]); // 免得取模
- int ans = 0;
- for(int u = 0, v = 1; u < n; u++) {
- // 一条直线贴住边p[u]-p[u+1]
- for(;;) {
- int diff = Cross(p[u+1]-p[u], p[v+1]-p[v]);
- if(diff <= 0) {
- ans = max(ans, Dist2(p[u], p[v])); // u和v是对踵点
- if(diff == 0) ans = max(ans, Dist2(p[u], p[v+1])); // diff == 0时u和v+1也是对踵点
- break;
- }
- v = (v + 1) % n;
- }
- }
- return ans;
- }
- //两凸包最近距离
- double RC_Distance(Point *ch1, Point *ch2, int n, int m)
- {
- int q=0, p=0;
- REP(i, n) if(ch1[i].y-ch1[p].y < -eps) p=i;
- REP(i, m) if(ch2[i].y-ch2[q].y > eps) q=i;
- ch1[n]=ch1[0]; ch2[m]=ch2[0];
- double tmp, ans=1e100;
- REP(i, n)
- {
- while((tmp = Cross(ch1[p+1]-ch1[p], ch2[q+1]-ch1[p]) - Cross(ch1[p+1]-ch1[p], ch2[q]- ch1[p])) > eps)
- q=(q+1)%m;
- if(tmp < -eps) ans = min(ans,DistanceToSegment(ch2[q],ch1[p],ch1[p+1]));
- else ans = min(ans,dis_pair_seg(ch1[p],ch1[p+1],ch2[q],ch2[q+1]));
- p=(p+1)%n;
- }
- return ans;
- }
- //凸包最大内接三角形
- double RC_Triangle(Point* res,int n)// 凸包最大内接三角形
- {
- if(n<3) return 0;
- double ans=0, tmp;
- res[n] = res[0];
- int j, k;
- REP(i, n)
- {
- j = (i+1)%n;
- k = (j+1)%n;
- while((j != k) && (k != i))
- {
- while(Cross(res[j] - res[i], res[k+1] - res[i]) > Cross(res[j] - res[i], res[k] - res[i])) k= (k+1)%n;
- tmp = Cross(res[j] - res[i], res[k] - res[i]);if(tmp > ans) ans = tmp;
- j = (j+1)%n;
- }
- }
- return ans;
- }
- //模拟退火求费马点 保存在ptres中
- double fermat_point(Point *pt, int n, Point& ptres)
- {
- Point u, v;
- double step = 0.0, curlen, explen, minlen;
- int i, j, k, idx;
- bool flag;
- u.x = u.y = v.x = v.y = 0.0;
- REP(i, n)
- {
- step += fabs(pt[i].x) + fabs(pt[i].y);
- u.x += pt[i].x;
- u.y += pt[i].y;
- }
- u.x /= n;
- u.y /= n;
- flag = 0;
- 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;
- curlen = explen = 0.0;
- REP(idx, n)
- {
- curlen += dist(u, pt[idx]);
- explen += dist(v, pt[idx]);
- }
- if(curlen > explen)
- {
- u = v;
- minlen = explen;
- flag = 1;
- }
- }
- }
- ptres = u;
- return flag ? minlen : curlen;
- }
- //最近点对
- bool cmpxy(const Point& a, const Point& b)
- {
- if(a.x != b.x)
- return a.x < b.x;
- return a.y < b.y;
- }
- bool cmpy(const int& a, const int& b)
- {
- return point[a].y < point[b].y;
- }
- double Closest_Pair(int left, int right)
- {
- double d = INF;
- if(left==right)
- return d;
- if(left + 1 == right)
- return dis(left, right);
- int mid = (left+right)>>1;
- double d1 = Closest_Pair(left,mid);
- double d2 = Closest_Pair(mid+1,right);
- d = min(d1,d2);
- int i,j,k=0;
- //分离出宽度为d的区间
- for(i = left; i <= right; i++)
- {
- if(fabs(point[mid].x-point[i].x) <= d)
- tmpt[k++] = i;
- }
- sort(tmpt,tmpt+k,cmpy);
- //线性扫描
- for(i = 0; i < k; i++)
- {
- for(j = i+1; j < k && point[tmpt[j]].y-point[tmpt[i]].y<d; j++)
- {
- double d3 = dis(tmpt[i],tmpt[j]);
- if(d > d3)
- d = d3;
- }
- }
- return d;
- }
- /************圆************/
- struct Circle
- {
- Point c;
- double r;
- Circle(){}
- Circle(Point c, double r):c(c), r(r){}
- Point point(double a) //根据圆心角求点坐标
- {
- return Point(c.x+cos(a)*r, c.y+sin(a)*r);
- }
- };
- //直线与圆交点 返回个数
- int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector<Point>& sol){
- double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
- double e = a*a + c*c, f = 2*(a*b + c*d), g = b*b + d*d - C.r*C.r;
- double delta = f*f - 4*e*g; // 判别式
- if(dcmp(delta) < 0) return 0; // 相离
- if(dcmp(delta) == 0) { // 相切
- t1 = t2 = -f / (2 * e); sol.push_back(L.point(t1));
- return 1;
- }
- // 相交
- t1 = (-f - sqrt(delta)) / (2 * e); sol.push_back(L.point(t1));
- t2 = (-f + sqrt(delta)) / (2 * e); sol.push_back(L.point(t2));
- return 2;
- }
- //两圆交点 返回个数
- int getCircleCircleIntersection(Circle C1, Circle C2, vector<Point>& sol) {
- double d = Length(C1.c - C2.c);
- if(dcmp(d) == 0) {
- if(dcmp(C1.r - C2.r) == 0) return -1; // 重合,无穷多交点
- return 0;
- }
- if(dcmp(C1.r + C2.r - d) < 0) return 0;
- if(dcmp(fabs(C1.r-C2.r) - d) > 0) return 0;
- double a = angle(C2.c - C1.c);
- double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2*C1.r*d));
- Point p1 = C1.point(a-da), p2 = C1.point(a+da);
- sol.push_back(p1);
- if(p1 == p2) return 1;
- sol.push_back(p2);
- return 2;
- }
- //三角形外接圆
- Circle CircumscribedCircle(Point p1, Point p2, Point p3) {
- double Bx = p2.x-p1.x, By = p2.y-p1.y;
- double Cx = p3.x-p1.x, Cy = p3.y-p1.y;
- double D = 2*(Bx*Cy-By*Cx);
- double cx = (Cy*(Bx*Bx+By*By) - By*(Cx*Cx+Cy*Cy))/D + p1.x;
- double cy = (Bx*(Cx*Cx+Cy*Cy) - Cx*(Bx*Bx+By*By))/D + p1.y;
- Point p = Point(cx, cy);
- return Circle(p, Length(p1-p));
- }
- //三角形内切圆
- Circle InscribedCircle(Point p1, Point p2, Point p3) {
- double a = Length(p2-p3);
- double b = Length(p3-p1);
- double c = Length(p1-p2);
- Point p = (p1*a+p2*b+p3*c)/(a+b+c);
- return Circle(p, DistanceToLine(p, p1, p2));
- }
- // 过点p到圆C的切线。v[i]是第i条切线的向量。返回切线条数
- int getTangents(Point p, Circle C, Vector* v) {
- Vector u = C.c - p;
- double dist = Length(u);
- if(dist < C.r) return 0;
- else if(dcmp(dist - C.r) == 0) { // p在圆上,只有一条切线
- v[0] = Rotate(u, PI/2);
- return 1;
- } else {
- double ang = asin(C.r / dist);
- v[0] = Rotate(u, -ang);
- v[1] = Rotate(u, +ang);
- return 2;
- }
- }
- //所有经过点p 半径为r 且与直线L相切的圆心
- vector<Point> CircleThroughPointTangentToLineGivenRadius(Point p, Line L, double r) {
- vector<Point> ans;
- double t1, t2;
- getLineCircleIntersection(L.move(-r), Circle(p, r), t1, t2, ans);
- getLineCircleIntersection(L.move(r), Circle(p, r), t1, t2, ans);
- return ans;
- }
- //半径为r 与a b两直线相切的圆心
- vector<Point> CircleTangentToLinesGivenRadius(Line a, Line b, double r) {
- vector<Point> ans;
- Line L1 = a.move(-r), L2 = a.move(r);
- Line L3 = b.move(-r), L4 = b.move(r);
- ans.push_back(GetLineIntersection(L1, L3));
- ans.push_back(GetLineIntersection(L1, L4));
- ans.push_back(GetLineIntersection(L2, L3));
- ans.push_back(GetLineIntersection(L2, L4));
- return ans;
- }
- //与两圆相切 半径为r的所有圆心
- vector<Point> CircleTangentToTwoDisjointCirclesWithRadius(Circle c1, Circle c2, double r) {
- vector<Point> ans;
- Vector v = c2.c - c1.c;
- double dist = Length(v);
- int d = dcmp(dist - c1.r -c2.r - r*2);
- if(d > 0) return ans;
- getCircleCircleIntersection(Circle(c1.c, c1.r+r), Circle(c2.c, c2.r+r), ans);
- return ans;
- }
- //多边形与圆相交面积
- Point GetIntersection(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;
- }
- bool OnSegment(Point p, Point a1, Point a2)
- {
- return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;
- }
- bool InCircle(Point x, Circle c) { return dcmp(sqr(c.r) - sqr(Length(c.c - x))) >= 0;}
- bool OnCircle(Point x, Circle c) { return dcmp(sqr(c.r) - sqr(Length(c.c - x))) == 0;}
- //线段与圆的交点
- int getSegCircleIntersection(Line L, Circle C, Point* sol)
- {
- Vector nor = normal(L.v);
- Line pl = Line(C.c, nor);
- Point ip = GetIntersection(pl, L);
- double dis = Length(ip - C.c);
- if (dcmp(dis - C.r) > 0) return 0;
- Point dxy = vecunit(L.v) * sqrt(sqr(C.r) - sqr(dis));
- int ret = 0;
- sol[ret] = ip + dxy;
- if (OnSegment(sol[ret], L.p, L.point(1))) ret++;
- sol[ret] = ip - dxy;
- if (OnSegment(sol[ret], L.p, L.point(1))) ret++;
- return ret;
- }
- double SegCircleArea(Circle C, Point a, Point b) //线段切割圆
- {
- double a1 = angle(a - C.c);
- double a2 = angle(b - C.c);
- double da = fabs(a1 - a2);
- if (da > PI) da = PI * 2.0 - da;
- return dcmp(Cross(b - C.c, a - C.c)) * da * sqr(C.r) / 2.0;
- }
- double PolyCiclrArea(Circle C, Point *p, int n)//多边形与圆相交面积
- {
- double ret = 0.0;
- Point sol[2];
- p[n] = p[0];
- REP(i, n)
- {
- double t1, t2;
- int cnt = getSegCircleIntersection(Line(p[i], p[i+1]-p[i]), C, sol);
- if (cnt == 0)
- {
- if (!InCircle(p[i], C) || !InCircle(p[i+1], C)) ret += SegCircleArea(C, p[i], p[i+1]);
- else ret += Cross(p[i+1] - C.c, p[i] - C.c) / 2.0;
- }
- if (cnt == 1)
- {
- if (InCircle(p[i], C) && !InCircle(p[i+1], C)) ret += Cross(sol[0] - C.c, p[i] - C.c) / 2.0, ret += SegCircleArea(C, sol[0], p[i+1]);
- else ret += SegCircleArea(C, p[i], sol[0]), ret += Cross(p[i+1] - C.c, sol[0] - C.c) / 2.0;
- }
- if (cnt == 2)
- {
- if ((p[i] < p[i + 1]) ^ (sol[0] < sol[1])) swap(sol[0], sol[1]);
- ret += SegCircleArea(C, p[i], sol[0]);
- ret += Cross(sol[1] - C.c, sol[0] - C.c) / 2.0;
- ret += SegCircleArea(C, sol[1], p[i+1]);
- }
- }
- return fabs(ret);
- }
二维计算几何模板整理
最新推荐文章于 2020-07-11 15:08:04 发布