/*----------------------基础----------------------*/
const double pi = acos(-1);
const double eps = 1e-13;
const int maxp = 100 + 10; //多边形最多点数
int equal(double a, double b) //判断double相等
{
if (fabs(a - b) <= eps)
return 1;
else
return 0;
}
int sign(double num) //判断符号,负数返回-1,正数返回1,零返回0
{
if (fabs(num) < eps) return 0;
if (num > 0) return 1;
return -1;
}
int dcmp(double a, double b) //浮点数比较,小于-1,等于0,大于1
{
if (fabs(a - b) <= eps) return 0;
if (a < b) return -1;
return 1;
}
/*---------------------------------------------------*/
/*----------------------点,向量----------------------*/
struct point
{
double x, y;
point(double x1 = 0, double y1 = 0) { x = x1, y = y1; }
point operator+(point a) { return point(x + a.x, y + a.y); }
point operator-(point a) { return point(x - a.x, y - a.y); }
point operator*(double t) { return point(x * t, y * t); }
point operator/(double t) { return point(x / t, y / t); }
bool operator!=(point a) { return (!equal(x, a.x)) || (!(equal(y, a.y))); }
bool operator==(point a) { return equal(x, a.x) && equal(y, a.y); }
};
point rotate(point p, double theta) //将点p旋转θ弧度
{
return point(p.x * cos(theta) - p.y * sin(theta), p.x * sin(theta) + p.y * cos(theta));
}
double dot(point a, point b) { return a.x * b.x + a.y * b.y; } //点乘
double cross(point a, point b) { return a.x * b.y - a.y * b.x; } //叉乘
double len(point a) { return sqrt(dot(a, a)); } //向量模长
double dis(point a, point b) { return len(b - a); } //求两点距离
double dis2(point a, point b) { return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y); } //返回两点距离平方
double area(point a, point b, point c) { return fabs(cross(b - a, c - a) / 2.0); } //三点组成的三角形面积
double angle(point a, point b) //两向量夹角
{
double num = (dot(a, b) / len(a)) / len(b);
if (num > 1.0)
num = 1.0;
else if (num < -1.0)
num = -1.0;
return acos(num);
}
point projection(point p0, point p1, point p2) //求p0到p1p2上的投影
{
point v0 = p2 - p1;
point v1 = p0 - p1;
double num = dot(v0, v1) / len(v0);
point pro = p1 + v0 * (num / len(v0));
return pro;
}
double project(point p0, point p1, point p2)
{
point v0 = p2 - p1;
point v1 = p0 - p1;
return dot(v0, v1) / len(v0);
}
point reflection(point p0, point p1, point p2) //求p0关于p1p2上的对称点
{
point pro = projection(p0, p1, p2);
point ans = p0 + (pro - p0) * 2.0;
return ans;
}
point clp[maxn]; //平面最近点对
double merge1(int l, int r) //平面最近点对区间合并
{
double res = 0x3f3f3f3f;
if (r - l <= 3)
{
for (int i = l; i <= r; i++)
for (int j = i + 1; j <= r; j++)
res = min(res, dis(clp[i], clp[j]));
return res;
}
int mid = (l + r) >> 1;
res = min(res, merge1(l, mid));
res = min(res, merge1(mid + 1, r));
vector<point> tmp;
for (int i = l; i <= r; i++)
{
if (fabs(clp[i].x - clp[mid].x) >= res) continue;
tmp.push_back(clp[i]);
}
sort(tmp.begin(), tmp.end(), [](point a, point b) { return a.y < b.y; });
for (int i = 0; i < tmp.size(); i++)
{
for (int j = i + 1; j < tmp.size() && fabs(tmp[i].y - tmp[j].y) < res; j++)
{
res = min(res, dis(tmp[i], tmp[j]));
}
}
return res;
}
double cloest_pair(int cnt1) //求平面最近点对
{
sort(clp, clp + cnt1, [](point a, point b) {
if (equal(a.x, b.x)) return a.y < b.y;
return a.x < b.x;
});
double res = merge1(0, cnt1 - 1);
return res;
}
/*---------------------------------------------------*/
/*----------------------直线,线段----------------------*/
struct line //直线
{
point s, v;
line(point s1 = {0, 0}, point v1 = {1, 0}) { s = s1, v = v1; }
};
struct seg //线段
{
point u, v;
seg(point u1 = {0, 0}, point v1 = {1, 0}) { u = u1, v = v1; }
};
double angle(line a) //求直线与x轴夹角
{
return atan2(a.v.y, a.v.x);
}
point intersection(line a, line b) //求直线交点
{
point u = a.s - b.s;
double t = cross(b.v, u) / cross(a.v, b.v);
return a.s + a.v * t;
}
point intersection(seg a, seg b) //求线段交点,记得先判是否相交
{
line l1(a.u, a.v - a.u), l2(b.u, b.v - b.u);
return intersection(l1, l2);
}
int relation(line a, line b) //两直线关系,重合3,平行2,垂直1,否则0
{
if (sign(cross(a.v, b.v)) == 0)
{
point v = b.s - a.s;
if (sign(cross(v, a.v)) == 0)
return 3;
else
return 2;
}
else
{
if (dcmp(angle(a.v, b.v), pi / 2.0) == 0)
return 1;
else
return 0;
}
}
int is_intersection(seg a, seg b) //判断线段是否相交,相交返回1,否则返回0
{
if (min(b.u.x, b.v.x) > max(a.u.x, a.v.x) || min(a.u.x, a.v.x) > max(b.u.x, b.v.x)) return 0;
if (min(b.u.y, b.v.y) > max(a.u.y, a.v.y) || min(a.u.y, a.v.y) > max(b.u.y, b.v.y)) return 0;
double num1, num2, num3, num4;
num1 = cross(b.u - a.u, b.v - a.u);
num2 = cross(b.u - a.v, b.v - a.v);
num3 = cross(a.u - b.u, a.v - b.u);
num4 = cross(a.u - b.v, a.v - b.v);
if (sign(num1 * num2) <= 0 && sign(num3 * num4) <= 0) return 1;
return 0;
}
int on_seg(point a, seg b) //判断点是否在线段上
{
point v1, v2;
v1 = b.u - a;
v2 = b.v - a;
if (a == b.u || a == b.v) return 1;
if (equal(dot(v1, v2) / len(v1) / len(v2), -1.0)) return 1;
return 0;
}
double dis(point a, seg b) //求点到线段距离
{
if (on_seg(a, b)) return 0.0;
point v1, v2;
v1 = b.v - b.u;
v2 = a - b.u;
double num = dot(v1, v2) / (len(v1) * len(v1));
if (dcmp(num, 0.0) <= 0)
return dis(a, b.u);
else if (dcmp(num, 1.0) >= 0)
return dis(a, b.v);
return dis(a, projection(a, b.u, b.v));
}
double dis(seg a, seg b) //求两线段距离
{
if (is_intersection(a, b)) return 0.0;
double num1, num2;
num1 = min(dis(a.u, b), dis(a.v, b));
num2 = min(dis(b.u, a), dis(b.v, a));
return min(num1, num2);
}
/*---------------------------------------------------*/
/*----------------------多边形----------------------*/
struct polygon
{
point p[maxp];
int cnt;
polygon() { cnt = 0; }
void add(point pp) { p[cnt++] = pp; } //加点
};
int stk[maxn]; //求凸包栈
bool used[maxn]; //求凸包标记已使用点
polygon an_poly, ans_an_poly;
void andrew() //求凸包
{
int top = 0;
sort(an_poly.p, an_poly.p + an_poly.cnt, [](point a, point b) {
if (equal(a.x, b.x)) return a.y < b.y;
return a.x < b.x;
});
for (int i = 0; i < an_poly.cnt; i++)
{
while (top >= 2 && cross(an_poly.p[stk[top - 1]] - an_poly.p[stk[top - 2]], an_poly.p[i] - an_poly.p[stk[top - 1]]) <= 0)
{
if (cross(an_poly.p[stk[top - 1]] - an_poly.p[stk[top - 2]], an_poly.p[i] - an_poly.p[stk[top - 1]]) < 0)
used[stk[--top]] = false;
else
top--;
}
stk[top++] = i;
used[i] = true;
}
int tmp = top;
used[0] = false;
for (int i = an_poly.cnt - 1; i >= 0; i--)
{
if (used[i]) continue;
while (top >= tmp && cross(an_poly.p[stk[top - 1]] - an_poly.p[stk[top - 2]], an_poly.p[i] - an_poly.p[stk[top - 1]]) <= 0)
top--;
stk[top++] = i;
used[i] = true;
}
for (int i = 0; i < top; i++) ans_an_poly.add(an_poly.p[stk[i]]);
}
int is_convex(polygon poly) //判断是否是凸包(逆时针)
{
if (poly.cnt < 3) return 0;
point p[maxp];
int tot = 0;
p[tot++] = poly.p[0];
p[tot++] = poly.p[1];
for (int i = 2; i < poly.cnt; i++)
{
point v1, v2;
v1 = p[tot - 1] - p[tot - 2];
v2 = poly.p[i] - p[tot - 1];
if (sign(cross(v1, v2)) == -1) return 0;
p[tot++] = poly.p[i];
}
point v1, v2;
v1 = p[tot - 1] - p[tot - 2];
v2 = poly.p[0] - p[tot - 1];
if (sign(cross(v1, v2)) == -1) return 0;
return 1;
}
double area(polygon &poly) //求多边形面积
{
double res = 0;
if (poly.cnt < 3) return 0;
for (int i = 1; i < poly.cnt; i++)
{
res += cross(poly.p[i], poly.p[i - 1]);
}
res += cross(poly.p[0], poly.p[poly.cnt - 1]);
return fabs(res / 2.0);
}
double len(polygon &poly) //求多边形周长
{
double res = 0;
for (int i = 1; i < poly.cnt; i++)
{
res += len(poly.p[i] - poly.p[i - 1]);
}
res += len(poly.p[0] - poly.p[poly.cnt - 1]);
return res;
}
int point_in_poly(point p, polygon poly) //点与多边形关系,2多边形内 ,1线段上,0多边形外
{
seg seg1[110]; //多边形边的数量
int tot = 0;
for (int i = 1; i < poly.cnt; i++)
{
point u, v;
u = poly.p[i - 1];
v = poly.p[i];
seg1[tot++] = seg(u, v);
}
seg1[tot++] = seg(poly.p[poly.cnt - 1], poly.p[0]);
int flag = 0;
point u, v;
for (int i = 0; i < tot; i++)
{
if (on_seg(p, seg1[i])) return 1;
u = seg1[i].u;
v = seg1[i].v;
if (((sign(u.y - p.y) > 0) != (sign(v.y - p.y) > 0)) && sign(p.x - (p.y - u.y) * (u.x - v.x) / (u.y - v.y) - u.x) < 0)
flag ^= 2;
}
return flag;
}
int on_right(line a, line b, line c) //判断bc交点是否在a的右边
{
point o = intersection(b, c);
if (sign(cross(o - a.s, a.v)) == 1) return 1;
return 0;
}
line hpi_li[maxn];
int hpi_q[maxn]; //半平面交双端队列
polygon half_plane_intersection(int cnt) //求半平面交
{
sort(hpi_li, hpi_li + cnt, [&](line a, line b) {
double num1, num2;
num1 = atan2(a.v.y, a.v.x);
num2 = atan2(b.v.y, b.v.x);
if (equal(num1, num2)) return sign(cross(b.s - a.s, a.v)) == 1;
return num1 < num2;
});
int hh = 0, tt = -1;
for (int i = 0; i < cnt; i++)
{
if (i > 0 && dcmp(angle(hpi_li[i]), angle(hpi_li[i - 1])) == 0) continue;
while (hh + 1 <= tt && on_right(hpi_li[i], hpi_li[hpi_q[tt - 1]], hpi_li[hpi_q[tt]])) tt--;
while (hh + 1 <= tt && on_right(hpi_li[i], hpi_li[hpi_q[hh]], hpi_li[hpi_q[hh + 1]])) hh++;
hpi_q[++tt] = i;
}
while (hh + 1 <= tt && on_right(hpi_li[hpi_q[hh]], hpi_li[hpi_q[tt - 1]], hpi_li[hpi_q[tt]])) tt--;
while (hh + 1 <= tt && on_right(hpi_li[hpi_q[tt]], hpi_li[hpi_q[hh]], hpi_li[hpi_q[hh + 1]])) hh++;
hpi_q[++tt] = hpi_q[hh];
int k = 0;
polygon ans;
for (int i = hh; i < tt; i++)
ans.add(intersection(hpi_li[hpi_q[i]], hpi_li[hpi_q[i + 1]]));
return ans;
}
polygon rc_poly;
double rotating_calipers() //旋转卡壳求凸包最远点对距离平方，最后一个点和第一个点相同
{
double res = 0;
int top1 = rc_poly.cnt - 1;
for (int i = 0, j = 2; i < top1; i++)
{
point s = rc_poly.p[i], e = rc_poly.p[i + 1];
while (area(s, e, rc_poly.p[j]) <= area(s, e, rc_poly.p[j + 1])) j = (j + 1) % top1;
res = max(res, max(dis2(rc_poly.p[j], s), dis2(rc_poly.p[j], e)));
}
return res;
}
polygon mrc_poly, ans_mrc_poly; //最小矩形覆盖用
void mini_rectangle_coverage() //最小矩形覆盖
{
double minarea = 1e20;
for (int i = 0, j = 2, k = 2, l = 2; i < mrc_poly.cnt; i++)
{
point s, e;
s = mrc_poly.p[i], e = mrc_poly.p[i + 1];
double num1, num2;
num1 = area(s, e, mrc_poly.p[j]);
num2 = area(s, e, mrc_poly.p[j + 1]);
while (dcmp(area(s, e, mrc_poly.p[j]), area(s, e, mrc_poly.p[j + 1])) < 0) j = (j + 1) % mrc_poly.cnt;
while (dcmp(project(mrc_poly.p[k], s, e), project(mrc_poly.p[k + 1], s, e)) < 0) k = (k + 1) % mrc_poly.cnt;
if (i == 0) l = j;
while (dcmp(project(mrc_poly.p[l], s, e), project(mrc_poly.p[l + 1], s, e)) > 0) l = (l + 1) % mrc_poly.cnt;
point right = mrc_poly.p[k], up = mrc_poly.p[j], left = mrc_poly.p[l];
double h = 2.0 * area(s, e, up) / len(e - s);
double w = fabs(dot(right - left, e - s)) / len(e - s);
if (h * w < minarea)
{
minarea = h * w;
ans_mrc_poly.cnt = 4;
ans_mrc_poly.p[0] = s + ((e - s) / len(e - s)) * project(right, s, e);
ans_mrc_poly.p[3] = s + ((e - s) / len(e - s)) * project(left, s, e);
point ahh = rotate(e - s, pi / 2.0) / len(e - s);
ans_mrc_poly.p[1] = ans_mrc_poly.p[0] + ahh * h;
ans_mrc_poly.p[2] = ans_mrc_poly.p[3] + ahh * h;
}
}
}
/*---------------------------------------------------*/
/*----------------------圆----------------------*/
struct circle
{
point c;
double r;
circle(point c1 = {0, 0}, double r1 = 0) { c = c1, r = r1; }
};
double area(circle ci) { return pi * ci.r * ci.r; } //圆面积
double len(circle ci) { return 2.0 * pi * ci.r; } //圆周长
circle get_circle(point a, point b, point c) //求三点外接圆
{
circle res;
line l1((a + b) / 2.0, rotate(b - a, pi / 2.0));
line l2((a + c) / 2.0, rotate(c - a, pi / 2.0));
res.c = intersection(l1, l2);
res.r = dis(res.c, a);
return res;
}
circle mini_circle_coverage(point p[], int cnt)
{
srand(time(0));
random_shuffle(p, p + cnt);
circle res(p[0], 0);
for (int i = 1; i < cnt; i++)
{
if (dcmp(dis(res.c, p[i]), res.r) <= 0) continue;
circle now(p[i], 0);
for (int j = 0; j < i; j++)
{
if (dcmp(dis(now.c, p[j]), now.r) <= 0) continue;
now.c = (p[i] + p[j]) / 2.0;
now.r = dis(now.c, p[i]);
for (int k = 0; k < j; k++)
{
if (dcmp(dis(now.c, p[k]), now.r) <= 0) continue;
now = get_circle(p[i], p[j], p[k]);
}
}
res = now;
}
return res;
}
/*---------------------------------------------------*/
/*----------------------三维----------------------*/
struct point3 //三维点
{
double x, y, z;
point3(double x1 = 0, double y1 = 0, double z1 = 0) { x = x1, y = y1, z = z1; }
point3 operator+(point3 a) { return point3(x + a.x, y + a.y, z + a.z); }
point3 operator-(point3 a) { return point3(x - a.x, y - a.y, z - a.z); }
point3 operator*(double t) { return point3(x * t, y * t, z * t); }
point3 operator/(double t) { return point3(x / t, y / t, z / t); }
bool operator!=(point3 a) { return (!equal(x, a.x)) || (!(equal(y, a.y))) || (!equal(z, a.z)); }
bool operator==(point3 a) { return equal(x, a.x) && equal(y, a.y) && equal(z, a.z); }
} pt3[maxn];
struct line3
{
point3 s, e;
line3() {}
line3(point3 _s, point3 _e)
{
s = _s;
e = _e;
}
};
struct plane
{
int p[3];
plane(int a = 0, int b = 0, int c = 0) { p[0] = a, p[1] = b, p[2] = c; }
};
double dis(point3 a, point3 b) //两点距离
{
double dx = a.x - b.x, dy = a.y - b.y, dz = a.z - b.z;
return sqrt(dx * dx + dy * dy + dz * dz);
}
double rand_eps() { return ((double)rand() / RAND_MAX - 0.5) * eps; } //随机抖动
point3 shake(point3 a) { return a + point3(rand_eps(), rand_eps(), rand_eps()); } //抖动
double len(point3 a) { return sqrt(a.x * a.x + a.y * a.y + a.z * a.z); } //向量长度
point3 cross(point3 a, point3 b) { return point3(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 dot(point3 a, point3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; } //点积
point3 norm(plane a) { return cross(pt3[a.p[1]] - pt3[a.p[0]], pt3[a.p[2]] - pt3[a.p[0]]); } //法向量
double area(plane a) { return fabs(len(cross(pt3[a.p[1]] - pt3[a.p[0]], pt3[a.p[2]] - pt3[a.p[0]]))) / 2.0; } //面积
int above(point3 p, plane pl) //点在面上和面内1，点在面下0
{
point3 nor = norm(pl);
double pro = dot(p - pt3[pl.p[0]], nor);
return pro >= 0;
}
plane pl[maxn * 10], np[maxn * 10]; //三维凸包结果,当前凸包
int gg[1][1]; //边是否被照射
int convex3(int nn) //三维凸包，返回凸包面数
{ //记得加抖动
int tot = 0;
pl[tot++] = {0, 1, 2};
pl[tot++] = {2, 1, 0};
for (int i = 3; i < nn; i++)
{
int cnt = 0;
for (int j = 0; j < tot; j++)
{
int x = above(pt3[i], pl[j]);
if (x == 0) np[cnt++] = pl[j];
for (int k = 0; k < 3; k++)
{
if (x == 0)
gg[pl[j].p[k]][pl[j].p[(k + 1) % 3]] = 0;
else
gg[pl[j].p[k]][pl[j].p[(k + 1) % 3]] = 1;
}
}
for (int j = 0; j < tot; j++)
{
for (int k = 0; k < 3; k++)
{
int u, v;
u = pl[j].p[k];
v = pl[j].p[(k + 1) % 3];
if (gg[u][v] == 1 && gg[v][u] == 0) np[cnt++] = plane(u, v, i);
}
}
tot = cnt;
for (int j = 0; j < tot; j++) pl[j] = np[j];
}
return tot;
}
point3 trunc(point3 a, double r)
{
double l = len(a);
if (!sign(l)) return a;
r /= l;
return point3(a.x * r, a.y * r, a.z * r);
}
point3 rotate(point3 p, line3 vec, double ang) //点p绕vec旋转ang弧度
{
if (sign(len(cross((vec.s - p), (vec.e - p)))) == 0) return p;
point3 f1 = cross((vec.e - vec.s), (p - vec.s));
point3 f2 = cross((vec.e - vec.s), (f1));
double len1 = len(cross((vec.s - p), (vec.e - p))) / dis(vec.s, vec.e);
f1 = trunc(f1, len1);
f2 = trunc(f2, len1);
point3 h = p + f2;
point3 pp = h + f1;
return h + ((p - h) * cos(ang)) + ((pp - h) * sin(ang));
}
/*---------------------------------------------------*/
浙公网安备 33010602011771号