南京华夏天成建设有限公司网站,自己做简单网站,网站开发多久,怎么授权小说做游戏网站文章目录前言精度点/向量相关表示向量基本运算角度相关向量夹角旋转直线/线段相关表示点与线求点到直线垂足求点关于直线的对称点点与直线的位置关系点与直线的距离线与线直线与直线的位置关系共线与垂直判断线段与线段是否相交求直线与直线的交点角平分线中垂线多边形表示求多…
文章目录前言精度点/向量相关表示向量基本运算角度相关向量夹角旋转直线/线段相关表示点与线求点到直线垂足求点关于直线的对称点点与直线的位置关系点与直线的距离线与线直线与直线的位置关系共线与垂直判断线段与线段是否相交求直线与直线的交点角平分线中垂线多边形表示求多边形面积判断多边形凹凸性判定点在多边形内判定点在凸包内圆表示位置关系点与圆的位置关系线与圆的位置关系圆与圆的位置关系三角形和圆三角形内切圆三角形外接圆交点相关圆与直线交点圆与圆的交点切点/线相关点到圆的切点圆与圆公切线的切点面积相关圆、扇形和弓形圆与多边形面积交圆与圆的面积交二维凸包定义维护代码旋转卡壳代码半平面交半平面的表示流程细节代码其他皮克定理距离曼哈顿距离和切比雪夫距离的转化source前言
强烈推荐AOJ的入门几何题库。 本文是一篇入门全解需要的前置知识只有基础的高中的解析几何和向量相关。 计算几何的一些原则
精度即生命。炸精度是一个糟糕的计算几何实现经常出现的问题就像我一开始尝试自己造轮子常常遇到的为了避免过大的精度误差应该尽量减少正反三角函数、除法、根号的使用次数但似乎根号函数的精度还是不错的。实在不行long double 或许能救我们一命但是会变慢简洁即美德。这一条是建立在不过分影响精度的前提之上的。很多时候对于某个问题存在多种情况但对于一种较简单情况给出的某种解法归纳发现在其他情况下也是正确的就可以大大减少分类讨论比如线段求交问题。此外代码合理的模块化也能使代码变得更加简洁。
所以尽管我们自己往往也能乱搞出对于某个问题的一种做法但是为了精度更高、更加简洁的实现学习计算几何是必要的。 后面的函数里可能会直接调用前面讲完的一些函数这也是代码模块化简化代码的一个体现。 注意任何时候用除法都要考虑考虑会不会除0 精度
计算几何常用浮点数由于浮点数本身难以避免、只能减小的精度误差常常不使用通常的、、 运算。 设 epsepseps 为一个极小量。通常在 [10−12,10−8][10^{-12},10^{-8}][10−12,10−8] 之间 那么三种运算符再浮点运算下就可以表示为 ab→abs(a−b)epsab \to\operatorname{abs}(a-b)epsab→abs(a−b)eps ab→ab−epsab \to ab-epsab→ab−eps ab→abepsab \to abepsab→abeps 这不是死记硬背的。其意义就是假设当 ababab 时表达式不应该成立那么类似的 和 应该在 ababab 时成立所以应该把对应的 epsepseps 之前的符号反过来。
点/向量相关
表示
由于在计算几何中点和向量常常是等价的所以归为一类。 只需要记录横纵坐标即可。
//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
V ans[10];//declared for other functions
int tot;
inline void input(V a){a.xread();a.yread();}
void print(const V a,int op1){printf(%.10lf %.10lf,a.x,a.y);putchar(op?10:32);}
//op:endl or space向量基本运算
基本的加、减、乘、除、点积、叉积、模长、中点、垂直向量、向量单位化、三角形面积。 拓展 高维点积定义不变但在投影方面几何意义只对二维和三维空间有效。 高维叉积叉积在k维下的定义是三行k列的行列式第一行为各维度的单位向量二三行为两个点各自在对应维度的坐标。 注意
严谨来说向量叉积是一个向量而不是一个数它的模长等于 x1y2−x2y1x_1y_2-x_2y_1x1y2−x2y1方向垂直与纸面遵循右手定则这也是叉积可以表示负面积的原因。这里的垂直向量没有关注方向。
//basic operation
inline V operator (const V a,const V b){return (V){a.xb.x,a.yb.y};}
inline V operator - (const V a,const V b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double x,const V a){return (V){a.x*x,a.y*x};}
inline V operator * (const V a,const double x){return (V){a.x*x,a.y*x};}
inline V operator / (const V a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator (const V a,const V b){return abs(a.x-b.x)epsabs(a.y-b.y)eps;}
inline bool operator ! (const V a,const V b){return !(ab);}
inline double operator * (const V a,const V b){return a.x*b.xa.y*b.y;}
inline double operator ^ (const V a,const V b){return a.x*b.y-a.y*b.x;}
inline double len(const V a){return sqrt(a.x*a.xa.y*a.y);}
inline V mid(const V a,const V b){return (V){(a.xb.x)/2,(a.yb.y)/2};}
inline V cui(const V a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V a){return a/len(a);}
inline double tri_S(const V a,const V b,const V c){return abs((b-a)^(c-a))/2;}//always be non-negative
inline bool operator (const V a,const V b){return a.xb.x-eps||(abs(a.x-b.x)epsa.yb.y-eps);
}角度相关
向量夹角
利用点积易得。 注意这个角的值域是0到π。
inline double angle(const V a,const V b) { return acos(a * b / len(a) / len(b)); }
inline bool dun(const V a,const V b,const V c){return ((b-a)*(c-a))-eps;}//angle:BAC
inline bool rui(const V a,const V b,const V c){return ((b-a)*(c-a))eps;}
inline bool zhi(const V a,const V b,const V c){return abs(((b-a)*(c-a)))eps;}旋转
若原来向量为 (x,y)(x,y)(x,y)旋转角为 θ\thetaθ那么旋转后的向量就是 (xcosθ−ysinθ,xsinθycosθ)(x\cos\theta-y\sin\theta,x\sin\thetay\cos\theta)(xcosθ−ysinθ,xsinθycosθ)。 可以用虚数乘法(xyi)×(cosθ,sinθi)(xyi)\times (\cos\theta,\sin\theta \space i)(xyi)×(cosθ,sinθ i)结合几何意义简单证明。
inline V rotate(const V o,double t){double ssin(t),ccos(t);return (V){o.x*c-o.y*s,o.x*so.y*c};
}直线/线段相关
表示
直线常用方向向量加上线上任意一点表示。 线段常用两个端点表示。 把方向向量转化为单位向量会比较方便下面直线的方向向量默认为单位向量。 下文如非特别说明在大多数时候不在特别区分线段和直线两者的算法是大同小异的。
//definition
struct line{V d,a,b;
};
inline line trans(double a,double b,double c,double d){//given coordinatesV dd(c-a,d-b),x(a,b),y(c,d);dddd/len(dd);return (line){dd,x,y};
}
inline line trans(const V a,const V b){//given pointsV dd(b.x-a.x,b.y-a.y);dddd/len(dd);return (line){dd,a,b};
}
inline void input(line l){double aread(),bread(),cread(),dread();ltrans(a,b,c,d);return;
}点与线
求点到直线垂足
把线上的代表点移动投影的距离。 这里的投影是带符号的符号体现移动的方向。
inline V cui(const line l,const V o){//directedreturn l.a((o-l.a)*l.d)*l.d;
}求点关于直线的对称点
求出垂足后用中点公式即可。
inline V duichen(const line l,const V o){V ucui(l,o);return (V){2*u.x-o.x,2*u.y-o.y};
}点与直线的位置关系
共线可以通过叉积等于0来判断。 不共线时可以通过叉积的正负判断是在直线的哪一侧。 共线时如果是线段或射线还可以利用点积的正负判断与端点的位置关系。 对于线段而言可以直接通过模长判断是否在线段上
inline bool on_line(const line l,const V o){return abs((l.d^(o-l.a)))eps;}
inline bool on_seg(const line l,const V o){return abs(len(o-l.a)len(o-l.b)-len(l.a-l.b))eps;
}
inline int pos(const line l,const V o){if(!on_line(l,o)){if((l.d^(o-l.a))eps) return 1;//counter clockwiseelse return 2;//clockwise}else{if(((o-l.a)*(o-l.b))eps) return 5;//on segmentelse if(((o-l.a)*l.d)-eps) return 3;//online backelse return 4;//online front}
}
点与直线的距离
直线时可以利用用叉积等面积法求出距离。 线段时若垂足在线段上还是点到直线的距离否则是到两端点距离的最小值。 垂足是否在线段上可以利用是否有钝角来判断。
inline double dis(const line l,const V o,int op0){//op0:dis on line,op1:dis on segmentif(op(dun(l.a,o,l.b)||dun(l.b,o,l.a))) return min(len(o-l.a),len(o-l.b));else return abs((o-l.a)^(o-l.b))/len(l.a-l.b);
}线与线
直线与直线的位置关系
共线与垂直
如果方向向量叉积是0则共线注意还要分平行和重合如果点积是0则垂直。
inline bool gongxian(const line a,const line b){return abs(a.d^b.d)eps;}
inline bool cuizhi(const line a,const line b){return abs(a.d*b.d)eps;}判断线段与线段是否相交
需要两个步骤快速排斥实验和跨立实验。
快速排斥实验定义一条线段所在的矩形为边与坐标轴平行、能完全包含该线段的最小的矩形如果两条线段所在的矩形没有公共部分显然不可能有交点否则进入下一步跨立实验。跨立实验如果l1的两端点在l2两侧同时l2的两端点也在l1两侧也可以在线段上则得出两线段有公共交点否则则没有公共交点。
实现上快速排斥实验可以用坐标 min、max 的比较轻松解决跨立实验则可以用前面讲的叉积求点与线的关系来实现注意叉积为0端点在线段上是符合条件的。 直观上感觉似乎跨立实验已经很充要了但是由于判断的时候叉积等于0认为合法会把没有公共部分的共线线段判为有交点所以快速排斥实验是必要的。也可以用特判共线来取代 inline bool jiao(const line u,const line v){if(min(u.a.x,u.b.x)max(v.a.x,v.b.x)eps||max(u.a.x,u.b.x)min(v.a.x,v.b.x)-eps||min(u.a.y,u.b.y)max(v.a.y,v.b.y)eps||max(u.a.y,u.b.y)min(v.a.y,v.b.y)-eps) return false;//rapid rejection testreturn ((u.a-v.a)^v.d)*((u.b-v.a)^v.d)eps((v.a-u.a)^u.d)*((v.b-u.a)^u.d)eps;//straddle test
}求直线与直线的交点
首先应保证两直线不共线。 设两直线为 l1,l2l1,l2l1,l2线上各有一点 x1,x2x1,x2x1,x2方向向量为 d1,d2d1,d2d1,d2。 假设交点 px1kd1px1kd1px1kd1就有 (p−x2)×d20(p-x2)\times d20(p−x2)×d20。 叉积是满足分配率的拆开化简得到 k(x2−x1)×d2d1×d2k\frac{(x2-x1)\times d2}{d1\times d2}kd1×d2(x2−x1)×d2 带回去即可。
inline V jiaodian(line u,line v){double k((v.a-u.a)^v.d)/(u.d^v.d);return u.a(k*u.d);
}角平分线
求出两边的单位向量d1,d2d1d2就是角平分线的方向向量。
inline line pingfen(const V a,const V b,const V c){//angle:BACV d1(b-a)/len(b-a),d2(c-a)/len(c-a),d(d1d2)/len(d1d2);return (line){d,a,ad};
}中垂线
找到中点和与原线段垂直的方向向量即可。
inline line zhongcui(const V a,const V b){V xmid(a,b),ddanwei(cui(b-a));return (line){d,x,xd};
}多边形
表示
用一个V数组顺时针或逆时针存储多边形的所有顶点即可。 由于封装的话老打点会比较烦就直接传数组和顶点数n了。
求多边形面积
取任意一个点 OOO则多边形面积为 S12∑i1nOAi→×OA→i%n1S\frac{1}{2}\sum_{i1}^n\overrightarrow{OA_{i}}\times \overrightarrow{OA}_{i\%n1}S21i1∑nOAi×OAi%n1 为了方便OOO 直接取原点即可。 时间复杂度 O(n)O(n)O(n)。
inline double S(const V *a,const int n){double res(0);for(int i1;in;i) res(a[i]^a[i%n1]);return res/2;
}判断多边形凹凸性
取相邻的三个点 (a,b,c)(a,b,c)(a,b,c)作为一个三元组。 判断多边形凹凸的充要条件对于多边形所有的 nnn 个这样三元组ccc 和 ababab 的位置关系指顺时针或逆时针必然全部都是相同的。 叉积判断符号即可。 注意如果没有保证相邻三个点不共线的话需要特判。 时间复杂度 O(n)O(n)O(n)。
bool is_convex(const V *a,const int n){a[0]a[n];a[n1]a[1];int op0;for(int i1;in;i){double o(a[i1]-a[i])^(a[i]-a[i-1]);if(abs(o)eps) continue;//three neighbouring points on the same lineint nopo0?1:-1;if(!op) opnop;else if(op!nop) return false;}return true;
}判定点在多边形内
一个看起来很好写但写起来很恶心的东西。 推荐使用点数判别法。 向左水平引出一条射线计算其与多边形产生多少个交点奇内偶外然后还有亿点点细节
如果当前边经过了该点直接返回相应结果由于可能出现射线经过多边形顶点导致一个点被计算两次的情况强制令这个点在下面的边被统计。实现上当 max(y1,y2)y0\max(y1,y2)y0max(y1,y2)y0 时统计该次答案而当 min(y1,y2)y0\min(y1,y2)y0min(y1,y2)y0 时不统计该次答案。出现水平边时直接忽略即可结合第二条的处理原则就依然是正确的。
int in_poly(const V *a,const int n,const V o){int res(0);for(int i1;in;i){V ua[i],va[i1];if(on_seg(trans(u,v),o)) return 1;//on the edgeif(abs(u.y-v.y)eps) continue;//ignore horizontalif(max(u.y,v.y)o.y-eps) continue;//equal will not continueif(min(u.y,v.y)o.y-eps) continue;//equal will continuedouble xu.x(o.y-u.y)/(v.y-u.y)*(v.x-u.x);if(xo.x) res^1;}return res?2:0;//odd:in even:out
}判定点在凸包内
首先需要保证 a1(0,0)a_1(0,0)a1(0,0)。如果不满足整体平移即可 lowerbound 找到极角夹在询问点两侧的相邻点对判断是否在二者连成的边内部即可。 注意需要对于在a1,ana_1,a_na1,an 外侧的点需要特判。
bool cmp2(V a,V b){return (a^b)0;
}
bool in(const V *c,const V o){if(((c[n]^o)eps)||((o^c[2])eps)) return 0;int pllower_bound(c1,c1n,o,cmp2)-c-1;return ((c[pl1]-c[pl])^(o-c[pl]))-eps;
}圆
表示
记录圆心和半径即可表示。
struct cir{V o;double r;cir(double a0,double b0,double c0):o(a,b),r(c){}
};
inline void input(cir cc){cc.o.xread();cc.o.yread();cc.rread();}位置关系
点与圆的位置关系
判断距离即可。
inline bool incir(const cir c,const V p){return len(p-c.o)c.r-eps;}
inline bool oncir(const cir c,const V p){return (len(p-c.o)-c.r)eps;}
inline bool outcir(const cir c,const V p){return len(p-c.o)c.reps;}线与圆的位置关系
按照线圆距和半径的关系判断即可。
inline double dis(const cir c,const line l){return dis(l,c.o);}
inline int pos(const cir c,const line l){double ddis(c,l);if(dc.r-eps) return 1;//intersectelse if(abs(d-c.r)eps) return 2;//tangentelse if(dc.reps) return 3;//disjoint
}圆与圆的位置关系
分为外离、外切、相交、内切、内含五种。 公切线的数量对应为4、3、2、1、0。 通过圆心距和半径之间的关系判断即可。
inline double dis(const cir c1,const cir c2){return len(c1.o-c2.o);}
inline int pos(const cir c1,const cir c2){double ddis(c1,c2);if(dc1.rc2.reps) return 4;//do not crosselse if(abs(d-c1.r-c2.r)eps) return 3;//circumscribedelse if(max(c1.r,c2.r)min(c1.r,c2.r)d -eps) return 2;//intersectelse if(abs(max(c1.r,c2.r)-min(c1.r,c2.r)-d)eps) return 1;//inscribedelse return 0;//include
}三角形和圆
三角形内切圆
圆心就是两条平分线的交点即可。 半径可以用公式 2Sabc\dfrac{2S}{abc}abc2S也可以调用点到直线距离。
inline cir incir(const V a,const V b,const V c){V x(jiaodian(pingfen(a,b,c),pingfen(b,a,c)));return cir(x,dis(trans(a,b),x));
}三角形外接圆
我使用的是直接求中垂线的交点。 但是这么做有点炸精度…开longdouble才能过。
inline cir outcir(const V a,const V b,const V c){V x(jiaodian(zhongcui(a,b),zhongcui(a,c)));return cir(x,len(x-a));
}
inline void input(cir cc){in(cc.o);cc.rread();}法二 设三边边长为 a,b,ca,b,ca,b,c则 Rabc(abc)(ab−c)(ac−b)(bc−a)R\dfrac{abc}{\sqrt{(abc)(ab-c)(ac-b)(bc-a)}}R(abc)(ab−c)(ac−b)(bc−a)abc。
交点相关
圆与直线交点
求出圆心到直线距离后勾股定理即可。
inline double calc(double xie,double zhi){return sqrt(xie*xie-zhi*zhi);}//Pythagorean Theorem
inline void cross_line_cir(const cir c,const line l){tot0;V xchui(l,c.o);double ddlen(x-c.o);if(c.rdd-eps) return;//none cross pointif(abs(c.r-dd)eps){//one cross pointans[tot]x;return;}double discalc(c.r,dd);V axdis*l.d,bx-dis*l.d;//two cross pointsans[tot]a;ans[tot]b;return;
}圆与圆的交点
首先得保证存在交点。 考虑两圆心和任意一个交点组成的三角形可以发现这个三角形的三边都是已知的。 我们就可以用余弦定理求出夹角然后转上去即可。
inline void cross_cir_cir(const cir c1,const cir c2){int oppos(c1,c2);if(op4||op0) return;//none cross pointV Lc2.o-c1.o;double ac2.r,bc1.r,clen(L);double tacos((b*bc*c-a*a)/(2*b*c));//find the angleV xc1.oc1.r*rotate((L)/len(L),t),yduichen(trans(c1.o,c2.o),x);if(xy) print(x,0),print(y,1);else print(y,0),print(x,1);
}切点/线相关
点到圆的切点
切点、圆心、给定点自然而然形成了一个直角三角形。 求出角后转一下即可。
inline void qiedian(const cir c,const V p){tot0;line Ltrans(c.o,p);double tacos(c.r/len(c.o-p));V ac.o(c.r*rotate(L.d,t)),bduichen(L,a);ans[tot]a;if(a!b) ans[tot]b;return;
}圆与圆公切线的切点
分为外侧切线和内测切线。 两者大致的思路类似都是试图求出切点圆心连线与两圆心连线的夹角然后转上去。 符号的正负很和谐不必特殊考虑。
inline void common_qie(const cir c1,const cir c2){tot0;int oppos(c1,c2);if(op){//outside tangent linesdouble dc1.r-c2.r,tacos(d/len(c1.o-c2.o));V uc1.o(c1.r*danwei(rotate(c2.o-c1.o,t))),vc1.o(c1.r*danwei(rotate(c2.o-c1.o,-t)));ans[tot]u;if(u!v) ans[tot]v; }if(op2){//inside tangent linesdouble dlen(c2.o-c1.o)/(c1.rc2.r)*c1.r,tacos(c1.r/d);V uc1.o(c1.r*danwei(rotate(c2.o-c1.o,t))),vc1.o(c1.r*danwei(rotate(c2.o-c1.o,-t)));ans[tot]u;if(u!v) ans[tot]v;}
}面积相关
圆、扇形和弓形
比较基础套初中公式即可。
inline double cir_S(const cir c){return pi*c.r*c.r;}
inline double shan(const cir c,const V a,const V b){//S of sectordouble tang(a-c.o,b-c.o);return t/2*c.r*c.r;
}
inline double gong(const cir c,const V a,const V b){//S of bowreturn shan(c,a,b)-tri_S(c.o,a,b);
}圆与多边形面积交
和求多边形面积类似的问题可以转化为求多边形 n 对相邻顶点与圆心组成的三角形与圆的有向面积交。 然后需要亿点点大力分讨…具体可见代码和注释不太难理解 一个实现的 trick 是先把有向面积的符号求出来后面求面积的绝对值就行了。
inline double O_tri(const cir c,V a,V b){
//directed S of Intersection between circle O and triangle OABif(on(trans(a,b),c.o)) return 0.0;int f(((a-c.o)^(b-c.o))0)?1:-1;//directionline ltrans(a,b);int f1incir(c,a),f2incir(c,b);if(f1f2) return f*tri_S(a,b,c.o);//both incircle:the S of triangleelse if(!f1!f2){//both outcircle:a sector(possibly minus a bow)V u(c.o(c.r*danwei(a-c.o))),v(c.o(c.r*danwei(b-c.o)));double Sshan(c,u,v); if(dis(l,c.o,1)c.r-eps){cross_line_cir(c,l);assert(tot2);S-gong(c,ans[1],ans[2]);}return f*S;}else{//one in and one out:a traingle and a sectorif(!f1) swap(a,b);double saSin(b-a,c.o-a),susa/c.r*len(c.o-a),Aasin(sa),Uasin(su);if((b-a)*(c.o-a)-eps) Api-A;double tpi-A-U,dissin(t)/sa*c.r; V ua(dis*danwei(b-a)),vc.r*danwei(b-c.o);return f*(tri_S(c.o,a,u)shan(c,u,v));}
}
double common_cir_poly(const V *a,const cir c,int n){double res(0);for(int i1;in;i) resO_tri(c,a[i],a[i1]);return res;
}圆与圆的面积交
首先把较为简单的外离和内含的情况判掉重点考虑相交。 先考虑比较简单的情形。交点与两圆心连线夹角都是锐角时相交部分是两个弓形利用余弦定理求出相交部分的圆心角然后用扇形面积减三角形面积即可。 然后发现如果变成了钝角三角形面积会变成一个负面积减完等于加上三角形面积的绝对值结果还是对的十分和谐不需要特判。
inline double common_cir_cir(const cir c1,const cir c2){int oppos(c1,c2);if(op3) return 0;//common area0else if(op1) return min(cir_S(c1),cir_S(c2));//completly includeelse{//partly includedouble Llen(c1.o-c2.o);double t12*acos((L*Lc1.r*c1.r-c2.r*c2.r)/(2*L*c1.r));double t22*acos((L*Lc2.r*c2.r-c1.r*c1.r)/(2*L*c2.r));return c1.r*c1.r*t1/2-c1.r*c1.r*sin(t1)/2c2.r*c2.r*t2/2-c2.r*c2.r*sin(t2)/2;}
}二维凸包
解析集合的第一课。 关键特征周长最小。此时一定是凸包。
定义 凸包在平面上能包含所有给定点的最小凸多边形叫做凸包。 性质凸包的周长是所有能包含给定点的多边形中最小的。
维护
凸包的主流求法有 Andrew 和 Graham两者的复杂度瓶颈都在于排序为 O(nlogn)O(n\log n)O(nlogn)。 这里介绍较为简单的 Graham 扫描法。
首先找出所有点中按照 X−YX-YX−Y 排序最小的点 OOO。 然后以 OOO 作为原点把所有其它点按照极角排序极角相同时按距离升序排序实测这里也可以降序但是绝对不能无序。 可以用快排重载 cmp 函数的方法实现利用叉积判断。
bool cmp(V a,V b){double d(a-p[1])^(b-p[1]);if(abs(d)eps) return d0;else return len(a-p[1])len(b-p[1]);
}排好序后开一个栈维护当前凸包中的点。 如果新点破坏了凸性则不断退栈。 最终栈内的元素就是凸包中的点。
是否破坏凸性可以用叉积判断如果新点和栈顶形成的向量向右拐了叉积小于0则说明破坏了凸性。 注意这里判断退栈的条件最好带等一方面共线时在边上的顶点没有什么意义另一方面当有重合点时不带等会导致程序错误。 void ConvexHull(V *p,int n){int top0;sort(p1,p1n);sort(p2,p1n,cmp);top0;for(int i1;in;i){while((top1((zhan[top]-zhan[top-1])^(p[i]-zhan[top]))0)) --top;zhan[top]p[i];}memcpy(p,zhan,sizeof(zhan));ntop;return;
}代码
P2742 [USACO5.1]圈奶牛Fencing the Cows /【模板】二维凸包
#includebits/stdc.h
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
inline ll read(){ll x(0),f(1);char cgetchar();while(!isdigit(c)){if(c-)f-1;cgetchar();}while(isdigit(c)){x(x1)(x3)c-0;cgetchar();}return x*f;
}//basic declare
//#define double long double
const double eps1e-10;
const double piacos(-1.0);//---about vectors (or points)//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
V ans[10];//declared for other functions
int tot;
inline void input(V a){scanf(%lf%lf,a.x,a.y);}
void print(const V a,int op1){printf(%.10lf %.10lf,a.x,a.y);putchar(op?10:32);}
//op:endl or space//basic operation
inline V operator (const V a,const V b){return (V){a.xb.x,a.yb.y};}
inline V operator - (const V a,const V b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double x,const V a){return (V){a.x*x,a.y*x};}
inline V operator * (const V a,const double x){return (V){a.x*x,a.y*x};}
inline V operator / (const V a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator (const V a,const V b){return abs(a.x-b.x)epsabs(a.y-b.y)eps;}
inline bool operator ! (const V a,const V b){return !(ab);}
inline double operator * (const V a,const V b){return a.x*b.xa.y*b.y;}
inline double operator ^ (const V a,const V b){return a.x*b.y-a.y*b.x;}
inline double len(const V a){return sqrt(a.x*a.xa.y*a.y);}
inline V mid(const V a,const V b){return (V){(a.xb.x)/2,(a.yb.y)/2};}
inline V chui(const V a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V a){return a/len(a);}
inline double tri_S(const V a,const V b,const V c){return abs((b-a)^(c-a))/2;}//always be non-negative
inline bool operator (const V a,const V b){return a.xb.x-eps||(abs(a.x-b.x)epsa.yb.y-eps);
}
inline double ang(const V a,const V b){return acos((a*b)/len(a)/len(b));}
inline V rotate(const V o,double t){//COUNTER_CLOCKWISEdouble ssin(t),ccos(t);return (V){o.x*c-o.y*s,o.x*so.y*c};
}const int N1e5100;
const int M505;
int n,m;
V p[N],zhan[N];
bool cmp(V a,V b){double d(a-p[1])^(b-p[1]);if(abs(d)eps) return d0;else return len(a-p[1])len(b-p[1]);
}
void ConvexHull(V *p,int n){int top0;sort(p1,p1n);sort(p2,p1n,cmp);top0;for(int i1;in;i){while((top1((zhan[top]-zhan[top-1])^(p[i]-zhan[top]))0)) --top;zhan[top]p[i];}memcpy(p,zhan,sizeof(zhan));ntop;return;
}
signed main(){
#ifndef ONLINE_JUDGE//freopen(a.in,r,stdin);//freopen(a.out,w,stdout);
#endifnread();for(int i1;in;i) input(p[i]);ConvexHull(p,n);zhan[n1]p[1];double res(0);for(int i1;in;i) reslen(zhan[i1]-zhan[i]);//print(zhan[i]);printf(%.2lf\n,res);return 0;
}
/*
3 5
0 -2
-5 3
0 -7
*/
旋转卡壳
前置知识凸包 个人感觉很像 two-pointers 算法。 能够在优秀的线性时间复杂度内完成总多求最值周长、面积…的神奇操作。
给出情境 给出平面内的 nnn 个点求出所有点中的最远点对。 n≤105n\le 10^5n≤105 首先有一个结论最远点对一定都是点集的凸包的顶点。 较为显然证明可以考虑把凸包内的点延长到凸包一条边上边两边的顶点一定有一个更优。
那么我们就转化成了求凸包上的最远点对这个问题也叫做凸包的直径问题。
给出一些定义 凸包的切线若一条直线过凸包上的一点或一边且整个凸包都在直线的同侧或在线上那么我们就称这条直线为凸包的切线。 对踵点如果经过凸包的两个顶点可以作两条平行的凸包的切线那么就称这两个点是一对对踵点。 不难发现最远点对一定是一对对踵点。 然而个人感觉旋转卡壳这个知识点完全不需要这个概念。
考虑换一个角度每次枚举边然后用到边距离最远的点和边的两端点的距离来更新答案。每次更新答案的点其实都是对踵点 显然最优答案一定会被枚举到。 不难发现如果我们逆时针枚举边最远点的位置也是在逆时针旋转。 那么我们利用类似 two-pointers 的思想就可以线性的求出答案。 问题得以解决。
实现的细节上我比较喜欢的方法是一开始先扫一遍暴力找到指针的起始位置而不是倍长野蛮或者每次移动指针都更新答案玄学。
代码
P1452 [USACO03FALL]Beauty Contest G /【模板】旋转卡壳
#includebits/stdc.h
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
inline ll read(){ll x(0),f(1);char cgetchar();while(!isdigit(c)){if(c-)f-1;cgetchar();}while(isdigit(c)){x(x1)(x3)c-0;cgetchar();}return x*f;
}//basic declare
//#define double long double
const double eps1e-10;
const double piacos(-1.0);//---about vectors (or points)//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
V ans[10];//declared for other functions
int tot;
inline void input(V a){scanf(%lf%lf,a.x,a.y);}
void print(const V a,int op1){printf(%.10lf %.10lf,a.x,a.y);putchar(op?10:32);}
//op:endl or space//basic operation
inline V operator (const V a,const V b){return (V){a.xb.x,a.yb.y};}
inline V operator - (const V a,const V b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double x,const V a){return (V){a.x*x,a.y*x};}
inline V operator * (const V a,const double x){return (V){a.x*x,a.y*x};}
inline V operator / (const V a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator (const V a,const V b){return abs(a.x-b.x)epsabs(a.y-b.y)eps;}
inline bool operator ! (const V a,const V b){return !(ab);}
inline double operator * (const V a,const V b){return a.x*b.xa.y*b.y;}
inline double operator ^ (const V a,const V b){return a.x*b.y-a.y*b.x;}
inline double len(const V a){return sqrt(a.x*a.xa.y*a.y);}
inline V mid(const V a,const V b){return (V){(a.xb.x)/2,(a.yb.y)/2};}
inline V chui(const V a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V a){return a/len(a);}
inline double tri_S(const V a,const V b,const V c){return abs((b-a)^(c-a))/2;}//always be non-negative
inline bool operator (const V a,const V b){return a.xb.x-eps||(abs(a.x-b.x)epsa.yb.y-eps);
}
inline double ang(const V a,const V b){return acos((a*b)/len(a)/len(b));}
inline V rotate(const V o,double t){//COUNTER_CLOCKWISEdouble ssin(t),ccos(t);return (V){o.x*c-o.y*s,o.x*so.y*c};
}const int N1e5100;
const int M505;
int n,m;
V p[N],zhan[N];
bool cmp(V a,V b){double d(a-p[1])^(b-p[1]);if(abs(d)eps) return d0;else return len(a-p[1])len(b-p[1]);
}
void graham(V *p,int n){int top0;sort(p1,p1n);sort(p2,p1n,cmp);top0;for(int i1;in;i){while((top1((zhan[top]-zhan[top-1])^(p[i]-zhan[top]))0)) --top;zhan[top]p[i];}memcpy(p,zhan,sizeof(zhan));ntop;return;
}
inline ll calc(const V a,const V b){return (a.x-b.x)*(a.x-b.x)(a.y-b.y)*(a.y-b.y)eps;
}
ll rotate_calipers(V *p,int n){ll res(0);int pl1;for(int i2;in;i){if(((p[2]-p[1])^(p[pl]-p[2]))((p[2]-p[1])^(p[i]-p[2]))) pli;}for(int i1;in;i){while(((p[i1]-p[i])^(p[pl]-p[i1]))((p[i1]-p[i])^(p[pl1]-p[i1]))){pl(pl1)%n;resmax(res,max(calc(p[i],p[pl]),calc(p[i1],p[pl])));}resmax(res,max(calc(p[i],p[pl]),calc(p[i1],p[pl])));}return res;
}
signed main(){
#ifndef ONLINE_JUDGEfreopen(a.in,r,stdin);freopen(a.out,w,stdout);
#endifnread();for(int i1;in;i) input(p[i]);graham(p,n);p[0]p[n];p[n1]p[1];//for(int i1;in;i) print(p[i]);//putchar(\n);printf(%lld\n,rotate_calipers(p,n));return 0;
}
/*
3 5
0 -2
-5 3
0 -7
*/
半平面交
感觉应用很灵活的一个算法一切有两个变量的线性规划问题都可以转化为半平面交。 有时可能要注意取等问题指射箭但多数时候似乎并不用。
半平面的表示
首先我们需要一个简洁的方式表达半平面。 不禁想到我们可以利用叉积的符号来判断点是在一条线的顺时针还是逆时针方向。那么类似的我们可以用一个点x和一个方向向量d表示一个半平面一个点p属于这个半平面当且仅当 d→×XP→0\overrightarrow{d}\times \overrightarrow{XP}0d×XP0也就是点在方向向量的左侧是否取等取决于半平面本身的开闭。
inline bool in(const line l,const V o){return (l.d^(o-l.a))eps;}
inline bool out(const line l,const V o){return (l.d^(o-l.a))-eps;}流程
有了一个合适的表示方法后我们就可以开始尝试解决半平面交问题了顾名思义这个问题就是 给出若干个半平面求它们的交。 显然得到的一定是一个凸包。 怎么做呢 首先我们把所有点半平面按照极角排序极角可以利用 atan2(y,x) 函数很方便的求出。 当两个向量共线时显然只有靠左的那个是有用的我们让有用的放在前面。这样在后面的算法中遇到没用的时候可以直接跳过避免求平行线的交点而出现错误
bool cmp(const line l1,const line l2){double t1atan2(l1.d.y,l1.d.x),t2atan2(l2.d.y,l2.d.x);if(abs(t1-t2)eps) return t1t2;return in(l2,l1.a);
}考虑维护一个双端队列来维护当前对答案有用的半平面并记录每个半平面和前一个半平面的交点。 维护分为如下步骤
每加入一个半平面时 若队尾与队尾前一个的交点在当前半平面外则不断弹出队尾。类似的若队首和后一个的交点在当前半平面外也不断弹出队首。加入当前半平面并计算其与上一个半平面的交点。插入所有半平面后若队尾与上一个的交点在队首半平面外那么队尾也是无用的弹出队尾所有没有用的半平面。最后插入队尾与队首的交点得到的封闭的凸包就是答案。
void half_plane(){sort(l1,l1tot,cmp);st1,ed0;for(int i1;itot;i){if(stedabs(q[ed].d^l[i].d)eps) continue;//和上一个共线必然是无用向量 while(ed-st12out(l[i],pt[ed])) --ed;//先弹队尾 while(ed-st12out(l[i],pt[st1])) st;q[ed]l[i];if(ed-st12) pt[ed]jiaodian(q[ed-1],q[ed]);}while(ed-st12out(q[st],pt[ed])) --ed;//弹出队尾无用向量 if(ed-st12) pt[ed1]jiaodian(q[st],q[ed]),ed;//加入队尾与队首的交点 med-st;for(int i1;im;i) a[i]pt[sti];a[m1]a[1];return;
}细节 可不可以先弹队首再弹队尾 不可以 考虑这两种写法不一样的情况就是当我判断在半平面外的那个交点恰好就是队首和队尾的交点时。 也就是 图片来自OI-wiki 半平面交
显然这个时候弹队首会出错。 如何判断交集为空或无边界 考虑在正无穷远处加入四条与坐标轴平行的边把整个平面框成一个大矩形在进行半平面交算法。 这样就不会出现无边界情况要判断也可以通过凸包内是否出现了那四条边来进行。 而且这样之后凸包元素不足三个就变成了答案无解的充要条件。 非常优美简洁。
代码
P4196 [CQOI2006]凸多边形 /【模板】半平面交
#includebits/stdc.h
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
inline ll read(){ll x(0),f(1);char cgetchar();while(!isdigit(c)){if(c-)f-1;cgetchar();}while(isdigit(c)){x(x1)(x3)c-0;cgetchar();}return x*f;
}//basic declare
//#define double long double
const double eps1e-10;
const double piacos(-1.0);//---about vectors (or points)//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
inline void input(V a){scanf(%lf%lf,a.x,a.y);}
void print(const V a,int op1){printf(%g %g,a.x,a.y);putchar(op?10:32);}
//op:endl or space//basic operation
inline V operator (const V a,const V b){return (V){a.xb.x,a.yb.y};}
inline V operator - (const V a,const V b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double x,const V a){return (V){a.x*x,a.y*x};}
inline V operator * (const V a,const double x){return (V){a.x*x,a.y*x};}
inline V operator / (const V a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator (const V a,const V b){return abs(a.x-b.x)epsabs(a.y-b.y)eps;}
inline bool operator ! (const V a,const V b){return !(ab);}
inline double operator * (const V a,const V b){return a.x*b.xa.y*b.y;}
inline double operator ^ (const V a,const V b){return a.x*b.y-a.y*b.x;}
inline double len(const V a){return sqrt(a.x*a.xa.y*a.y);}
inline V mid(const V a,const V b){return (V){(a.xb.x)/2,(a.yb.y)/2};}
inline V chui(const V a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V a){return a/len(a);}
inline double ang(const V a,const V b){return acos((a*b)/len(a)/len(b));}
inline V rotate(const V o,double t){//COUNTER_CLOCKWISEdouble ssin(t),ccos(t);return (V){o.x*c-o.y*s,o.x*so.y*c};
}
struct line{V d,a;
};
void print(const line l,int op1){printf(point: );print(l.a,0);printf(dir: );print(l.d,op);
}
inline line trans(const V a,const V b){return (line){b-a,a};
}
inline line trans(const double a,const double b,const double c,const double d){
//given one point and directionreturn (line){(V){c,d},(V){a,b}};
}
inline V jiaodian(const line u,const line v){if(u.dv.d){print(u);print(v);exit(0);}double k((v.a-u.a)^v.d)/(u.d^v.d);return u.ak*u.d;
}inline bool in(const line l,const V o){return (l.d^(o-l.a))0;}
inline bool out(const line l,const V o){return (l.d^(o-l.a))-eps;}const int N1005;
int n,m;int tot,st,ed;
V a[N],pt[N];
line l[N],q[N];bool cmp(const line l1,const line l2){double t1atan2(l1.d.y,l1.d.x),t2atan2(l2.d.y,l2.d.x);if(abs(t1-t2)eps) return t1t2;return in(l2,l1.a);
}const double inf1e9;
void half_plane(){l[tot]trans(0,0,1,0);//插入边界l[tot]trans(X,0,0,1);l[tot]trans(0,Y,-1,0);l[tot]trans(0,0,0,-1);sort(l1,l1tot,cmp);st1,ed0;for(int i1;itot;i){if(stedabs(q[ed].d^l[i].d)eps) continue;//和上一个共线必然是无用向量 while(ed-st12out(l[i],pt[ed])) --ed;//先弹队尾 while(ed-st12out(l[i],pt[st1])) st;q[ed]l[i];if(ed-st12) pt[ed]jiaodian(q[ed-1],q[ed]);}while(ed-st12out(q[st],pt[ed])) --ed;//弹出队尾无用向量 if(ed-st12) pt[ed1]jiaodian(q[st],q[ed]),ed;//加入队尾与队首的交点 med-st;for(int i1;im;i) a[i]pt[sti];a[m1]a[1];return;
}
signed main(){
#ifndef ONLINE_JUDGEfreopen(a.in,r,stdin);freopen(a.out,w,stdout);
#endifnread();for(int i1;in;i){mread();for(int j1;jm;j) input(a[j]);a[m1]a[1];for(int j1;jm;j) l[tot]trans(a[j],a[j1]);}half_plane();if(m3){printf(0.000);return 0;}double res(0);for(int i1;im;i) res(a[i]^a[i1])/2;if(res0) res-res;printf(%.3lf\n,res);return 0;
}
/*
4
2 8
4 2
6 5
-8 2
*/
其他
皮克定理
条件整点任意多边形。 内容面积整点数1/2*边点数-1。 为什么减一可以考虑极限情况多边形退化成单格点时面积是0 拓展对于平行四边形的格点依然成立对于三角形的格点面积要乘二即2面积整点数1/2边点数-1可以理解成变成三角形后面积减少了一半所以要乘二才能相等。 距离 规定以下的 ddd 表示维度xi,dx_{i,d}xi,d 表示点 iii 在第 ddd 维的坐标。 欧几里得距离 ∑d(xi,d−xj,d)2\sqrt{\sum_d (x_{i,d}-x_{j,d})^2 }∑d(xi,d−xj,d)2 曼哈顿距离∑d∣xi,d−xj,d∣\sum_d|x_{i,d-}x_{j,d}|∑d∣xi,d−xj,d∣ 切比雪夫距离maxd∣xi,d−xj,d∣\max_d|x_{i,d-}x_{j,d}|maxd∣xi,d−xj,d∣
曼哈顿距离和切比雪夫距离的转化
曼哈顿距离下的点 (x,y)(x,y)(x,y) 等价于切比雪夫距离下的点 (xy,x−y)(xy,x-y)(xy,x−y)。 反过来切比雪夫下的点 (x,y)(x,y)(x,y) 等价于切比雪夫距离下的点 (xy2,x−y2)(\dfrac{xy}{2},\dfrac{x-y}{2})(2xy,2x−y)。 证明暴力拆分曼哈顿的定义或者结合单位正方形。
source
//basic declare
//#define double long double
const double eps1e-10;
const double piacos(-1.0);//---about vectors (or points)//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
V ans[10];//declared for other functions
int tot;
inline void input(V a){a.xread();a.yread();}
void print(const V a,int op1){printf(%.10lf %.10lf,a.x,a.y);putchar(op?10:32);}
//op:endl or space//basic operation
inline V operator (const V a,const V b){return (V){a.xb.x,a.yb.y};}
inline V operator - (const V a,const V b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double x,const V a){return (V){a.x*x,a.y*x};}
inline V operator * (const V a,const double x){return (V){a.x*x,a.y*x};}
inline V operator / (const V a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator (const V a,const V b){return abs(a.x-b.x)epsabs(a.y-b.y)eps;}
inline bool operator ! (const V a,const V b){return !(ab);}
inline double operator * (const V a,const V b){return a.x*b.xa.y*b.y;}
inline double operator ^ (const V a,const V b){return a.x*b.y-a.y*b.x;}
inline double len(const V a){return sqrt(a.x*a.xa.y*a.y);}
inline V mid(const V a,const V b){return (V){(a.xb.x)/2,(a.yb.y)/2};}
inline V chui(const V a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V a){return a/len(a);}
inline double tri_S(const V a,const V b,const V c){return abs((b-a)^(c-a))/2;}//always be non-negative
inline bool operator (const V a,const V b){return a.xb.x-eps||(abs(a.x-b.x)epsa.yb.y-eps);
}//operations about angle and vectors
inline double ang(const V a,const V b){return acos((a*b)/len(a)/len(b));}
inline bool dun(const V a,const V b,const V c){return ((b-a)*(c-a))-eps;}//angle:BAC
inline bool rui(const V a,const V b,const V c){return ((b-a)*(c-a))eps;}
inline bool zhi(const V a,const V b,const V c){return abs(((b-a)*(c-a)))eps;}
inline V rotate(const V o,double t){//COUNTER_CLOCKWISEdouble ssin(t),ccos(t);return (V){o.x*c-o.y*s,o.x*so.y*c};
}//---about lines//definition
struct line{V d,a,b;
};
inline line trans(double a,double b,double c,double d){//given coordinatesV dd(c-a,d-b),x(a,b),y(c,d);dddd/len(dd);return (line){dd,x,y};
}
inline line trans(const V a,const V b){//given pointsV dd(b.x-a.x,b.y-a.y);dddd/len(dd);return (line){dd,a,b};
}
inline void input(line l){double aread(),bread(),cread(),dread();ltrans(a,b,c,d);return;
}//functions about points and lines
inline V chui(const line l,const V o){//directedreturn l.a((o-l.a)*l.d)*l.d;
}
inline V duichen(const line l,const V o){V uchui(l,o);return (V){2*u.x-o.x,2*u.y-o.y};
}
inline bool on_line(const line l,const V o){return abs((l.d^(o-l.a)))eps;}
inline bool on_seg(const line l,const V o){return abs(len(o-l.a)len(o-l.b)-len(l.a-l.b))eps;
}
inline int pos(const line l,const V o){if(!on_line(l,o)){if((l.d^(o-l.a))0) return 1;//counter clockwiseelse return 2;//clockwise}else{if(((o-l.a)*(o-l.b))0) return 5;//on segmentelse if(((o-l.a)*l.d)0) return 3;//online backelse return 4;//online front}
}
inline double dis(const line l,const V o,int op0){//op0:dis on line,op1:dis on segmentif(op(dun(l.a,o,l.b)||dun(l.b,o,l.a))) return min(len(o-l.a),len(o-l.b));else return abs((o-l.a)^(o-l.b))/len(l.a-l.b);
}//functions about lines and lines
inline bool gongxian(const line a,const line b){return abs(a.d^b.d)eps;}
inline bool chuizhi(const line a,const line b){return abs(a.d*b.d)eps;}
inline bool jiao(const line u,const line v){if(min(u.a.x,u.b.x)max(v.a.x,v.b.x)eps||max(u.a.x,u.b.x)min(v.a.x,v.b.x)-eps||min(u.a.y,u.b.y)max(v.a.y,v.b.y)eps||max(u.a.y,u.b.y)min(v.a.y,v.b.y)-eps) return false;//rapid rejection testreturn ((u.a-v.a)^v.d)*((u.b-v.a)^v.d)eps((v.a-u.a)^u.d)*((v.b-u.a)^u.d)eps;//straddle test
}
inline double Sin(const V a,const V b){return abs(a^b)/len(a)/len(b);}
inline V jiaodian(const line u,const line v){if(on_line(v,u.a)) return u.a;if(on_line(v,u.b)) return u.b;double s1tri_S(u.a,v.a,v.b),s2tri_S(u.b,v.a,v.b);return u.a((s1/(s1s2))*(u.b-u.a));
}
inline line pingfen(const V a,const V b,const V c){//angle:BACV d1danwei(b-a),d2danwei(c-a),ddanwei(d1d2);return (line){d,a,ad};
}
inline line zhongchui(const V a,const V b){V xmid(a,b),ddanwei(chui(b-a));return (line){d,x,xd};
}//---about polygons//definition
//V a[N];
//int n;//funtions
double S(const V *a,const int n){double res(0);for(int i1;in;i) res(a[i]^a[i%n1]);return res/2;
}
bool is_convex(const V *a,const int n){int op0;for(int i1;in;i){double o(a[i1]-a[i])^(a[i]-a[i-1]);if(abs(o)eps) continue;//three neighbouring points on the same lineint nopo0?1:-1;if(!op) opnop;else if(op!nop) return false;}return true;
}
int in_poly(const V *a,const int n,const V o){int res(0);for(int i1;in;i){V ua[i],va[i1];if(on_seg(trans(u,v),o)) return 1;//on the edgeif(abs(u.y-v.y)eps) continue;//ignore horizontalif(max(u.y,v.y)o.y-eps) continue;//equal will not continueif(min(u.y,v.y)o.y-eps) continue;//equal will continuedouble xu.x(o.y-u.y)/(v.y-u.y)*(v.x-u.x);if(xo.x) res^1;}return res?2:0;//odd:in even:out
}//---about circles//definition
struct cir{V o;double r;cir(double a0,double b0,double c0):o(a,b),r(c){}cir(V O,double c0):o(O),r(c){}
};
inline void input(cir cc){input(cc.o);cc.rread();}//functions about position relation
inline bool incir(const cir c,const V p){return len(p-c.o)c.r-eps;}
inline bool oncir(const cir c,const V p){return (len(p-c.o)-c.r)eps;}
inline bool outcir(const cir c,const V p){return len(p-c.o)c.reps;}
inline double dis(const cir c,const line l){return dis(l,c.o);}
inline int pos(const cir c,const line l){double ddis(c,l);if(dc.r-eps) return 1;//intersectelse if(abs(d-c.r)eps) return 2;//tangentelse if(dc.reps) return 3;//disjoint
}
inline double dis(const cir c1,const cir c2){return len(c1.o-c2.o);}
inline int pos(const cir c1,const cir c2){double ddis(c1,c2);if(dc1.rc2.reps) return 4;//do not crosselse if(abs(d-c1.r-c2.r)eps) return 3;//circumscribedelse if(max(c1.r,c2.r)min(c1.r,c2.r)d -eps) return 2;//intersectelse if(abs(max(c1.r,c2.r)-min(c1.r,c2.r)-d)eps) return 1;//inscribedelse return 0;//include
}//functions about circles and triangles
inline cir incir(const V a,const V b,const V c){V x(jiaodian(pingfen(a,b,c),pingfen(b,a,c)));return cir(x,dis(trans(a,b),x));
}
inline cir outcir(const V a,const V b,const V c){V x(jiaodian(zhongchui(a,b),zhongchui(a,c)));return cir(x,len(x-a));
}//functions about cross points
inline double calc(double xie,double zhi){return sqrt(xie*xie-zhi*zhi);}//Pythagorean Theorem
inline void cross_line_cir(const cir c,const line l){tot0;V xchui(l,c.o);double ddlen(x-c.o);if(c.rdd-eps) return;//none cross pointif(abs(c.r-dd)eps){//one cross pointans[tot]x;return;}double discalc(c.r,dd);V axdis*l.d,bx-dis*l.d;//two cross pointsans[tot]a;ans[tot]b;return;
}
inline void cross_cir_cir(const cir c1,const cir c2){int oppos(c1,c2);if(op4||op0) return;//none cross pointV Lc2.o-c1.o;double ac2.r,bc1.r,clen(L);double tacos((b*bc*c-a*a)/(2*b*c));//find the angleV xc1.oc1.r*rotate((L)/len(L),t),yduichen(trans(c1.o,c2.o),x);if(xy) print(x,0),print(y,1);else print(y,0),print(x,1);
}//functions about tangent points and lines
inline void qiedian(const cir c,const V p){tot0;line Ltrans(c.o,p);double tacos(c.r/len(c.o-p));V ac.o(c.r*rotate(L.d,t)),bduichen(L,a);ans[tot]a;if(a!b) ans[tot]b;return;
}
inline void common_qie(const cir c1,const cir c2){tot0;int oppos(c1,c2);if(op){//outside tangent linesdouble dc1.r-c2.r,tacos(d/len(c1.o-c2.o));V uc1.o(c1.r*danwei(rotate(c2.o-c1.o,t))),vc1.o(c1.r*danwei(rotate(c2.o-c1.o,-t)));ans[tot]u;if(u!v) ans[tot]v; }if(op2){//inside tangent linesdouble dlen(c2.o-c1.o)/(c1.rc2.r)*c1.r,tacos(c1.r/d);V uc1.o(c1.r*danwei(rotate(c2.o-c1.o,t))),vc1.o(c1.r*danwei(rotate(c2.o-c1.o,-t)));ans[tot]u;if(u!v) ans[tot]v;}
}//functions about area
inline double cir_S(const cir c){return pi*c.r*c.r;}
inline double shan(const cir c,const V a,const V b){//S of sectordouble tang(a-c.o,b-c.o);return t/2*c.r*c.r;
}
inline double gong(const cir c,const V a,const V b){//S of bowreturn shan(c,a,b)-tri_S(c.o,a,b);
}
inline double O_tri(const cir c,V a,V b){
//directed S of Intersection between circle O and triangle OABif(on_line(trans(a,b),c.o)) return 0.0;int f(((a-c.o)^(b-c.o))0)?1:-1;//directionline ltrans(a,b);int f1incir(c,a),f2incir(c,b);if(f1f2) return f*tri_S(a,b,c.o);//both incircle:the S of triangleelse if(!f1!f2){//both outcircle:a sector(possibly minus a bow)V u(c.o(c.r*danwei(a-c.o))),v(c.o(c.r*danwei(b-c.o)));double Sshan(c,u,v); if(dis(l,c.o,1)c.r-eps){cross_line_cir(c,l);assert(tot2);S-gong(c,ans[1],ans[2]);}return f*S;}else{//one in and one out:a traingle and a sectorif(!f1) swap(a,b);double saSin(b-a,c.o-a),susa/c.r*len(c.o-a),Aasin(sa),Uasin(su);if((b-a)*(c.o-a)-eps) Api-A;double tpi-A-U,dissin(t)/sa*c.r; V ua(dis*danwei(b-a)),vc.r*danwei(b-c.o);return f*(tri_S(c.o,a,u)shan(c,u,v));}
}
double common_cir_poly(const V *a,const cir c,int n){double res(0);for(int i1;in;i) resO_tri(c,a[i],a[i1]);return res;
}
inline double common_cir_cir(const cir c1,const cir c2){int oppos(c1,c2);if(op3) return 0;//common area0else if(op1) return min(cir_S(c1),cir_S(c2));//completly includeelse{//partly includedouble Llen(c1.o-c2.o);double t12*acos((L*Lc1.r*c1.r-c2.r*c2.r)/(2*L*c1.r));double t22*acos((L*Lc2.r*c2.r-c1.r*c1.r)/(2*L*c2.r));return c1.r*c1.r*t1/2-c1.r*c1.r*sin(t1)/2c2.r*c2.r*t2/2-c2.r*c2.r*sin(t2)/2;}
}
