struct node{
double x,y;
};
node a,b,c;
//求两个点之间的长度
double len(node a,node b) {
double tmp = sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
return tmp;
}
//给出三个点，求三角形的面积 海伦公式: p=(a+b+c)/2，S=sqrt(p(p-a)(p-b)(p-c))
double Area(node a,node b,node c){
double lena = len(a,b);
double lenb = len(b,c);
double lenc = len(a,c);
double p = (lena + lenb + lenc) / 2.0;
double S = sqrt(p * (p - lena) * (p - lenb) * (p - lenc));
return S;
}
//三角形求每条边对应的圆心角
void Ran() {
double lena = len(a,b);
double lenb = len(b,c);
double lenc = len(a,c);
double A = acos((lenb * lenb + lenc * lenc - lena * lena) / (2 * lenb * lenc));
double B = acos((lena * lena + lenc * lenc - lenb * lenb) / (2 * lena * lenc));
double C = acos((lena * lena + lenb * lenb - lenc * lenc) / (2 * lena * lenb));
}
//求外接圆半径r = a * b * c / 4S
double R(node a,node b,node c) {
double lena = len(a,b);
double lenb = len(b,c);
double lenc = len(a,c);
double S = Area(a,b,c);
double R = lena * lenb * lenc / (4.0 * S);
}
const double eps = 1e-8, pi = acos(-1);
const int N = 1e5 + 10;
struct point {
double x, y;
point() {}
point(double _x, double _y) {
x = _x;
y = _y;
}
point operator-(const point &p) const {
return point(x - p.x, y - p.y);
}
point operator+(const point &p) const {
return point(x + p.x, y + p.y);
}
double operator*(const point &p) const {
return x * p.y - y * p.x;
}
double operator/(const point &p) const {
return x * p.x + y * p.y;
}
}a[N], b[N << 1];
struct line {
point p1, p2;
double ang;
line() {}
line(point p01, point p02) {
p1 = p01;
p2 = p02;
}
void getang() {
ang = atan2(p2.y - p1.y, p2.x - p1.x);
}
}c[N], q[N];
struct circle
{
double x, y, r;
}cc;
struct point3 {
double x, y, z;
point3() {}
point3(double _x, double _y, double _z) {
x = _x;
y = _y;
z = _z;
}
point3 operator-(const point3 &p) const {
return point3(x - p.x, y - p.y, z - p.z);
}
point3 operator+(const point3 &p) const {
return point3(x + p.x, y + p.y, z - p.z);
}
point3 operator*(const point3 &p) const {
return point3(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x);
}
double operator/(const point3 &p) const {
return x * p.x + y * p.y + z * p.z;
}
};
int n;
bool cmp(const point &aa, const point &bb) {
return aa.x < bb.x || (aa.x == bb.x && aa.y < bb.y);
}
bool cmpl(const line &aa, const line &bb) {
if (abs(aa.ang - bb.ang) > eps)
return aa.ang < bb.ang;
return (bb.p2 - aa.p1) * (bb.p1 - bb.p2) > 0;
}
int sgn(double x) {
if (x > eps)
return 1;
if (x < -eps)
return -1;
return 0;
}
double getdist2(const point &aa) {
return (aa.x * aa.x + aa.y * aa.y);
}
bool online(const point &aa, const point &bb, const point &cc) { //三点共线判断
if (!(min(bb.x, cc.x) <= aa.x && aa.x <= max(bb.x, cc.x)))
return 0;
if (!(min(bb.y, cc.y) <= aa.y && aa.y <= max(bb.y, cc.y)))
return 0;
return !sgn((bb - aa) * (cc - aa));
}
bool checkline(const point &aa, const point &bb, const point &cc, const point &dd) { //线段是否相交
int c1 = sgn((bb - aa) * (cc - aa)), c2 = sgn((bb - aa) * (dd - aa)),
c3 = sgn((dd - cc) * (aa - cc)), c4 = sgn((dd - cc) * (bb - cc));
if (c1 * c2 < 0 && c3 * c4 < 0)
return 1;//规范相交
if (c1 == 0 && c2 == 0) {
if (min(aa.x, bb.x) > max(cc.x, dd.x) || min(cc.x, dd.x) > max(aa.x, bb.x)
|| min(aa.y, bb.y) > max(cc.y, dd.y) || min(cc.y, dd.y) > max(aa.y, bb.y))
return 0;
return 1;//共线有重叠
}
if (online(cc, aa, bb) || online(dd, aa, bb)
|| online(aa, cc, dd) || online(bb, cc, dd))
return 1;//端点相交
return 0;
}
point getpoint(const point &aa, const point &bb,
const point &cc, const point &dd) { //两直线交点
double a1, b1, c1, a2, b2, c2;
point re;
a1 = aa.y - bb.y;
b1 = bb.x - aa.x;
c1 = aa * bb;
a2 = cc.y - dd.y;
b2 = dd.x - cc.x;
c2 = cc * dd;
//以下为交点横纵坐标
re.x = (c1 * b2 - c2 * b1) / (a2 * b1 - a1 * b2);
re.y = (a2 * c1 - a1 * c2) / (a1 * b2 - a2 * b1);
return re;
}
bool protrusion(point aa[]){ //判断多边形凹(凸)
int flag = sgn((aa[1] - aa[0]) * (aa[2] - aa[0]));
for (int i = 1; i < n; ++i)
if (sgn((aa[(i + 1) % n] - aa[i]) * (aa[(i + 2) % n] - aa[i])) != flag)
return 1;
return 0;
}
bool inpolygon(const point &aa) { //点在凸多边形内(外)
int flag = sgn((a[0] - aa) * (a[1] - aa));
for (int i = 1; i < n; ++i)
if (sgn((a[i] - aa) * (a[(i + 1) % n] - aa)) != flag)
return 0;
return 1;
}
void transline(line &aa, double dist) { //逆时针方向平移线段
double d = sqrt((aa.p1.x - aa.p2.x) * (aa.p1.x - aa.p2.x) +
(aa.p1.y - aa.p2.y) * (aa.p1.y - aa.p2.y));
point ta;
ta.x = aa.p1.x + dist / d * (aa.p1.y - aa.p2.y);
ta.y = aa.p1.y - dist / d * (aa.p1.x - aa.p2.x);
aa.p2 = ta + aa.p2 - aa.p1;
aa.p1 = ta;
}
void convexhull(point aa[], point bb[]) { //凸包
int len = 0;
sort(aa, aa + n, cmp);
bb[len++] = aa[0];//bb数组要开2倍(防止出现直线)
bb[len++] = aa[1];
for (int i = 2; i < n; ++i) {
while (len > 1 && (aa[i] - bb[len - 2]) * (bb[len - 1] - bb[len - 2]) > 0)
//若严格则加上等于
--len;
bb[len++] = aa[i];
}
int t = len;
for (int i = n - 2; i >= 0; --i) {
while (len > t && (aa[i] - bb[len - 2]) * (bb[len - 1] - bb[len - 2]) > 0)
//同上
--len;
bb[len++] = aa[i];
}
--len;
n = len;
}
bool checkout(const line &aa, const line &bb, const line &cc) { //检查交点是否在向量顺时针侧
point p = getpoint(aa.p1, aa.p2, bb.p1, bb.p2);
return (cc.p1 - p) * (cc.p2 - p) < 0;//如果不允许共线或算面积 则此处不取等
}
double halfplane(point aa[], line bb[]) { //半平面交
sort(bb, bb + n, cmpl);
int n2 = 1;
for (int i = 1; i < n; ++i) {
if (bb[i].ang - bb[i - 1].ang > eps)
++n2;
bb[n2 - 1] = bb[i];
}
n = n2;
int front = 0, tail = 0;
q[tail++] = bb[0], q[tail++] = bb[1];
for (int i = 2; i < n; ++i) {
while (front + 1 < tail && checkout(q[tail - 2], q[tail - 1], bb[i]))
--tail;
while (front + 1 < tail && checkout(q[front], q[front + 1], bb[i]))
++front;
q[tail++] = bb[i];
}
while (front + 2 < tail && checkout(q[tail - 2], q[tail - 1], q[front]))
--tail;
while (front + 2 < tail && checkout(q[front], q[front + 1], q[tail - 1]))
++front;
if (front + 2 >= tail)
return 0;
int j = 0;
for (int i = front; i < tail; ++i, ++j) {
aa[j] = getpoint(q[i].p1, q[i].p2, q[(i != tail - 1 ? i + 1 : front)].p1,
q[(i != tail - 1 ? i + 1 : front)].p2);
}
double re = 0;
for (int i = 1; i < j - 1; ++i)
re += (aa[i] - aa[0]) * (aa[i + 1] - aa[0]);
return abs(re * 0.5);
}
double calipers(point aa[]) { //旋转卡壳
double re = 0;
aa[n] = aa[0];
int now = 1;
for (int i = 0; i < n; ++i) {
while ((aa[i + 1] - aa[i]) * (aa[now + 1] - aa[i]) >
(aa[i + 1] - aa[i]) * (aa[now] - aa[i]))
now = (now == n - 1 ? 0 : now + 1);
re = max(re, getdist2(aa[now] - aa[i]));
}
return re;
}
line bisector(const point &aa, const point &bb) { //中垂线(正方形)
double mx, my;
mx = (aa.x + bb.x) / 2;
my = (aa.y + bb.y) / 2;
line cc;
cc.p1.x = mx - (aa.y - bb.y) / 2;
cc.p1.y = my + (aa.x - bb.x) / 2;
cc.p2 = aa + bb - cc.p1;
cc.getang();
return cc;
}
point rotate(const point &p, double cost, double sint) { //逆时针向量旋转
double x = p.x, y = p.y;
return point(x * cost - y * sint, x * sint + y * cost);
}
void getpoint(circle c1, circle c2) { //已确保两圆有交点时求出两圆交点
long double dab = sqrt(getdist2(point(c1.x, c1.y) -
point(c2.x, c2.y)));
if (c1.r > c2.r)
swap(c1, c2);
long double cost = (c1.r * c1.r + dab * dab - c2.r * c2.r) /
(c1.r * dab * 2);
long double sint = sqrt(1 - cost * cost);
point re = rotate(point(c2.x, c2.y) - point(c1.x, c1.y), cost, sint);
re.x = c1.x + re.x * (c1.r / dab);
re.y = c1.y + re.y * (c1.r / dab);
point re2 = rotate(point(c2.x, c2.y) - point(c1.x, c1.y), cost, -sint);
re2.x = c1.x + re2.x * (c1.r / dab);
re2.y = c1.y + re2.y * (c1.r / dab);
}
double mycos(double B, double C, double A) { //余弦定理 给定边长
return (B * B + C * C - A * A) / (B * C * 2);
}
double mycos2(const point &aa, const point &bb, const point &cc) {//余弦定理 给定点坐标
double C2 = getdist2(aa - bb), A2 = getdist2(bb - cc), B2 = getdist2(cc - aa);
return (B2 + C2 - A2) / (sqrt(B2 * C2) * 2);
}
point mirror_point(const point &aa, const point &bb, const point &cc) { //轴对称点 也可用来求垂足
double cost, sint;
cost = mycos2(bb, cc, aa);
sint = sqrt(1 - cost * cost);
if ((cc - bb) * (aa - bb) > 0)
sint = -sint;
point re;
re = rotate(aa - bb, cost, sint) + bb;
re = rotate(re - bb, cost, sint) + bb;
return re;
}
double area(const point &aa, const point &bb, const point &cc) {//求三角形面积
return abs((bb - aa) * (cc - aa) / 2);
}
void max_triangle(const point aa[]) { //求凸包上最大三角形
double ans = 0;
int p1, p2, p3;
for (int i = 0; i < n; ++i) {
int j = (i + 1) % n;
int k = (j + 1) % n;
while (k != i && area(aa[i], aa[j], aa[k]) < area(aa[i], aa[j], aa[(k + 1) % n]))
k = (k + 1) % n;
if (k == i)
continue;
int kk = (k + 1) % n;
while (j != kk && k != i) {
if (ans < area(aa[i], aa[j], aa[k])) {
ans = area(aa[i], aa[j], aa[k]);//三角形面积
p1 = i;
p2 = j;
p3 = k;
}
while (k != i && area(aa[i], aa[j], aa[k]) < area(aa[i], aa[j], aa[(k + 1) % n]))
k = (k + 1) % n;
j = (j + 1) % n;
}
}
point q1, q2, q3;
q1 = b[p2] + b[p3] - b[p1];
q2 = b[p1] + b[p3] - b[p2];
q3 = b[p1] + b[p2] - b[p3];
//三角形上三个点
}
void get_panel(const point3 &p1, const point3 &p2, const point3 &p3,
double &A, double &B, double &C, double &D) { //由三点确定一个平面方程
A = ((p2.y - p1.y) * (p3.z - p1.z) - (p2.z - p1.z) * (p3.y - p1.y));
B = ((p2.z - p1.z) * (p3.x - p1.x) - (p2.x - p1.x) * (p3.z - p1.z));
C = ((p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x));
D = (-(A * p1.x + B * p1.y + C * p1.z));
}
double dist_panel(const point3 &pt, double A, double B, double C, double D) { //点到平面距离
return abs(A * pt.x + B * pt.y + C * pt.z + D) /
sqrt(A * A + B * B + C * C);
}
