/***********三分法求函数极值*************/
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; }
/*
向量叉积
若 P × Q > 0 , 则P在Q的顺时针方向。
若 P × Q < 0 , 则P在Q的逆时针方向。
若 P × Q = 0 , 则P与Q共线，但可能同向也可能反向。
*/
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); }
Point point(double a){
return p+(v*a);
}
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;
}
//多边形重心 点集逆时针给出
Point PolyGravity(Point *p, int n) {
Point tmp, g = Point(0, 0);
double sumArea = 0, area;
for (int i=2; i<n; ++i) {
area = Cross(p[i-1]-p[0], p[i]-p[0]);
sumArea += area;
tmp.x = p[0].x + p[i-1].x + p[i].x;
tmp.y = p[0].y + p[i-1].y + p[i].y;
g.x += tmp.x * area;
g.y += tmp.y * area;
}
g.x /= (sumArea * 3.0); g.y /= (sumArea * 3.0);
return g;
}
// 点集凸包
// 如果不希望在凸包的边上有输入点，把两个 <= 改成 <
// 注意：输入点集会被修改
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 isPointInPolygon(Point p, Polygon poly)
{
int wn = 0;
int n = poly.size();
for (int i = 0; i < n; i++)
{
if (OnSegment(p, poly[i], poly[(i + 1) % n])) return -1; //边界
int k = dcmp(Cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
int d1 = dcmp(poly[i].y - p.y);
int d2 = dcmp(poly[(i + 1) % n].y - p.y);
if (k > 0 && d1 <= 0 && d2 > 0) wn++;
if (k < 0 && d2 <= 0 && d1 > 0) wn--;
}
if (wn != 0) return 1; //内部
return 0; //外部
}
// 凸包直径返回 点集直径的平方
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);
}
};
//求a点到b点(逆时针)在的圆上的圆弧长度
double D(Point a,Point b,int id)
{
double ang1,ang2;
Vector v1,v2;
v1=a-Point(C[id].c.x,C[id].c.y);
v2=b-Point(C[id].c.x,C[id].c.y);
ang1=atan2(v1.y,v1.x);
ang2=atan2(v2.y,v2.x);
if(ang2<ang1) ang2+=2*pi;
return C[id].r*(ang2-ang1);
}
//直线与圆交点 返回个数
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;
}
//P到圆的切线
//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;
}
}
//两圆的公切线, -1表示无穷条切线
int getTangents(Circle A, Circle B, Point* a, Point* b)
{
int cnt = 0;
if (A.r < B.r) swap(A, B), swap(a, b);
///****************************
int d2 = (A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y);
int rdiff = A.r - B.r;
int rsum = A.r + B.r;
if (d2 < rdiff * rdiff) return 0; //内含
///***************************************
double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
if (d2 == 0 && A.r == B.r) return -1; //无线多条切线
if (d2 == rdiff * rdiff) //内切, 1条切线
{
///**********************
a[cnt] = A.point(base); b[cnt] = B.point(base); cnt++;
return 1;
}
//有外公切线
double ang = acos((A.r - B.r) / sqrt(d2));
a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;
a[cnt] = A.point(base - ang); b[cnt] = B.point(base - ang); cnt++;
if (d2 == rsum * rsum) //一条内公切线
{
a[cnt] = A.point(base); b[cnt] = B.point(PI + base); cnt++;
}
else if (d2 > rsum * rsum) //两条内公切线
{
double ang = acos((A.r + B.r) / sqrt(d2));
a[cnt] = A.point(base + ang); b[cnt] = B.point(PI + base + ang); cnt++;
a[cnt] = A.point(base - ang); b[cnt] = B.point(PI + base - ang); cnt++;
}
return cnt;
}
//三角形外接圆
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 InCircle(Point x, Circle c) { return dcmp(c.r - Length(c.c - x)) >= 0;}
bool OnCircle(Point x, Circle c) { return dcmp(c.r - 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);
}
/*********其他模板*********/
//以下模板来自网上，都未使用过
//pick定理
LL x_mult(cpoint a,cpoint b,cpoint p)
{
return (a.x-p.x)*(b.y-p.y)-(a.y-p.y)*(b.x-p.x);
}
LL pick()
{
LL s =0, e = 0;
for(int i=0;i<n;i++)
{
s += x_mult(re[i],re[i+1],re[0]);
e += gcd(abs(re[i].y-re[i+1].y),abs(re[i].x-re[i+1].x));
}e/=2;s/=2;
return Abs(s)+1-e;
}
//快速判断点是否在凸包内
struct POINT{
double x,y;
POINT(double _x = 0, double _y = 0):x(_x),y(_y){};
void show(){
cout<<x<<" "<<y<<endl;
}
};
POINT p[MAXN],wp[MAXN];
double multiply(POINT sp,POINT ep,POINT op){ //叉积
return (sp.x-op.x) * (ep.y-op.y) - (ep.x-op.x) * (sp.y-op.y);
}
bool onseg(POINT a,POINT s,POINT e){ // 判断点是否在线段上
if(multiply(a,s,e) == 0 && a.x <= max(s.x,e.x) && a.x >= min(s.x,e.x)
&& a.y <= max(s.y,e.y) && a.y >= min(s.y,e.y))
return true;
return false;
}
bool inside(POINT pp,POINT sp,POINT ep,POINT op){ //判断点pp是否在三角形中(极角序)
if(onseg(pp,sp,ep) || onseg(pp,sp,op) || onseg(pp,ep,op)) //如果在三角形上
return true;
if(multiply(sp,ep,pp) > 0 && multiply(ep,op,pp) > 0
&& multiply(sp,op,pp) < 0) //如果在三角形内
return true;
return false;
}
bool bsearch(POINT a,int len)
{ //二分所构造的三角形
int l = 1,r = len,m;
while(l < r){
m = (l + r) / 2;
if(inside(a,p[0],p[m],p[m+1]) == true) return true;
if(multiply(p[0],p[m],a) >= 0 && multiply(p[0],p[m+1],a) <= 0
&& multiply(p[m],p[m+1],a) < 0) return false;
if(multiply(p[0],p[m],a) > 0 && multiply(p[0],p[m+1],a) > 0)
l = m + 1;
else r = m;
}
return false;
}
int main()
{
int n,m,k,tmp = 0,cnt = 0;
scanf("%d%d%d",&n,&m,&k);
for(int i = 0 ; i < n ; i++)
scanf("%lf%lf",&p[i].x,&p[i].y);
p[n] = p[0];
for(int i = 0 ; i < m ; i++){
scanf("%lf%lf",&wp[i].x,&wp[i].y);
if(bsearch(wp[i],n-1) == true) cnt++;
}
if(cnt >= k) printf("YES\n");
else printf("NO\n");
return 0;
}
//圆的面积并
//圆的面积并
//
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#define sqr(x) ((x)*(x))
using namespace std;
const int N = 1010;
const double eps = 1e-8;
const double pi = acos(-1.0);
double area[N];
int n;
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[N], tp[N * 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]);
}
}
}
void solve()
{
scanf("%d", &n);
for (int i = 0; i < n; ++i)
cir[i].get();
memset(area, 0, sizeof(area));
CirUnion(cir, n);
//去掉重复计算的
for (int i = 1; i <= n; ++i) {
area[i] -= area[i + 1];
}
//area[i]为重叠了i次的面积
//tot 为总面积
double tot = 0;
for(int i=1; i<=n; i++) tot += area[i];
printf("%f\n", tot);
}
//多边形面积并
#define PDI pair<double,int>
#define point pair<double,double>
#define mp make_pair
#define pb push_back
#define x first
#define y second
#define zero 1e-8
#define maxN 502
#define maxp 5
point operator +(point a,point b) { return mp(a.x+b.x,a.y+b.y); }
point operator -(point a,point b) { return mp(a.x-b.x,a.y-b.y); }
double operator *(point a,point b) { return a.x*b.y-b.x*a.y; }
double operator ^(point a,point b) { return a.x*b.x+a.y*b.y; }
inline double cross(point o,point a,point b) { return (a-o)*(b-o); }
inline int cmp(double x) { if (fabs(x)<zero) return 0; return x>0? 1:-1; }
class Polygon
{
private: int i; double s;
public: int n; point p[maxp];
point& operator[] (int idx) { return p[idx]; }
void input() { for (i=0;i<n;i++) scanf("%lf %lf",&p[i].x,&p[i].y); p[n]=p[0]; }
double Area() { for (s=0,i=0;i<n;i++) s+=p[i]*p[i+1]; return s/2; }
};
PDI s[maxN*maxp*2];
Polygon P[maxN];
double S,ts;
int N;
inline double seg(point o,point a,point b)
{
if (cmp(b.x-a.x)==0) return (o.y-a.y)/(b.y-a.y);
return (o.x-a.x)/(b.x-a.x);
}
double PolygonUnion()
{
int M,c1,c2; double s1,s2,ret=0;
for (int i=0;i<N;i++)
for (int ii=0;ii<P[i].n;ii++)
{
M=0;
s[M++]=mp(0.00,0);
s[M++]=mp(1.00,0);
for (int j=0;j<N;j++) if (i!=j)
for (int jj=0;jj<P[j].n;jj++)
{
c1=cmp(cross(P[i][ii],P[i][ii+1],P[j][jj]));
c2=cmp(cross(P[i][ii],P[i][ii+1],P[j][jj+1]));
if (c1==0 && c2==0)
{
if (((P[i][ii+1]-P[i][ii])^(P[j][jj+1]-P[j][jj]))>0 && i>j)
{
s[M++]=mp(seg(P[j][jj],P[i][ii],P[i][ii+1]),1);
s[M++]=mp(seg(P[j][jj+1],P[i][ii],P[i][ii+1]),-1);
}
}
else
{
s1=cross(P[j][jj],P[j][jj+1],P[i][ii]);
s2=cross(P[j][jj],P[j][jj+1],P[i][ii+1]);
if (c1>=0 && c2<0) s[M++]=mp(s1/(s1-s2),1);
else if (c1<0 && c2>=0) s[M++]=mp(s1/(s1-s2),-1);
}
}
sort(s,s+M);
double pre=min(max(s[0].x,0.0),1.0),now;
double sum=0;
int cov=s[0].y;
for (int j=1;j<M;j++)
{
now=min(max(s[j].x,0.0),1.0);
if (!cov) sum+=now-pre;
cov+=s[j].y;
pre=now;
}
ret+=P[i][ii]*P[i][ii+1]*sum;
}
return ret/2;
}
int main()
{
scanf("%d\n",&N);
for (int i=0;i<N;i++)
{
P[i].n=4;
P[i].input();
ts=P[i].Area();
if (cmp(ts<0))
{
reverse(P[i].p,P[i].p+P[i].n);
P[i][P[i].n]=P[i][0];
ts=-ts;
}
S+=ts;
}
printf("%.9lf\n",S/PolygonUnion());
}
//二维平面一个n个节点的简单多边形，多边形内有一个灯泡，求照明面积。
const double eps = 1e-8;
const double pi = acos(-1.);
using namespace std;
int dblcmp( double x )
{
if( fabs(x) < eps ) return 0;
return x > 0 ? 1 : -1;
}
double nowAng;
struct point
{
double x, y, a;
point(){}
point( double _x, double _y ) : x(_x), y(_y)
{
a = atan2(y, x);
}
bool operator<( const point p ) const
{
if( dblcmp(a-p.a) == 0 )
return x*x+y*y < p.x*p.x+p.y*p.y;
return a < p.a;
}
} p[60000], O;
inline double dis( point a, point b )
{
double dx = a.x-b.x;
double dy = a.y-b.y;
return sqrt(dx*dx+dy*dy);
}
inline double cross( point k, point a, point b )
{
return (a.x-k.x)*(b.y-k.y) - (a.y-k.y)*(b.x-k.x);
}
point inter( point a1, point a2, double ang )
{
point b1 = point(0, 0), b2 = point(cos(ang), sin(ang));
double u = cross(a1, a2, b1), v = cross(a2, a1, b2);
return point((b1.x*v+b2.x*u)/(v+u), (b1.y*v+b2.y*u)/(v+u));
}
struct line
{
point a, b;
line(){};
line( point _a, point _b ) : a(_a), b(_b){};
bool operator<( const line p ) const
{
if( fabs(a.x-p.a.x) < eps && fabs(a.y-p.a.y) < eps )
return cross(a, b, p.b) < 0;
point d1 = inter(a, b, nowAng);
point d2 = inter(p.a, p.b, nowAng);
return d1.x*d1.x+d1.y*d1.y < d2.x*d2.x+d2.y*d2.y;
}
};
struct Event
{
double ang;
int id, st;
line L;
bool operator<( const Event p ) const
{
if( dblcmp(ang-p.ang) == 0 )
return st > p.st;
return ang < p.ang;
}
} E[200000];
int c;
void add( point a, point b, int k )
{
if( b < a )
swap(a, b);
E[c].ang = a.a, E[c].st = 1, E[c].L = line(a, b), E[c++].id = k;
E[c].ang = b.a, E[c].st = 0, E[c++].id = k;
}
multiset<line> S;
multiset<line>::iterator itArr[100000];
inline double cal( line L, double d1, double d2 )
{
point a = inter(L.a, L.b, d1);
point b = inter(L.a, L.b, d2);
return fabs(0.5*cross(O, a, b));
}
int main()
{
int i, j, k, n;
double d, t, ans, pre;
while( scanf("%lf %lf", &O.x, &O.y) != EOF )
{
c = 0;
scanf("%d", &n);
for( i = 0; i < n; ++i )
{
scanf("%lf %lf", &p[i].x, &p[i].y);
p[i].x -= O.x, p[i].y -= O.y;
p[i].a = atan2(p[i].y, p[i].x);
}
O.x = O.y = 0;
p[n] = p[0];
for( i = k = 0; i < n; ++i )
{
d = fabs(p[i+1].a-p[i].a);
if( d < pi )
add(p[i], p[i+1], k++);
else
{
point tmp = inter(p[i], p[i+1], pi);
tmp.a = pi*dblcmp(p[i].a);
add(p[i], tmp, k++);
tmp.a = pi*dblcmp(p[i+1].a);
add(p[i+1], tmp, k++);
}
}
sort(E, E+c);
S.clear();
ans = 0; pre = -pi;
for( i = 0; i < c; ++i )
{
nowAng = E[i].ang;
if( E[i].st )
{
if( S.size() > 0 )
ans += cal(*S.begin(), pre, E[i].ang);
itArr[E[i].id] = S.insert(E[i].L);
}
else
{
ans += cal(*S.begin(), pre, E[i].ang);
S.erase(itArr[E[i].id]);
}
pre = E[i].ang;
}
printf("%.10lf\n", ans);
}
return 0;
}
//二维平面有n(0 < n <= 50000)条线段，要求判断n条线段是否存在交点，如果有，输出相交线段的编号。
#define MP make_pair
#define PI pair
#define FI first
#define SE second
#define PB push_back
#define SZ size()
const double eps = 1e-10;
const double pi = acos(-1.);
const int mod = 1000000007;
const int maxn = 50100;
const int INF = 99999999;
struct point
{
int x, y;
bool operator<( const point p ) const
{
if( x == p.x ) return y < p.y;
return x < p.x;
}
} L[maxn][2];
int tim;
struct eve
{
int x, id, st;
eve(){}
eve( int _x, int _id, int _st ) : x(_x), id(_id), st(_st) {}
bool operator<( const eve p ) const
{
if( p.x == x ) return st > p.st;
return x < p.x;
}
} E[maxn*2];
void getLine( point x, point y, double& a, double& b, double& c )
{
a = y.y - x.y;
b = x.x - y.x;
c = y.x*x.y - x.x*y.y;
}
struct ele
{
int id;
double k, c;
ele( int _id )
{
id = _id;
if( L[id][0].x == L[id][1].x )
k = 0, c = L[id][0].y;
else
{
double A, B, C;
getLine(L[id][0], L[id][1], A, B, C);
k = -A/B, c = -C/B;
}
}
bool operator<( const ele p ) const
{
return tim*k+c < tim*p.k+p.c;
}
};
set<ele> S;
set<ele>::iterator itArr[maxn];
inline set<ele>::iterator preIt( set<ele>::iterator it )
{
return it == S.begin() ? S.end() : --it;
}
inline set<ele>::iterator nxtIt( set<ele>::iterator it )
{
return it == S.end() ? S.end() : ++it;
}
int cross( point& k, point& a, point& b )
{
return (a.x-k.x)*(b.y-k.y) - (a.y-k.y)*(b.x-k.x);
}
int dot( point& k, point& a, point& b )
{
return (a.x-k.x)*(b.x-k.x) + (a.y-k.y)*(b.y-k.y);
}
inline int sgn( int x )
{
if( x > 0 ) return 1;
if( x < 0 ) return -1;
return x;
}
bool inter( int a, int b )
{
int d1 = sgn(cross(L[a][0], L[a][1], L[b][0]));
int d2 = sgn(cross(L[a][0], L[a][1], L[b][1]));
int d3 = sgn(cross(L[b][0], L[b][1], L[a][0]));
int d4 = sgn(cross(L[b][0], L[b][1], L[a][1]));
if( (d1^d2)==-2 && (d3^d4)==-2 ) return 1;
if( d1 == 0 && dot(L[b][0], L[a][0], L[a][1]) <= 0 ) return 1;
if( d2 == 0 && dot(L[b][1], L[a][0], L[a][1]) <= 0 ) return 1;
if( d3 == 0 && dot(L[a][0], L[b][0], L[b][1]) <= 0 ) return 1;
if( d4 == 0 && dot(L[a][1], L[b][0], L[b][1]) <= 0 ) return 1;
return 0;
}
bool solve( int n )
{
sort(E, E+n);
S.clear();
for( int i = 0; i < n; ++i )
{
tim = E[i].x;
int id = E[i].id;
if( E[i].st == 1 )
{
ele t = ele(id);
set<ele>::iterator nxt = S.lower_bound(t), pre = preIt(nxt);
if( nxt != S.end() && inter((*nxt).id, id) )
{
printf("YES\n%d %d\n", (*nxt).id, id);
return 1;
}
if( pre != S.end() && inter((*pre).id, id) )
{
printf("YES\n%d %d\n", (*pre).id, id);
return 1;
}
itArr[id] = S.insert(nxt, t);
}
else
{
set<ele>::iterator pre = preIt( itArr[id] ), nxt = nxtIt( itArr[id] );
if( pre != S.end() && nxt != S.end() && inter((*pre).id, (*nxt).id) )
{
printf("YES\n%d %d\n", (*pre).id, (*nxt).id);
return 1;
}
S.erase(itArr[id]);
}
}
return 0;
}
int main()
{
srand(4);
int T, cases = 1;
int i, j, k, e;
double A, B, C;
int N, M;
scanf("%d", &N);
for( i = 1; i <= N; ++i )
{
scanf("%d %d %d %d", &L[i][0].x, &L[i][0].y, &L[i][1].x, &L[i][1].y);
if( L[i][1] < L[i][0] )
swap(L[i][0], L[i][1]);
}
for( i = 1, e = 0; i <= N; ++i )
{
E[e++] = eve(L[i][0].x, i, +1);
E[e++] = eve(L[i][1].x, i, -1);
}
if( !solve(e) )
puts("NO");
return 0;
}
/*******三维几何*******/
#include <iostream>
#include <cstdio>
#include <cmath>
#define zero( x ) ( ( ( x ) > 0 ? ( x ) : -( x ) ) < eps )
const double eps = 1e-8;
struct point3
{
double x, y, z;
};
struct line3
{
point3 a, b;
};
struct plane3
{
point3 a, b, c;
};
//计算 Cross Product U * V
point3 xmult( point3 u, point3 v )
{
point3 ret;
ret.x = u.y * v.x - v.y * u.z;
ret.y = u.z * v.x - u.x * v.z;
ret.z = u.x * v.y - u.y * v.x;
return ret;
}
//计算 Dot Product U ? V
double dmult( point3 u, point3 v )
{
return u.x * v.x + u.y * v.y + u.z * v.z;
}
//矢量差 U - V
point3 subt( point3 u, point3 v )
{
point3 ret;
ret.x = u.x - v.x;
ret.y = u.y - v.y;
ret.z = u.z - v.z;
return ret;
}
//取平面法向量
point3 pvec( plane3 s )
{
return xmult( subt( s.a, s.b ), subt( s.b, s.c ) );
}
point3 pvec( point3 s1, point3 s2, point3 s3 )
{
return xmult( subt( s1, s2 ), subt( s2, s3 ) );
}
//两点距离，单参数取向量大小
double distance( point3 p1, point3 p2 )
{
return sqrt( ( p1.x - p2.x ) * ( p1.x - p2.x ) + ( p1.y - p2.x ) * ( p1.y - p2.y ) + ( p1.z - p2.z ) * ( p1.z - p2.z ) );
}
//向量大小
double vlen( point3 p )
{
return sqrt( p.x * p.x + p.y * p.y + p.z * p.z );
}
//判定三点共线
int dots_inline( point3 p1, point3 p2, point3 p3 )
{
return vlen( xmult( subt( p1, p2 ), subt( p2, p3 ) ) ) < eps;
}
//判定四点共面
int dots_onplane( point3 a, point3 b, point3 c, point3 d )
{
return zero( dmult( pvec( a, b, c ), subt( d, a ) ) );
}
//判定点是否在线段上，包括端点和共线
int dot_online_in( point3 p, line3 l )
{
return zero( vlen( xmult( subt( p, l.a ), subt( p, l.b ) ) ) ) && ( l.a.x - p.x ) * ( l.b.x - p.x ) < eps && ( l.a.y - p.y ) * ( l.b.y - p.y ) < eps && ( l.a.z - p.z ) * ( l.b.z - p.z ) < eps;
}
int dot_online_in( point3 p, point3 l1, point3 l2 )
{
return zero( vlen( xmult( subt( p, l1 ), subt( p, l2 ) ) ) ) && ( l1.x - p.x ) * ( l2.x - p.x ) < eps && ( l1.y - p.y ) * ( l2.y - p.y ) < eps && ( l1.z - p.z ) * ( l2.z - p.z ) < eps;
}
//判定点是否在线段上，不包括端点
int dot_online_ex( point3 p, line3 l )
{
return dot_online_in( p, l ) && ( !zero( p.x - l.a.x ) || ! zero( p.y - l.a.y ) || !zero( p.z - l.a.z ) ) && ( !zero( p.x - l.b.x ) || !zero( p.y - l.b.y ) || !zero( p.z - l.b.z ) );
}
int dot_online_ex( point3 p, point3 l1, point3 l2 )
{
return dot_online_in( p, l1, l2 ) && ( !zero( p.x - l1.x ) || !zero( p.y - l1.y ) || !zero( p.z - l1.z ) ) && ( !zero( p.x - l2.x ) || !zero( p.y - l2.y ) || !zero( p.z - l2.z ) );
}
//判定点是否在空间三角形上，包括边界，三点共线无意义
int dot_inplane_in( point3 p, plane3 s )
{
return zero( vlen( xmult( subt( s.a, s.b ), subt( s.a, s.c ) ) ) - vlen( xmult( subt( p, s.a ), subt( p, s.b ) ) ) - vlen( xmult( subt( p, s.b ), subt( p, s.c ) ) ) - vlen( xmult( subt( p, s.c ), subt( p, s.a ) ) ) );
}
int dot_inplane_in( point3 p, point3 s1, point3 s2, point3 s3 )
{
return zero( vlen( xmult( subt( s1, s2 ), subt( s1, s3 ) ) ) - vlen( xmult( subt( p, s1 ), subt( p, s2 ) ) ) - vlen( xmult( subt( p, s2 ), subt( p, s3 ) ) ) - vlen( xmult( subt( p, s3 ), subt( p, s1 ) ) ) );
}
//判定点是否在空间三角形上，不包括边界，三点共线无意义
int dot_inplane_ex( point3 p, plane3 s )
{
return dot_inplane_in( p, s ) && vlen( xmult( subt( p, s.a ), subt( p, s.b ) ) ) > eps && vlen( xmult( subt( p, s.b ), subt( p, s.c ) ) ) > eps && vlen( xmult( subt( p, s.c ), subt( p, s.a ) ) ) > eps;
}
int dot_inplane_ex( point3 p, point3 s1, point3 s2, point3 s3 )
{
return dot_inplane_in( p, s1, s2, s3 ) && vlen( xmult( subt( p, s1 ), subt( p, s2 ) ) ) > eps && vlen( xmult( subt( p, s2 ), subt( p, s3 ) ) ) > eps && vlen( xmult( subt( p, s3 ), subt( p, s1 ) ) ) > eps;
}
//判定两点在线段同侧，点在线段上返回0，不共面无意义
int same_side( point3 p1, point3 p2, line3 l )
{
return dmult( xmult( subt( l.a, l.b ), subt( p1, l.b ) ), xmult( subt( l.a, l.b ), subt( p2, l.b ) ) ) > eps;
}
int same_side( point3 p1, point3 p2, point3 l1, point3 l2 )
{
return dmult( xmult( subt( l1, l2 ), subt( p1, l2 ) ), xmult( subt( l1, l2 ), subt( p2, l2 ) ) ) > eps;
}
//判定两点在线段异侧，点在平面上返回0
int opposite_side( point3 p1, point3 p2, line3 l )
{
return dmult( xmult( subt( l.a, l.b ), subt( p1, l.b ) ), xmult( subt( l.a, l.b ), subt( p2, l.b ) ) ) < -eps;
} int opposite_side( point3 p1, point3 p2, point3 l1, point3 l2 )
{
return dmult( xmult( subt( l1, l2 ), subt( p1, l2 ) ), xmult( subt( l1, l2 ), subt( p2, l2 ) ) ) < -eps;
}
//判定两点在平面同侧，点在平面上返回0
int same_side( point3 p1, point3 p2, plane3 s )
{
return dmult( pvec( s ), subt( p1, s.a ) ) * dmult( pvec( s ), subt( p2, s.a ) ) > eps;
}
int same_side( point3 p1, point3 p2, point3 s1, point3 s2, point3 s3 )
{
return dmult( pvec( s1, s2, s3 ), subt( p1, s1 ) ) * dmult( pvec( s1, s2, s3 ), subt( p2, s1 ) ) > eps;
}
//判定两点在平面异侧，点在平面上返回0
int opposite_side( point3 p1, point3 p2, plane3 s )
{
return dmult( pvec( s ), subt( p1, s.a ) ) * dmult( pvec( s ), subt( p2, s.a ) ) < -eps;
}
int opposite_side( point3 p1, point3 p2, point3 s1, point3 s2, point3 s3 )
{
return dmult( pvec( s1, s2, s3 ), subt( p1, s1 ) ) * dmult( pvec( s1, s2, s3 ), subt( p2, s1 ) ) < eps;
}
//判定两直线平行
int parallel( line3 u, line3 v )
{
return vlen( xmult( subt( u.a, u.b ), subt( v.a, v.b ) ) ) < eps;
}
int parallel( point3 u1, point3 u2, point3 u3, point3 v1, point3 v2, point3 v3 )
{
return vlen( xmult( pvec( u1, u2, u3 ), pvec( v1, v2, v3 ) ) ) < eps;
}
//判定直线与平面平行
int parallel( line3 l, plane3 s )
{
return zero( dmult( subt( l.a, l.b ), pvec( s ) ) );
}
int parallel( point3 l1, point3 l2, point3 s1, point3 s2, point3 s3 )
{
return zero( dmult( subt( l1, l2 ), pvec( s1, s2, s3 ) ) );
}
//判定两直线垂直
int perpendicular( line3 u, line3 v )
{
return zero( dmult( subt( u.a, u.b ), subt( v.a, v.b ) ) );
}
int perpendicular( point3 u1, point3 u2, point3 v1, point3 v2 )
{
return zero( dmult( subt( u1, u2 ), subt( v1, v2 ) ) );
}
//判定两平面垂直
int perpendicular( plane3 u, plane3 v )
{
return zero( dmult( pvec( u ), pvec( v ) ) );
}
int perpendicular( point3 u1, point3 u2, point3 u3, point3 v1, point3 v2, point3 v3 )
{
return zero( dmult( pvec( u1, u2, u3 ), pvec( v1, v2, v3 ) ) );
}
//判定直线与平面平行
int perpendicular( line l, plane3 s )
{
return vlen( xmult( subt( l.a, l.b ), pvec( s ) ) ) < eps;
}
int perpendicular( point3 l1, point3 l2, point3 s1, point3 s2, point3 s3 )
{
return vlen( xmult( subt( l1, l2 ), pvec( s1, s2, s3 ) ) ) < eps;
}
//判定两直线相交，包括端点和部分重合
int intersect_in( line3 u, line3 v )
{
if( !dots_onplane( u.a, u.b, v.a, v.b ) )
return 0;
if( !dotsinline( u.a, u.b, v.a ) || !dots_inline( u.a, u.b, v.b ) )
return !same_side( u.a, u.b, v ) && !same_side( v.a, v.b, u );
return dot_online_in( u.a, v ) || dot_online_in( u.b, v ) || dot_online_in( v.a, u ) || dot_online_in( v.b, u );
}
int intersect_in( point3 u1, point3 u2, opint3 v1, point3 v2 )
{
if( !dots_onplane( u1, u2, v1, v2 ) )
return 0;
if( !dots_inline( u1, u2, v1 ) || !dots_inline( u1, u2, v2 ) )
return !same_side( u1, u2, v1, v2 ) && !same_side( v1, v2, u1, u2 );
return dot_online_in( u1, v1, v2 ) || dot_online_in( u2, v1, v2 ) || dot_online_in( v1, u1, u2 ) || dot_online_in( v2, u1, u2 );
}
//判定两线段相交，不包括端点和部分重合
int intersect_ex( line3 u, line3 v )
{
return dots_onplane( u.a, u.b, v.a, v.b ) && opposite_side( u.a, u.b, v ) && opposite_side( v.a, v.b, u );
}
int intersect_ex( point3 u1, point3 u2, point3 v1, point3 v2 )
{
return dots_onplane( u1, u2, v1, v2 ) && opposite_side( u1, u2, v1, v2 ) && opposite_side( v1, v2, u1, u2 );
}
//判定线段与空间三角形相交，包括交于边界和（部分）包含
int intersect_in( line3 l, plane3 s )
{
return !same_side( l.a, l.b, s ) && !same_side( s.a, s.b, l.a, l.b, s.c ) && !same_side( s.b, s.c, l.a, l.b, s.a ) && !same_side( s.c, s.a, l.a, l.b, s.b );
}
int intersect_in( point3 l1, point3 l2, point3 s1, point3 s2, point3 s3 )
{
return !same_side( l1, l2, s1, s2, s3 ) && !same_side( s1, s2, l1, l2, s3 ) && !same_side( s2, s3, l1, l2, s1 ) && !same_side( s3, s1, l1, l2, s2 );
}
//判定线段与空间三角形相交，不包括交于边界和（部分）包含
int intersect_ex( line3 l, plane3 s )
{
return opposite_side( l.a, l.b, s ) && opposite_side( s.a, s.b, l.a, l.b, s.c ) && opposite_side( s.b, s.c, l.a, l.b, s.a ) && opposite_side( s.c, s.a, l.a, l.b, s.b );
}
int intersect_ex( point3 l1, point3 l2, point3 s1, point3 s2, point3 s3 )
{
return opposite_side( l1, l2, s1, s2, s3 ) && opposite_side( s1, s2, l1, l2, s3 ) && opposite_side( s2, s3, l1, l2, s1 ) && opposite_side( s3, s1, l1, l2, s2 );
}
//计算两直线交点，注意事先判断直线是否共面和平行
//线段交点另外判断线段相交（同时还是要判断是否平行）
point3 intersection( line3 u, line3 v )
{
point3 ret = u.a;
double t = ( ( u.a.x - v.a.x ) * ( v.a.y - v.b.y ) - ( u.a.y - v.a.y ) * ( v.a.x - v.b.x ) ) / ( ( u.a.x - u.b.x ) * ( v.a.y - v.b.y ) - ( u.a.y - u.b.y ) * ( v.a.x - v.b.x ) );
ret.x += ( u.b.x - u.a.x ) * t;
ret.y += ( u.b.y - u.a.y ) * t;
ret.z += ( u.b.z - u.a.z ) * t;
return ret;
}
point3 intersection( point3 u1, point3 u2, point3 v1, point3 v2 )
{
point3 ret = u1;
double t = ( ( u1.x - v1.x ) * ( v1.y - v2.y ) - ( u1.y - v1.y ) * ( v1.x - v2.x ) ) / ( ( u1.x - u2.x ) * ( v1.y - v2.y ) - ( u1.y - u2.y ) * ( v1.x - v2.x ) );
ret.x += ( u2.x - u1.x ) * t;
ret.y += ( u2.y - u1.y ) * t;
ret.z += ( u2.z - u1.z ) * t;
return ret;
}
//计算直线与平面交点，注意事先判断是否平行，并保证三点不共线
//线段和空间三角形交点另外判断
point3 intersection( line3 l, plane3 s )
{
point3 ret = pvec( s );
double t = ( ret.x * ( s.a.x - l.a.x ) + ret.y * ( s.a.y - l.a.y ) + ret.z * ( s.a.z - l.a.z ) ) / ( ret.x * ( l.b.x - l.a.x ) + ret.y * ( l.b.y - l.a.y ) + ret.z * ( l.b.z - l.a.z ) );
ret.x = l.a.x + ( l.b.x - l.a.x ) * t;
ret.y = l.a.y + ( l.b.y - l.a.y ) * t;
ret.z = l.a.z + ( l.b.z - l.a.z ) * t;
return ret;
}
point3 intersection( point3 l1, point3 l2, point3 s1, point3 s2, point3 s3 )
{
point3 ret = pvec( s1, s2, s3 );
double t = ( ret.x * ( s1.x - l1.x ) + ret.y * ( s1.y - l1.y ) + ret.z * ( s1.z - l1.z ) ) / ( ret.x * ( l2.x - l1.x ) + ret.y * ( l2.y - l1.y ) + ret.z * ( l2.z - l1.z ) );
ret.x = l1.x + ( l2.x - l1.x ) * t;
ret.y = l1.y + ( l2.y - l1.y ) * t;
ret.z = l1.z + ( l2.z - l1.z ) * t;
return ret;
}
//计算两平面交线，注意事先判断是否平行，并保证三点不共线
line3 intersection( plane3 u, plane3 v )
{
line3 ret;
ret.a = parallel( v.a, v.b, u.a, u.b, u.c ) ? intersection( v.b, v.c, u.a, u.b, u.c ) : intersection( v.a, v.b, u.a, u.b, u.c );
ret.b = parallel( v.c, v.a, u.a, u.b, u.c ) ? intersection( v.b, v.c, u.a, u.b, u.c ) : intersection( v.c, v.a, u.a, u.b, u.c );
return ret;
}
line3 intersection( point3 u1, point3 u2, point3 u3, point3 v1, point3 v2, point3 v3 )
{
line3 ret;
ret.a = parallel( v1, v2, u1, u2, u3 ) ? intersection( v2, v3, u1, u2, u3 ) : intersection( v1, v2, u1, u2, u3 );
ret.b = parallel( v3, v1, u1, u2, u3 ) ? intersection( v2, v3, u1, u2, u3 ) : intersection( v3, v1, u1, u2, u3 );
return ret;
}
//点到直线距离
double ptoline( point3 p, line3 l )
{
return vlen( xmult( subt( p, l.a ), subt( l.b, l.a ) ) ) / distance( l.a, l.b );
}
double ptoline( point3 p, opint3 l1, point3 l2 )
{
return vlen( xmult( subt( p, l1 ), subt( l2, l1 ) ) ) / distance( l1, l2 );
}
//点到平面距离
double ptoplane( point3 p, plane3 s )
{
return fabs( dmult( pvec( s ), subt( p, s.a ) ) ) / vlen( pvec( s ) );
}
double ptoplane( point3 p, point3 s1, point3 s2, point3 s3 )
{
return fabs( dmult( pvec( s1, s2, s3 ), subt( p, s1 ) ) ) / vlen( pvec( s1, s2, s3 ) );
}
//直线到直线距离
double linetoline( lien3 u, line3 v )
{
point3 n = xmult( subt( u.a, u.b ), subt( v.a, v.b ) );
return fabs( dmult( subt( u.a, v.a ), n ) ) / vlen( n );
}
double linetoline( point3 u1, point3 u2, point3 v1, point3 v2 )
{
point3 n = xmult( subt( u1, u2 ), subt( v1, v2 ) );
return fabs( dmult( subt( u1, v1 ), n ) ) / vlen( n );
}
//两直线夹角 cos 值
double angle_cos( line3 u, line3 v )
{
return dmult( subt( u.a, u.b ), subt( v.a, v.b ) ) / vlen( subt( u.a, u.b ) ) / vlen( subt( v.a, v.b ) );
}
double angle_cos( point3 u1, point3 u2, point3 v1, point3 v2 )
{
return dmult( subt( u1, u2 ), subt( v1, v2 ) ) / vlen( sut( u1, u2 ) ) / vlen( subt( v1, v2 ) );
}
//两平面夹角 cos 值
double angle_cos( plane3 u, plane3 v )
{
return dmult( pvec( u ), pvec( v ) ) / vlen( pvec( u ) ) / vlen( pvec( v ) );
}
double angle_cos( point3 u1, point3 u2, point3 v1, point3 v2, point3 v3 )
{
return dmult( pvec( u1, u2, u3 ), pvec( v1, v2, v3 ) ) / vlen( pvec( u1, u2, u3 ) ) / vlen( pvec( v1, v2, v3 ) );
}
//直线平面夹角 sin 值
double angle_sin( line3 l, plane3 s )
{
return dmult( subt( l.a, l.b ), pvec( s ) ) / vlen( subt( l.a, l.b ) ) / vlen( pvec( s ) );
}
double angle_sin( point3 l1, point3 l2, point3 s1, point3 s2, point3 s3 )
{
return dmult( subt( l1, l2 ), pvc( s1, s2, s3 ) ) / vlen( subt( l1, l2 ) ) / vlen( pvec( s1, s2, s3 ) );
}
扫一扫
04-08
2691
2691
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
