mirror of
https://github.com/Kiritow/OJ-Problems-Source.git
synced 2024-03-22 13:11:29 +08:00
425 lines
13 KiB
C++
425 lines
13 KiB
C++
#include <cstdio>
|
|
#include <cmath>
|
|
#include <vector>
|
|
#include <deque>
|
|
#include <algorithm>
|
|
const double eps = 1e-13;
|
|
const double pi = 3.14159265358979324;
|
|
struct point2 {
|
|
double x, y;
|
|
point2& operator += (point2 a) { x+=a.x, y+=a.y; return *this; }
|
|
point2& operator -= (point2 a) { x-=a.x, y-=a.y; return *this; }
|
|
point2& operator *= (double a) { x*=a, y*=a; return *this; }
|
|
point2& operator /= (double a) { x/=a, y/=a; return *this; }
|
|
};
|
|
point2 operator + (point2 a, point2 b) { point2 c(a); c += b; return c; }
|
|
point2 operator - (point2 a, point2 b) { point2 c(a); c -= b; return c; }
|
|
point2 operator * (point2 a, double b) { point2 c(a); c *= b; return c; }
|
|
point2 operator * (double a, point2 b) { point2 c(b); c *= a; return c; }
|
|
point2 operator / (point2 a, double b) { point2 c(a); c /= b; return c; }
|
|
double operator * (point2 a, point2 b) { return a.x*b.x+a.y*b.y; }
|
|
double operator % (point2 a, point2 b) { return a.x*b.y-a.y*b.x; }
|
|
double dis(point2 a) { return sqrt(a.x * a.x + a.y * a.y); }
|
|
double arg(point2 a) { return atan2(a.y, a.x); }
|
|
point2 rotate(point2 a, double th) {
|
|
point2 b;
|
|
b.x = a.x * cos(th) - a.y * sin(th);
|
|
b.y = a.x * sin(th) + a.y * cos(th);
|
|
return b;
|
|
}
|
|
int parallel(point2 a, point2 b) {
|
|
return a * a < eps * eps || b * b < eps * eps
|
|
|| (a % b) * (a % b) / ((a * a) * (b * b)) < eps * eps;
|
|
}
|
|
int perpend(point2 a, point2 b) {
|
|
return a * a < eps * eps || b * b < eps * eps
|
|
|| (a * b) * (a * b) / ((a * a) * (b * b)) < eps * eps;
|
|
}
|
|
struct line2 {
|
|
point2 a, s;
|
|
};
|
|
struct circle2 {
|
|
point2 a;
|
|
double r;
|
|
};
|
|
double point_line_dis(point2 a, line2 b, point2 *res = NULL) {
|
|
point2 p;
|
|
p = b.a + ((a - b.a) * b.s) / (b.s * b.s) * b.s;
|
|
if (res != NULL) *res = p;
|
|
return dis(a - p);
|
|
}
|
|
int line_line_cross(line2 a, line2 b, point2 *res = NULL) {
|
|
if (parallel(a.s, b.s))
|
|
if (parallel(b.a - a.a, a.s))
|
|
return -1;
|
|
else
|
|
return 0;
|
|
double k1 = (b.a - a.a) % b.s / (a.s % b.s);
|
|
if (res != NULL) *res = a.a + k1 * a.s;
|
|
return 1;
|
|
}
|
|
int segment_segment_cross(line2 a, line2 b, point2 *res = NULL) {
|
|
if (a.s * a.s < eps * eps && b.s * b.s < eps * eps)
|
|
if ((b.a - a.a) * (b.a - a.a) < eps * eps) {
|
|
if (res != NULL) *res = a.a;
|
|
return 1;
|
|
} else
|
|
return 0;
|
|
if (parallel(a.s, b.s) && parallel(b.a - a.a, a.s) && parallel(a.a - b.a, b.s)) {
|
|
double y1, y2, y3, y4;
|
|
point2 y1p = a.a, y2p = a.a + a.s, y3p = b.a, y4p = b.a + b.s;
|
|
if (std::abs(a.s.x) < std::abs(a.s.y) || std::abs(b.s.x) < std::abs(b.s.y))
|
|
y1 = y1p.y, y2 = y2p.y, y3 = y3p.y, y4 = y4p.y;
|
|
else
|
|
y1 = y1p.x, y2 = y2p.x, y3 = y3p.x, y4 = y4p.x;
|
|
if (y1 > y2) std::swap(y1, y2), std::swap(y1p, y2p);
|
|
if (y3 > y4) std::swap(y3, y4), std::swap(y3p, y4p);
|
|
if (y2 - y1 < y4 - y3) std::swap(y1, y3), std::swap(y1p, y3p), std::swap(y2, y4), std::swap(y2p, y4p);
|
|
if (y3 > y2 + (y2 - y1) * eps || y4 < y1 - (y2 - y1) * eps)
|
|
return 0;
|
|
else if (fabs(y3 - y2) < (y2 - y1) * eps || fabs(y3 - y4) < eps) {
|
|
if (res != NULL) *res = y3p;
|
|
return 1;
|
|
} else if (fabs(y4 - y1) < (y2 - y1) * eps || fabs(y1 - y2) < eps) {
|
|
if (res != NULL) *res = y1p;
|
|
return 1;
|
|
} else
|
|
return -1;
|
|
} else {
|
|
double k1 = (b.a - a.a) % a.s, k2 = (b.a + b.s - a.a) % a.s;
|
|
k1 /= a.s * a.s, k2 /= a.s * a.s;
|
|
double k3 = (a.a - b.a) % b.s, k4 = (a.a + a.s - b.a) % b.s;
|
|
k3 /= b.s * b.s, k4 /= b.s * b.s;
|
|
int ret = (k1 < eps && k2 > -eps || k1 > -eps && k2 < eps)
|
|
&& (k3 < eps && k4 > -eps || k3 > -eps && k4 < eps);
|
|
if (ret) line_line_cross(a, b, res);
|
|
return ret;
|
|
}
|
|
}
|
|
int line_circle_cross(line2 a, circle2 b, point2 *res1 = NULL, point2 *res2 = NULL) {
|
|
point2 p;
|
|
double d = point_line_dis(b.a, a, &p);
|
|
if (d / b.r > 1 + eps)
|
|
return 0;
|
|
else if (d / b.r > 1 - eps) {
|
|
if (res1 != NULL) *res1 = p;
|
|
return 1;
|
|
} else {
|
|
d = sqrt(b.r * b.r - d * d) / dis(a.s);
|
|
if (res1 != NULL) *res1 = p + d * a.s;
|
|
if (res2 != NULL) *res2 = p - d * a.s;
|
|
return 2;
|
|
}
|
|
}
|
|
int circle_circle_cross(circle2 a, circle2 b, point2 *res1 = NULL, point2 *res2 = NULL) {
|
|
double d = dis(a.a - b.a);
|
|
point2 u = (b.a - a.a) / d;
|
|
if (d / (a.r + b.r) > 1 + eps)
|
|
return 0;
|
|
else if (d / (a.r + b.r) > 1 - eps) {
|
|
if (res1 != NULL) *res1 = a.a + u * a.r;
|
|
return 1;
|
|
} else if ((d - fabs(a.r - b.r)) / (a.r + b.r) > eps) {
|
|
double th = acos((a.r * a.r + d * d - b.r * b.r) / (2 * a.r * d));
|
|
if (res1 != NULL) *res1 = a.a + rotate(u * a.r, th);
|
|
if (res2 != NULL) *res2 = a.a + rotate(u * a.r, -th);
|
|
return 2;
|
|
} else if ((d - fabs(a.r - b.r)) / (a.r + b.r) > -eps) {
|
|
if (a.r / b.r < 1 - eps) {
|
|
if (res1 != NULL) *res1 = b.a - u * b.r;
|
|
return 1;
|
|
} else if (a.r / b.r > 1 + eps) {
|
|
if (res1 != NULL) *res1 = a.a + u * a.r;
|
|
return 1;
|
|
} else return -1;
|
|
} else
|
|
return 0;
|
|
}
|
|
int point_circle_tangent(point2 a, circle2 b, point2 *res1 = NULL, point2 *res2 = NULL) {
|
|
double d = dis(a - b.a);
|
|
point2 u = (a - b.a) / d;
|
|
if (d / b.r > 1 + eps) {
|
|
double th = acos(b.r / d);
|
|
if (res1 != NULL) *res1 = b.a + rotate(u * b.r, th);
|
|
if (res2 != NULL) *res2 = b.a + rotate(u * b.r, -th);
|
|
return 2;
|
|
} else if (d / b.r > 1 - eps) {
|
|
if (res1 != NULL) *res1 = a;
|
|
return 1;
|
|
} else
|
|
return 0;
|
|
}
|
|
int circle_circle_tangent(circle2 a, circle2 b, line2 *reso1 = NULL, line2 *reso2 = NULL, line2 *resi1 = NULL, line2 *resi2 = NULL) {
|
|
double d = dis(a.a - b.a);
|
|
point2 u = (b.a - a.a) / d;
|
|
int cnt = 0;
|
|
if ((d - fabs(a.r - b.r)) / (a.r + b.r) > eps) {
|
|
double th = acos((a.r - b.r) / d);
|
|
if (reso1 != NULL) {
|
|
reso1->a = a.a + rotate(u * a.r, th);
|
|
reso1->s = b.a + rotate(u * b.r, th) - reso1->a;
|
|
}
|
|
if (reso2 != NULL) {
|
|
reso2->a = a.a + rotate(u * a.r, -th);
|
|
reso2->s = b.a + rotate(u * b.r, -th) - reso2->a;
|
|
}
|
|
cnt += 2;
|
|
} else if ((d - fabs(a.r - b.r)) / (a.r + b.r) > -eps) {
|
|
if (a.r / b.r < 1 - eps) {
|
|
if (reso1 != NULL) {
|
|
reso1->a = b.a - u * b.r;
|
|
reso1->s = rotate(u, pi / 2);
|
|
}
|
|
cnt++;
|
|
} else if (a.r / b.r > 1 + eps) {
|
|
if (reso1 != NULL) {
|
|
reso1->a = a.a + u * a.r;
|
|
reso1->s = rotate(u, pi / 2);
|
|
}
|
|
cnt++;
|
|
} else return -1;
|
|
}
|
|
if (d / (a.r + b.r) > 1 + eps) {
|
|
double th = acos((a.r + b.r) / d);
|
|
if (resi1 != NULL) {
|
|
resi1->a = a.a + rotate(u * a.r, th);
|
|
resi1->s = b.a - rotate(u * b.r, th) - resi1->a;
|
|
}
|
|
if (resi2 != NULL) {
|
|
resi2->a = a.a + rotate(u * a.r, -th);
|
|
resi2->s = b.a - rotate(u * b.r, -th) - resi2->a;
|
|
}
|
|
cnt += 2;
|
|
} else if (d / (a.r + b.r) > 1 - eps) {
|
|
if (resi1 != NULL) {
|
|
resi1->a = a.a + u * a.r;
|
|
resi1->s = rotate(u, pi / 2);
|
|
}
|
|
cnt++;
|
|
}
|
|
return cnt;
|
|
}
|
|
typedef std::vector<point2> convex2;
|
|
double area(const convex2 &a) {
|
|
double s = 0;
|
|
for (int i = 0; i < a.size(); i++)
|
|
s += a[i] % a[(i+1)%a.size()];
|
|
return fabs(s)/2;
|
|
}
|
|
int point_convex_inside(point2 a, const convex2 &b) {
|
|
double sum = 0;
|
|
for (int i = 0; i < b.size(); i++)
|
|
sum += fabs((b[i] - a) % (b[(i+1)%b.size()] - a));
|
|
return fabs(sum / (2*area(b)) - 1) < eps;
|
|
}
|
|
int line_convex_cross(line2 a, const convex2 &b, point2 *res1 = NULL, point2 *res2 = NULL) {
|
|
int cnt = 0;
|
|
for (int i = 0; i < b.size(); i++) {
|
|
line2 ltmp; point2 ptmp;
|
|
ltmp.a = b[i], ltmp.s = b[(i+1)%b.size()] - b[i];
|
|
int flag = line_line_cross(a, ltmp, &ptmp);
|
|
if (flag == -1) return -1;
|
|
if (flag == 0) continue;
|
|
double k = (ptmp - ltmp.a) * ltmp.s / (ltmp.s * ltmp.s);
|
|
if (k < eps || k > 1+eps) continue;
|
|
if (cnt == 0 && res1 != NULL) *res1 = ptmp;
|
|
else if (cnt == 1 && res2 != NULL) *res2 = ptmp;
|
|
cnt++;
|
|
}
|
|
return cnt;
|
|
}
|
|
int convex_gen_cmp(point2 a, point2 b) {
|
|
return a.y < b.y - eps || fabs(a.y - b.y) < eps && a.x < b.x - eps;
|
|
}
|
|
int convex_gen(const convex2 &a, convex2 &b) {
|
|
std::deque<point2> q;
|
|
convex2 t(a);
|
|
std::sort(t.begin(), t.end(), convex_gen_cmp);
|
|
q.push_back(t[0]), q.push_back(t[1]);
|
|
for (int i = 2; i < t.size(); i++) {
|
|
while (q.size() > 1) {
|
|
point2 p1 = t[i]-q[q.size()-1], p2 = q[q.size()-1]-q[q.size()-2];
|
|
if (p1 % p2 > 0 || parallel(p1, p2)) q.pop_back();
|
|
else break;
|
|
}
|
|
q.push_back(t[i]);
|
|
}
|
|
int pretop = q.size();
|
|
for (int i = t.size()-1; i >= 0; i--) {
|
|
while (q.size() > pretop) {
|
|
point2 p1 = t[i]-q[q.size()-1], p2 = q[q.size()-1]-q[q.size()-2];
|
|
if (p1 % p2 > 0 || parallel(p1, p2)) q.pop_back();
|
|
else break;
|
|
}
|
|
q.push_back(t[i]);
|
|
}
|
|
q.pop_back();
|
|
if (q.size() < 3) {
|
|
b.clear();
|
|
return 0;
|
|
}
|
|
b.clear();
|
|
for (int i = 0; i < q.size(); i++) b.push_back(q[i]);
|
|
return 1;
|
|
}
|
|
int halfplane_cross_cmp(line2 a, line2 b) {
|
|
double c1 = arg(a.s), c2 = arg(b.s);
|
|
return c1 < c2-eps || fabs(c1-c2) < eps && b.s % (a.a - b.a) / dis(b.s) > eps;
|
|
}
|
|
int halfplane_cross(const std::vector<line2> &a, convex2 &b) {
|
|
std::vector<line2> t(a);
|
|
std::sort(t.begin(), t.end(), halfplane_cross_cmp);
|
|
int j = 0;
|
|
for (int i = 0; i < t.size(); i++)
|
|
if (!i || arg(t[i].s) > arg(t[i-1].s) + eps) t[j++] = t[i];
|
|
if (j > 0 && arg(t[j].s) > arg(t[0].s) + 2*pi - eps) j--;
|
|
t.resize(j);
|
|
std::deque<line2> q;
|
|
q.push_back(t[0]), q.push_back(t[1]);
|
|
point2 p;
|
|
for (int i = 2, k = 0; i < t.size(); i++) {
|
|
for (; k < q.size() && t[i].s % q[k].s > 0; k++);
|
|
point2 s1 = q[q.size()-1].s, s2 = q[0].s;
|
|
if (k > 0 && k < q.size() && s1 % s2 > 0 && !parallel(s1, s2)) {
|
|
line_line_cross(q[k], q[k-1], &p);
|
|
double r1 = t[i].s % (p - t[i].a) / dis(t[i].s);
|
|
line_line_cross(q[0], q[q.size()-1], &p);
|
|
double r2 = t[i].s % (p - t[i].a) / dis(t[i].s);
|
|
if (r1 < eps && r2 < eps) {
|
|
b.clear();
|
|
return 0;
|
|
} else if (r1 > -eps && r2 > -eps)
|
|
continue;
|
|
}
|
|
while (q.size() > 1) {
|
|
line_line_cross(q[q.size()-1], q[q.size()-2], &p);
|
|
if (t[i].s % (p - t[i].a) / dis(t[i].s) < eps) {
|
|
q.pop_back();
|
|
if (k == q.size()) k--;
|
|
} else break;
|
|
}
|
|
while (q.size() > 1) {
|
|
line_line_cross(q[0], q[1], &p);
|
|
if (t[i].s % (p - t[i].a) / dis(t[i].s) < eps) {
|
|
q.pop_front();
|
|
k--; if (k < 0) k = 0;
|
|
} else break;
|
|
}
|
|
q.push_back(t[i]);
|
|
}
|
|
b.clear();
|
|
for (int i = 0; i < q.size(); i++) {
|
|
line_line_cross(q[i], q[(i+1)%q.size()], &p);
|
|
b.push_back(p);
|
|
}
|
|
return 1;
|
|
}
|
|
struct point3 {
|
|
double x, y, z;
|
|
point3& operator += (point3 a) { x+=a.x,y+=a.y,z+=a.z; return *this; }
|
|
point3& operator -= (point3 a) { x-=a.x,y-=a.y,z-=a.z; return *this; }
|
|
point3& operator *= (double a) { x*=a, y*=a, z*=a; return *this; }
|
|
point3& operator /= (double a) { x/=a, y/=a, z/=a; return *this; }
|
|
};
|
|
point3 operator + (point3 a, point3 b) { point3 c(a); c += b; return c; }
|
|
point3 operator - (point3 a, point3 b) { point3 c(a); c -= b; return c; }
|
|
point3 operator * (point3 a, double b) { point3 c(a); c *= b; return c; }
|
|
point3 operator * (double a, point3 b) { point3 c(b); c *= a; return c; }
|
|
point3 operator / (point3 a, double b) { point3 c(a); c /= b; return c; }
|
|
double operator * (point3 a, point3 b) { return a.x*b.x+a.y*b.y+a.z*b.z; }
|
|
point3 operator % (point3 a, point3 b) {
|
|
point3 c;
|
|
c.x = a.y * b.z - a.z * b.y;
|
|
c.y = a.z * b.x - a.x * b.z;
|
|
c.z = a.x * b.y - a.y * b.x;
|
|
return c;
|
|
}
|
|
double dis(point3 a) { return sqrt(a.x * a.x + a.y * a.y + a.z * a.z); }
|
|
int parallel(point3 a, point3 b) {
|
|
return a * a < eps * eps || b * b < eps * eps
|
|
|| (a % b) * (a % b) / ((a * a) * (b * b)) < eps * eps;
|
|
}
|
|
int perpend(point3 a, point3 b) {
|
|
return a * a < eps * eps || b * b < eps * eps
|
|
|| (a * b) * (a * b) / ((a * a) * (b * b)) < eps * eps;
|
|
}
|
|
struct line3 {
|
|
point3 a, s;
|
|
};
|
|
struct face3 {
|
|
point3 a, n;
|
|
};
|
|
struct circle3 {
|
|
point3 a;
|
|
double r;
|
|
};
|
|
point2 project(point3 a, face3 b, point3 xs) {
|
|
point3 ys; ys = b.n % xs;
|
|
point2 c; c.x = ((a - b.a) * xs) / dis(xs), c.y = ((a - b.a) * ys) / dis(ys);
|
|
return c;
|
|
}
|
|
line2 project(line3 a, face3 b, point3 xs) {
|
|
line2 c; c.a = project(a.a, b, xs), c.s = project(a.a + a.s, b, xs) - c.a;
|
|
return c;
|
|
}
|
|
point3 revproject(point2 a, face3 b, point3 xs) {
|
|
point3 ys; ys = b.n % xs;
|
|
return a.x * xs / dis(xs) + a.y * ys / dis(ys) + b.a;
|
|
}
|
|
line3 revproject(line2 a, face3 b, point3 xs) {
|
|
line3 c; c.a = revproject(a.a, b, xs), c.s = revproject(a.a + a.s, b, xs) - c.a;
|
|
return c;
|
|
}
|
|
double point_line_dis(point3 a, line3 b, point3 *res = NULL) {
|
|
point3 p;
|
|
p = b.a + ((a - b.a) * b.s) / (b.s * b.s) * b.s;
|
|
if (res != NULL) *res = p;
|
|
return dis(a - p);
|
|
}
|
|
double point_face_dis(point3 a, face3 b, point3 *res = NULL) {
|
|
point3 p;
|
|
p = ((a - b.a) * b.n) / (b.n * b.n) * b.n;
|
|
if (res != NULL) *res = a - p;
|
|
return dis(p);
|
|
}
|
|
double line_line_dis(line3 a, line3 b, point3 *resa = NULL, point3 *resb = NULL) {
|
|
point3 p;
|
|
if (parallel(a.s, b.s)) {
|
|
double d = point_line_dis(a.a, b, &p);
|
|
if (resa != NULL) *resa = a.a;
|
|
if (resb != NULL) *resb = p;
|
|
return d;
|
|
}
|
|
point3 n = a.s % b.s;
|
|
face3 f; f.a = b.a, f.n = n;
|
|
double d = point_face_dis(a.a, f, &p);
|
|
double k1 = ((b.a - p) % b.s) * n / (n * n);
|
|
if (resb != NULL) *resb = p + k1 * a.s;
|
|
if (resa != NULL) *resa = a.a + k1 * a.s;
|
|
return d;
|
|
}
|
|
int line_face_cross(line3 a, face3 b, point3 *res = NULL) {
|
|
if (perpend(a.s, b.n))
|
|
if (perpend(b.a - a.a, b.n))
|
|
return -1;
|
|
else
|
|
return 0;
|
|
double k = (b.a - a.a) * b.n / (a.s * b.n);
|
|
if (res != NULL) *res = a.a + k * a.s;
|
|
return 1;
|
|
}
|
|
int face_face_cross(face3 a, face3 b, line3 *res = NULL) {
|
|
if (parallel(a.n, b.n))
|
|
if (perpend(b.a - a.a, a.n))
|
|
return -1;
|
|
else
|
|
return 0;
|
|
point3 s = a.n % b.n;
|
|
point3 p;
|
|
line3 t; t.a = a.a, t.s = a.n % s;
|
|
line_face_cross(t, b, &p);
|
|
if (res != NULL)
|
|
res->a = p, res->s = s;
|
|
return 1;
|
|
}
|