#include <iostream> #include <cmath> #include <vector> #include <stdio.h> #include <string.h> #include <stdlib.h> #include <algorithm> using namespace std; #define MAX_N 110 /*------------------常量区-------------------*/ const double INF = 1e10; // 无穷大 const double EPS = 1e-6; // 计算精度 const double PI = acos(-1.0);// PI const int PINXING = 0; // 平行 const int XIANGJIAO = 1; // 相交 const int XIANGLI = 0; // 相离 const int GONGXIAN = 2; // 共线 const int CHONGDIE = -1; // 重叠 const int INSIDE = 1; // 点在图形内部 const int OUTSIDE = 0; // 点在图形外部 const int BORDER = 2; // 点在图形边界 /*-----------------类型定义区----------------*/ struct Point { // 二维点或矢量 double x, y; //double angle, dis; Point() {} Point(double x0, double y0): x(x0), y(y0) {} void read() { scanf("%lf%lf",&x,&y); } }; struct Point3D { //三维点或矢量 double x, y, z; Point3D() {} Point3D(double x0, double y0, double z0): x(x0), y(y0), z(z0) {} }; struct Line { // 二维的直线或线段 Point p1, p2; Line() {} Line(Point p10, Point p20): p1(p10), p2(p20) {} void read() { scanf("%lf%lf%lf%lf",&p1.x,&p1.y,&p2.x,&p2.y); } }; struct Line3D { // 三维的直线或线段 Point3D p1, p2; Line3D() {} Line3D(Point3D p10, Point3D p20): p1(p10), p2(p20) {} }; struct Rect { // 用长宽表示矩形的方法 w, h分别表示宽度和高度 double w, h; Rect() {} Rect(double _w,double _h) : w(_w),h(_h) {} }; struct Rect_2 { // 表示矩形，左下角坐标是(xl, yl)，右上角坐标是(xh, yh) double xl, yl, xh, yh; Rect_2() {} Rect_2(double _xl,double _yl,double _xh,double _yh) : xl(_xl),yl(_yl),xh(_xh),yh(_yh) {} }; struct Circle { //圆 Point c; double r; Circle() {} Circle(Point _c,double _r) :c(_c),r(_r) {} }; typedef vector<Point> Polygon; // 二维多边形 typedef vector<Point> Points; // 二维点集 /*-------------------基本函数区---------------------*/ inline double max(double x,double y) { return x > y ? x : y; } inline double min(double x, double y) { return x > y ? y : x; } inline bool ZERO(double x) // x == 0 { return (fabs(x) < EPS); } inline bool ZERO(Point p) // p == 0 { return (ZERO(p.x) && ZERO(p.y)); } inline bool ZERO(Point3D p) // p == 0 { return (ZERO(p.x) && ZERO(p.y) && ZERO(p.z)); } inline bool EQ(double x, double y) // eqaul, x == y { return (fabs(x - y) < EPS); } inline bool NEQ(double x, double y) // not equal, x != y { return (fabs(x - y) >= EPS); } inline bool LT(double x, double y) // less than, x < y { return ( NEQ(x, y) && (x < y) ); } inline bool GT(double x, double y) // greater than, x > y { return ( NEQ(x, y) && (x > y) ); } inline bool LEQ(double x, double y) // less equal, x <= y { return ( EQ(x, y) || (x < y) ); } inline bool GEQ(double x, double y) // greater equal, x >= y { return ( EQ(x, y) || (x > y) ); } // 输出浮点数前，防止输出-0.00调用该函数进行修正！ inline double FIX(double x) { return (fabs(x) < EPS) ? 0 : x; } /*------------------二维矢量运算重载区---------------------*/ bool operator==(Point p1, Point p2) { return ( EQ(p1.x, p2.x) && EQ(p1.y, p2.y) ); } bool operator!=(Point p1, Point p2) { return ( NEQ(p1.x, p2.x) || NEQ(p1.y, p2.y) ); } bool operator<(Point p1, Point p2) { if (NEQ(p1.x, p2.x)) { return (p1.x < p2.x); } else { return (p1.y < p2.y); } } Point operator+(Point p1, Point p2) { return Point(p1.x + p2.x, p1.y + p2.y); } Point operator-(Point p1, Point p2) { return Point(p1.x - p2.x, p1.y - p2.y); } double operator*(Point p1, Point p2) // 计算叉乘 p1 × p2 { return (p1.x * p2.y - p2.x * p1.y); } double operator&(Point p1, Point p2) { // 计算点积 p1·p2 return (p1.x * p2.x + p1.y * p2.y); } double Norm(Point p) // 计算矢量p的模 { return sqrt(p.x * p.x + p.y * p.y); } /*-------------------基本函数区------------------*/ //得到两点之间的距离 double Dis(Point p1,Point p2) { return sqrt( (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) ); } //求二维平面上点到直线的距离 double Dis(Point p, Line L) { return ( fabs((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1) ); } //得到两点之间距离的平方，为减少误差用 double Dis2(Point p1,Point p2) { return (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y); } //返回A关于B的对称点C，即(A+C)/2=b Point SymmetryPoint(Point A,Point B) { Point C; C.x = 2*B.x-A.x; C.y = 2*B.y-A.y; return C; } // 把矢量p旋转角度angle (弧度表示) // angle > 0表示逆时针旋转 // angle < 0表示顺时针旋转 Point Rotate(Point p, double angle) { Point result; result.x = p.x * cos(angle) - p.y * sin(angle); result.y = p.x * sin(angle) + p.y * cos(angle); return result; } //得到向量p与x正半轴的夹角[0,2PI) double GetAngle(Point p) { double tmp=atan2(p.y,p.x); if(tmp<0) tmp=2*PI+tmp; return tmp; } //得到两个向量之间的夹角[0,PI] //若p1按顺时针转到p2的角小于PI(p1在p2的逆时针方向)，则返回正数.否则返回负数. double GetAngle(Point p1,Point p2) { double tmp = GetAngle(p1) - GetAngle(p2); if( GT( tmp,PI) ) return -(2*PI-tmp); if( LEQ( tmp,-PI) ) return (tmp+2*PI); return tmp; } // 判断二维平面上点是否在线段上 // 输入：任意点p，和任意直线L // 输出：p在线段上返回1，否则返回0 bool OnSeg(Point p, Line L) { return ( ZERO( (L.p1 - p) * (L.p2 - p) ) && LEQ((p.x - L.p1.x)*(p.x - L.p2.x), 0) && LEQ((p.y - L.p1.y)*(p.y - L.p2.y), 0) ); } // 判断二维平面上点p是否在直线L上,在线段上返回1,否则返回0 bool OnLine(Point p, Line L) { return ZERO( (p - L.p1) * (L.p2 - L.p1) ); } bool OnCir(Point p,Circle cir) { return EQ( (p.x-cir.c.x)*(p.x-cir.c.x)+(p.y-cir.c.y)*(p.y-cir.c.y),cir.r*cir.r ); } //得到点p到直线L的距离，并返回p到直直线L的最近点rep double PointToLine(Point p,Line L,Point& rep) { if(L.p1==L.p2) { rep=L.p1; return Dis(p,L.p1); } Point a,b; a = L.p2-L.p1; b = p-L.p1; double dis12 = Dis(L.p1,L.p2); double dis = ( fabs(a*b) )/dis12; double k = (a&b)/(Norm(a)*dis12) ; rep.x = L.p1.x + k*(L.p2.x-L.p1.x); rep.y = L.p1.y + k*(L.p2.y-L.p1.y); return dis; } //得到点P到线段L的距离，并放回p到线段L的最近点rep double PointToSeg(Point P, Line L,Point& rep) { if(L.p1 == L.p2) { rep = L.p1; return Dis(rep,P);//如果线段是一个点，返回这个点。 } Point result; double a, b, t; a = L.p2.x - L.p1.x; b = L.p2.y - L.p1.y; t = ( (P.x - L.p1.x) * a + (P.y - L.p1.y) * b ) / (a * a + b * b);//线段上的投影 if ( GEQ(t, 0) && LEQ(t, 1) ) { result.x = L.p1.x + a * t;//值得学习的由比例求坐标的方法。 result.y = L.p1.y + b * t; } else { if ( Norm(P - L.p1) < Norm(P - L.p2) ) { result = L.p1; } else { result = L.p2; } } return Dis(result, P); } //返回点A关于直线L的对称点 Point SymmetryPonitToLine(Point A,Line L) { Point B; PointToLine(A, L, B);//A在L上的投影 return SymmetryPoint(A, B); } /*-------------------几何题面积计算（注意正负!）----------------------*/ // 根据三个顶点坐标计算三角形面积 // 面积的正负按照右手旋规则确定,向量AB->向量AC double Area(Point A, Point B, Point C) { return ((B-A)*(C-A) / 2.0); } // 根据三条边长计算三角形面积 double Area(double a, double b, double c) { double s = (a + b + c) / 2.0; return sqrt(s * (s - a) * (s - b) * (s - c)); } //求圆的面积 double Area(Circle C) { return PI * C.r * C.r; } // 计算多边形面积,复杂度:O(顶点数) // 面积的正负按照右手旋规则确定，顺时针为负 double Area(Polygon _poly) { int nsize=_poly.size(); double area=0; for(int i=0;i<nsize;i++) { area += _poly[i]*_poly[(i+1)%nsize]; } return area/2.0; } //求两条直线之间的关系（二维） //输入：两条不为点的直线 //输出：相交返回XIANGJIAO和交点p,平行返回PINGXING,共线返回GONGXIAN int LineAndLine(Line L1,Line L2,Point &p) { Point px,py; px = L1.p1 - L1.p2; py = L2.p1 - L2.p2; if( EQ(px*py,0) )//平行或者共线 { if( ZERO( (L2.p1-L1.p1)*py ) ) //共线 { return GONGXIAN; } return PINXING; } double xa,xb,xc,ya,yb,yc; xa=(L1.p2.y-L1.p1.y); xb=(L1.p1.x-L1.p2.x); xc=(L1.p1.y*L1.p2.x-L1.p1.x*L1.p2.y); ya=(L2.p2.y-L2.p1.y); yb=(L2.p1.x-L2.p2.x); yc=(L2.p1.y*L2.p2.x-L2.p1.x*L2.p2.y); p.y = (xa*yc-xc*ya)/(xb*ya-xa*yb); p.x = (xb*yc-xc*yb)/(xa*yb-xb*ya); return XIANGJIAO; } //判断两条线段是否相交，相交返回1 bool SegAndSeg(Line L1,Line L2) { return ( GEQ( max(L1.p1.x, L1.p2.x), min(L2.p1.x, L2.p2.x) ) && GEQ( max(L2.p1.x, L2.p2.x), min(L1.p1.x, L1.p2.x) ) && GEQ( max(L1.p1.y, L1.p2.y), min(L2.p1.y, L2.p2.y) ) && GEQ( max(L2.p1.y, L2.p2.y), min(L1.p1.y, L1.p2.y) ) && LEQ( ((L2.p1 - L1.p1) * (L1.p2 - L1.p1)) * ((L2.p2 - L1.p1) * (L1.p2 - L1.p1)), 0 ) && LEQ( ((L1.p1 - L2.p1) * (L2.p2 - L2.p1)) * ((L1.p2 - L2.p1) * (L2.p2 - L2.p1)), 0 ) ); } //求两条线段交点(二维) //输入：两条不为点的直线 //输出：相交返回XIANGJIAO和交点p，相离返回XIANGLI，重叠返回CHONGDIE int SegAndSeg(Line L1,Line L2,Point &p) { double signx,signy; //跨立实验 if( LEQ(signx=( ((L1.p2-L1.p1)*(L1.p1-L2.p1))*((L1.p2-L1.p1)*(L1.p1-L2.p2)) ),0) && LEQ(signy=( ((L2.p2-L2.p1)*(L2.p1-L1.p1))*((L2.p2-L2.p1)*(L2.p1-L1.p2)) ),0) ) { if( ZERO(signx) && ZERO(signy) ) { //线段共线 signx = min( max(L1.p1.x,L1.p2.x),max(L2.p1.x,L2.p2.x) )- max( min(L1.p1.x,L1.p2.x),min(L2.p1.x,L2.p2.x) ); signy = min( max(L1.p1.y,L1.p2.y),max(L2.p1.y,L2.p2.y) )- max( min(L1.p1.y,L1.p2.y),min(L2.p1.y,L2.p2.y) ); if( ZERO(signx) && ZERO(signy) ) //说明共线，且相交一点 { if(L1.p1==L2.p1||L1.p1==L2.p2) p=L1.p1; if(L1.p2==L2.p1||L1.p2==L2.p2) p=L1.p2; return XIANGJIAO; } else if( GEQ(signx, 0) && GEQ(signy, 0) ) { return CHONGDIE;//重叠 } else { return XIANGLI;//相离 } } return LineAndLine(L1, L2, p);//转化为直线相交 } return XIANGLI;//相离 } // 判断点p是否在简单多边形poly内, 多边形可以是凸的或凹的 // poly的顶点数目要大于等于3 // 返回值为： // INSIDE -- 点在poly内 // BORDER -- 点在poly边界上 // OUTSIDE -- 点在poly外 int InPolygon(const Polygon poly, Point p) { int i, n, count; Line ray, side; n = poly.size(); count = 0; ray.p1 = p; ray.p2.y = p.y; ray.p2.x = - INF;// 设定一个极大值 for (i = 0; i < n; i++) { side.p1 = poly[i]; side.p2 = poly[(i+1)%n]; if( OnSeg(p, side) ) { return BORDER; } // 如果side平行x轴则不作考虑 if ( EQ(side.p1.y, side.p2.y) ) { continue; } if (OnSeg(side.p1, ray)) { if ( GT(side.p1.y, side.p2.y) ) count++; } else if (OnSeg(side.p2, ray)) { if ( GT(side.p2.y, side.p1.y) ) count++; } else if ( SegAndSeg(ray, side) ) { count++; } } return ((count % 2 == 1) ? INSIDE : OUTSIDE); } //得到直线与圆的交点 int LineToCir(Line L,Circle R,Point p[2]) { if(L.p1 == L.p2)//当直线为1个点时 { if( EQ( Dis(L.p1, R.c),R.r ) ) { p[0]=L.p1; return 1; } else return 0;//相离 } Point tp;//表示圆心在直线L上的投影。 double dis=PointToLine(R.c, L,tp ); if( LT(R.r, dis) )//相离 { return 0; } if( EQ(dis,R.r) )//相切 { p[0]=tp; return 1; } double len=sqrt(R.r*R.r-dis*dis); Point onep=L.p2-L.p1; double _t = len/Norm(onep); p[0].x =tp.x + onep.x*_t; p[0].y =tp.y + onep.y*_t; onep=L.p1-L.p2; p[1].x =tp.x + onep.x*_t; p[1].y =tp.y + onep.y*_t; return 2; } //得到三角形外接圆 //注意：A,B,C三点不能共线 Circle OutCircle(Point A,Point B,Point C) { Circle tmp; double a, b, c, c1, c2; double xA, yA, xB, yB, xC, yC; a = Dis(A, B); b = Dis(B, C); c = Dis(C, A); //根据 S = a * b * c / R / 4;求半径 R tmp.r = (a*b*c)/( fabs(Area(A,B,C)) *4.0); xA = A.x; yA = A.y; xB = B.x; yB = B.y; xC = C.x; yC = C.y; c1 = (xA*xA+yA*yA - xB*xB-yB*yB) / 2; c2 = (xA*xA+yA*yA - xC*xC-yC*yC) / 2; tmp.c.x = (c1*(yA - yC)-c2*(yA - yB)) / ((xA - xB)*(yA - yC)-(xA - xC)*(yA - yB)); tmp.c.y = (c1*(xA - xC)-c2*(xA - xB)) / ((yA - yB)*(xA - xC)-(yA - yC)*(xA - xB)); return tmp; } //得到三角形内切圆 //注意：A,B,C三点不能共线 Circle InCircle(Point A,Point B,Point C) { Circle rec; double a=Dis(B,C); double b=Dis(A,C); double c=Dis(A,B); rec.c.x = (a*A.x+b*B.x+c*C.x)/(a+b+c); rec.c.y = (a*A.y+b*B.y+c*C.y)/(a+b+c); rec.r = 2*fabs( Area(A,B,C) )/(a+b+c); return rec; } //得到两圆的面积并 double CirArea(Circle _c1,Circle _c2) { if(_c2.r<_c1.r) swap(_c1,_c2); //保证_c2的半径大 double d12 = Dis(_c1.c,_c2.c); if( LEQ(_c1.r+_c2.r,d12) )//相离 { return 0; } if( LEQ(d12+_c1.r, _c2.r) )//包含 { return Area(_c1); } //相交 double _area=0; _area -= 2.0*Area(_c1.r,_c2.r,d12); double ang1 = acos( (d12*d12+_c1.r*_c1.r-_c2.r*_c2.r) / (2*d12*_c1.r) ); double ang2 = acos( (d12*d12+_c2.r*_c2.r-_c1.r*_c1.r) / (2*d12*_c2.r) ); _area += ang1*_c1.r*_c1.r+ang2*_c2.r*_c2.r; return _area; } //得到两个圆的交点p[2] //返回值为交点数，-1为两圆重叠。 int CirAndCir(Circle _c1,Circle _c2,Point p[2]) { if(_c2.r < _c1.r) swap(_c1,_c2); //保证_c2的半径大 double d12 = Dis(_c1.c,_c2.c); if( LT(_c1.r+_c2.r,d12) )//相离 { return 0; } if( LT(d12+_c1.r, _c2.r) )//包含 { return 0; } if(_c1.c == _c2.c)//两个圆重叠 { return -1; } Point u,v; double t; double r1=_c1.r,r2=_c2.r; t=( 1+(r1*r1-r2*r2)/(Dis(_c1.c,_c2.c)*Dis(_c1.c,_c2.c)) ) /2; u.x = _c1.c.x + (_c2.c.x-_c1.c.x)*t; u.y = _c1.c.y + (_c2.c.y-_c1.c.y)*t; v.x = u.x + _c1.c.y - _c2.c.y; v.y = u.y - _c1.c.x + _c2.c.x; Line _l(u,v); return LineToCir(_l,_c1,p); } //求凸包Graham-Scan（GS）算法，复杂度nlog(n),内附两种模式，最小点集凸包，与最大点集凸包。 //调用GetConvex_GS(Point ps[],int pn,Point cx[],int &cxn) //其中ps是输入的点集，pn为ps的大小，cx为返回的凸包集，cxn表示凸包的大小。 //注意：pn必须>=3,下标都是从0开始，返回的凸包集为逆时针。且如果点带标号，一样处理。 int GScmp(Point a,Point b) { double _tmp=a*b; Point zero; zero.x=0; zero.y=0; if( ZERO(_tmp) ) { return Dis2(a,zero)>Dis2(b,zero); } return _tmp > 0; } void GetConvex_GS(Point ps[],int pn,Point cx[],int &cxn) { Polygon pg; pg.clear(); //先找出最下面，如果有相同的找最靠左，也就找y轴最小的，然后y相同时选x最小。 Point minp=ps[0]; for(int i=1;i<pn;i++) { if( ps[i].y < minp.y ) { minp = ps[i]; } else if(EQ(ps[i].y, minp.y) && ps[i].x < minp.x) { minp = ps[i]; } } for(int i=0;i<pn;i++) { ps[i].x -= minp.x; ps[i].y -= minp.y; pg.push_back( ps[i] ); } //排序 sort(pg.begin(),pg.end(),GScmp); //在这一步，除去与minp共线的点。 long int pgsize=pg.size(); int pgcnt=1; for(int i=1;i<pgsize;i++) { if( !ZERO(pg[i]*pg[pgcnt-1]) ) { pg[ pgcnt++ ] = pg[i]; } } cxn = 0; cx[ cxn++ ] = minp; if( !(ZERO(pg[0].x) && ZERO(pg[0].y)) )//不为一点时 { pg[0].x += minp.x; pg[0].y += minp.y; cx[ cxn++ ] = pg[0];//必须保存id for(int i=1;i < pgcnt;i++) { pg[i].x += minp.x; pg[i].y += minp.y; double _tmp = (cx[cxn-1]-cx[cxn-2])*(pg[i]-cx[cxn-1]); // 无法处理重复点！ // nice ! while( LEQ(_tmp, 0) ) { cxn--; //if(cxn==1) break;//只剩下一个的时候，退出 _tmp = (cx[cxn-1]-cx[cxn-2])*(pg[i]-cx[cxn-1]); } cx[ cxn++ ] = pg[i]; } } //已经找到了，最小凸包集,且为逆时针排序。 /* //接下来步骤为找出所有在凸包边界上的点，并按逆时针排序。 sort(ps,ps+pn,GScmp); int cxn1=0; int tj=0; for(int i=0;i<pn;i++) { ps[i].x += minp.x; ps[i].y += minp.y; Line l(cx[tj],cx[(tj+1)%cxn]); if( OnSeg(ps[i], l) ) { ps[ cxn1++ ] = ps[i]; continue; } if( tj+1<cxn && GEQ( (cx[tj+1]-minp)*(ps[i]-minp),0 ) ) { tj++; } l.p1 = cx[ tj ]; l.p2 = cx[ (tj+1)%cxn ]; if(OnSeg(ps[i], l)) ps[cxn1++]=ps[i]; } cxn=cxn1; for(int i=0;i<cxn1;i++) cx[i]=ps[i]; //接下来的步骤是，将凸包点集调整为逆时针，其实只需要调整斜率最大的一条边。 tj=1; Line l(minp,cx[0]); while(tj<cxn && OnSeg(cx[tj], l)) tj++; tj--; for(int i=0;i<=tj;i++) cx[i] = ps[tj-i]; //////////////////////////////////// */ } Point BaryCenter(Point p[], int n){ Point ret = Point(0,0); double area = 0; for(int i = 1; i < n-1; i++){ double area2 =fabs( Area(p[0], p[i], p[i+1]) ); ret.x += (p[0].x+p[i].x + p[i+1].x) / 3 * area2; ret.y += (p[0].y+p[i].y + p[i+1].y) / 3 * area2; area += area2; } ret.x /= area; ret.y /= area; return ret; } /*---------------------代码区---------------------------*/