计算几何笔记
课堂笔记
线段与直线
- 求点到直线的垂足
$\begin{array}{l}
1|\overrightarrow{A C}|=\overrightarrow{A B} \cdot \hat{d} \\
2 C=A+\overrightarrow{A C}
\end{array}$
- 求点关于直线的对称点
- 求出垂足之后用中点公式
线段与直线
- 点与直线的位置关系
1 叉积 $\overrightarrow{A B} \times \hat{d}$ : 为零, 点在直线上; 不为零, 点在直线外。
2 点在直线外时: 若 $\overrightarrow{A B} \times \hat{d}<0$ , 则点在直线的左侧(逆时 针); 反之在右侧 (顺时针)。
3 点在直线上时: 若点积 $\overrightarrow{A B} \cdot \hat{d}<0$ , 则点在标志点 $A$ 的后方; 反之在前方。
用此可以判断某一个点是否在某条线段上。
线段与直线
- 线段与线段的位置关系
1 快速排斥试验: 分别找到包含两线段的最小矩形, 判断两矩 形是否有公共部分。
2 跨立试验: 判断是否每一线段的两端点位于另一线段的两侧。
3 同时通过上述两个试验, 则两线段相交; 反之不相交。
- 正确性:跨立试验:没有共线情况下是充要的(规范相交); 共线。
- 角平分线
- 一个顶点 A , 两个单位方向向量 $\hat{d}_{1}, \hat{d}_{2} $, 则其角平分线为 $\left(A, \frac{\hat{d}_{1}+\hat{d}_{2}}{\left(\hat{d}_{1}+\hat{d}_{2} \mid\right.}\right) $。
多边形
- 判断点在多边形内
- 夹角和判别法
- 面积和判别法
- 引射线法
- 夹角和判别法
1 计算夹角和: 多边形上任意相邻点与该点构成的有向角度之 和。设 P 为要判定的点:
$S_{\text {angle }}=\sum_{i=1}^{n} \operatorname{rrccos}\left(\frac{\overrightarrow{P A}_{i} \cdot \overrightarrow{P A}_{i \text { mod } n+1}}{\left|\overrightarrow{P A}_{i}\right|\left|\overrightarrow{P A}_{i \bmod n+1}\right|}\right)$
2 $S_{\text {angle }}$ 为 $0$ 或 $2 \pi $, 当 $S_{\text {angle }}=0$ 时, $P$ 在多边形外部, 反之在多边形内部。
3 $P$ 点在多边形边界上可能需要特判。
多边形
- 面积和判别法
1 计算面积和: 多边形上任意相邻两点与该点构成的无向面积 之和。设 $P$ 为要判定的点:
$S_{\text {area }}=\frac{1}{2} \sum_{i=1}^{n}\left|\overrightarrow{P A}_{i} \times \overrightarrow{P A}_{i \bmod n+1}\right|$
2 计算多边形面积 $S$ 。
3 若 $S_{\text {area }}=S$ , 则 $P$ 在多边形内部; 反之在外部, 此时 $S_{\text {area }}>S$ 。
4 $P$ 在多边形边界上可能需要特判。
多边形
- 引射线法
-从 $P$ 向任意方向引出一条射线。当这条射线从无穷远处驶向 $ P$ 时, 每穿过一次多边形的边界, 其是否在多边形内的状态 将发生一次转变。故 $P$ 点是否在多边形内可以由这条射线穿 过边界次数的奇偶性判断。
1 随机选取一个射线方向(目的是避免穿过多边形顶点)。
2 枚举多边形的每条线段, 判断射线是否与线段相交。
3 若相交次数为奇数, 则 $P$ 在多边形内; 反之在多边形外。
4 $P$ 在多边形边界上可能需要特判。
圆
- 圆的表示方法: 圆心 + 半径
- 点与圆的位置关系
- 根据点到圆心的距离直接判断(圆内、圆外、圆上)
- 线与圆的位置关系
- 根据圆心到直线的距离直接判断(相截、相切、相离)
- 圆与圆的位置关系
- 根据两圆心之间的距离直接判断(包含、截交、内切、外切、 相离)
圆
- 三角形的内接圆
- 圆心: 两条角平分线的交点
- 半径: $r=\frac{2 S \Delta}{a+b+c}$
或者调用圆心到某条边的距离
- 三角形的外接圆
- 圆心: 两条中垂线的交点
追求更高的精度直接解方程可能更好
- 半径: 调用圆心到顶点的距离
圆
- 点到圆的切点
1 点、圆心、切点构成一个直角三角形
2 求出夹角后, 沿两个方向旋转得到两个切点
- 圆与圆的公切线
- 内公切线与外公切线类似, 只有符号的区别
1 两个圆心、两个切点构成梯形 (不一定简单)
2 求出圆心到切点、另一个圆心的张角, 沿两个方向旋转得到两个切点。
凸包算法
- 凸包的定义: 能够包围平面上 $n$ 个点的周长最小的多边形。 可以证明它一定是凸的(将一个凹多边形修改成凸的一定不会使周长变长)。
- 凸包的常用求法
- 暴力求解法 $O\left(N^{\beta}\right)$
- Jarvis 步进法 $O(N M)$, $M$ 为凸包上点的个数
- Graham 扫描法 $O(N \log N)$
- Andrew 算法 $O(N \log N)$
- 分治法 $O(N \log N)-O\left(N^{2}\right)$
- Andrew 算法
- 核心: 维护上凸壳与下凸壳
1 将所有点以横坐标为第一关键字, 纵坐标为第二关键字排序。 最小和最大的点一定在凸包上。
2 沿逆时针方向看, 下凸壳从最小点出发一路左拐走到最大点, 上凸壳从最大点出发一路左拐走到最小的点。
3 先求下凸壳, 同样用栈来维护。新加入点的时候判断是否为 左拐, 是则加入; 否则就弹出栈顶的点直至构成左拐。
4 上凸壳同理, 只是加点的顺序. 改为倒序。
5 将上凸壳与下凸壳合并得到整个凸包。
6 总复杂度 $O(N \log N)$ , 瓶颈为排序。
旋转卡壳
- 在求出凸包之后, 沿着凸包用两条平行线, 在旋转的过程中 不断地将凸包 “卡住”。以此在线性的时间内求解一些和凸包性质相关的问题。
- 经典问题
- 最远点对(凸包直径)
- 最小矩形覆盖
Acwing整理知识点(yxcyyds): 1. 前置知识点 (1) pi = acos(-1); (2) 余弦定理 c^2 = a^2 + b^2 - 2abcos(t) 2. 浮点数的比较 const double eps = 1e-8; int sign(double x) // 符号函数 { if (fabs(x) < eps) return 0; if (x < 0) return -1; return 1; } int cmp(double x, double y) // 比较函数 { if (fabs(x - y) < eps) return 0; if (x < y) return -1; return 1; } 3. 向量 3.1 向量的加减法和数乘运算 3.2 内积(点积) A·B = |A||B|cos(C) (1) 几何意义:向量A在向量B上的投影与B的长度的乘积。 (2) 代码实现 double dot(Point a, Point b) { return a.x * b.x + a.y * b.y; } 3.3 外积(叉积) AxB = |A||B|sin(C) (1) 几何意义:向量A与B张成的平行四边形的有向面积。B在A的逆时针方向为正。 (2) 代码实现 double cross(Point a, Point b) { return a.x * b.y - b.x * a.y; } 3.4 常用函数 3.4.1 取模 double get_length(Point a) { return sqrt(dot(a, a)); } 3.4.2 计算向量夹角 double get_angle(Point a, Point b) { return acos(dot(a, b) / get_length(a) / get_length(b)); } 3.4.3 计算两个向量构成的平行四边形有向面积 double area(Point a, Point b, Point c) { return cross(b - a, c - a); } 3.4.5 向量A顺时针旋转C的角度: Point rotate(Point a, double angle) { return Point(a.x * cos(angle) + a.y * sin(angle), -a.x * sin(angle) + a.y * cos(angle)); } 4. 点与线 4.1 直线定理 (1) 一般式 ax + by + c = 0 (2) 点向式 p0 + vt (3) 斜截式 y = kx + b 4.2 常用操作 (1) 判断点在直线上 A x B = 0 (2) 两直线相交 // cross(v, w) == 0则两直线平行或者重合 Point get_line_intersection(Point p, Vector v, Point q, vector w) { vector u = p - q; double t = cross(w, u) / cross(v, w); return p + v * t; } (3) 点到直线的距离 double distance_to_line(Point p, Point a, Point b) { vector v1 = b - a, v2 = p - a; return fabs(cross(v1, v2) / get_length(v1)); } (4) 点到线段的距离 double distance_to_segment(Point p, Point a, Point b) { if (a == b) return get_length(p - a); Vector v1 = b - a, v2 = p - a, v3 = p - b; if (sign(dot(v1, v2)) < 0) return get_length(v2); if (sign(dot(v1, v3)) > 0) return get_length(v3); return distance_to_line(p, a, b); } (5) 点在直线上的投影 double get_line_projection(Point p, Point a, Point b) { Vector v = b - a; return a + v * (dot(v, p - a) / dot(v, v)); } (6) 点是否在线段上 bool on_segment(Point p, Point a, Point b) { return sign(cross(p - a, p - b)) == 0 && sign(dot(p - a, p - b)) <= 0; } (7) 判断两线段是否相交 bool segment_intersection(Point a1, Point a2, Point b1, Point b2) { double c1 = cross(a2 - a1, b1 - a1), c2 = cross(a2 - a1, b2 - a1); double c3 = cross(b2 - b1, a2 - b1), c4 = cross(b2 - b1, a1 - b1); return sign(c1) * sign(c2) <= 0 && sign(c3) * sign(c4) <= 0; } 5. 多边形 5.1 三角形 5.1.1 面积 (1) 叉积 (2) 海伦公式 p = (a + b + c) / 2; S = sqrt(p(p - a) * (p - b) * (p - c)); 5.1.2 三角形四心 (1) 外心,外接圆圆心 三边中垂线交点。到三角形三个顶点的距离相等 (2) 内心,内切圆圆心 角平分线交点,到三边距离相等 (3) 垂心 三条垂线交点 (4) 重心 三条中线交点(到三角形三顶点距离的平方和最小的点,三角形内到三边距离之积最大的点) 5.2 普通多边形 通常按逆时针存储所有点 5.2.1 定义 (1) 多边形 由在同一平面且不再同一直线上的多条线段首尾顺次连接且不相交所组成的图形叫多边形 (2) 简单多边形 简单多边形是除相邻边外其它边不相交的多边形 (3) 凸多边形 过多边形的任意一边做一条直线,如果其他各个顶点都在这条直线的同侧,则把这个多边形叫做凸多边形 任意凸多边形外角和均为360° 任意凸多边形内角和为(n−2)180° 5.2.2 常用函数 (1) 求多边形面积(不一定是凸多边形) 我们可以从第一个顶点除法把凸多边形分成n − 2个三角形,然后把面积加起来。 double polygon_area(Point p[], int n) { double s = 0; for (int i = 1; i + 1 < n; i ++ ) s += cross(p[i] - p[0], p[i + 1] - p[i]); return s / 2; } (2) 判断点是否在多边形内(不一定是凸多边形) a. 射线法,从该点任意做一条和所有边都不平行的射线。交点个数为偶数,则在多边形外,为奇数,则在多边形内。 b. 转角法 (3) 判断点是否在凸多边形内 只需判断点是否在所有边的左边(逆时针存储多边形)。 5.3 皮克定理 皮克定理是指一个计算点阵中顶点在格点上的多边形面积公式该公式可以表示为: S = a + b/2 - 1 其中a表示多边形内部的点数,b表示多边形边界上的点数,S表示多边形的面积。 6. 圆 (1) 圆与直线交点 (2) 两圆交点 (3) 点到圆的切线 (4) 两圆公切线 (5) 两圆相交面积 #include <iostream> #include <cstring> #include <algorithm> #define x first #define y second using namespace std; typedef long long LL; typedef pair<LL, LL> PLL; const int N = 5010; int n, m; PLL a[N], b[N]; int ans[N]; LL cross(LL x1, LL y1, LL x2, LL y2) { return x1 * y2 - x2 * y1; } LL area(PLL a, PLL b, PLL c) { return cross(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y); } int find(LL x, LL y) { int l = 0, r = n; while (l < r) { int mid = l + r >> 1; if (area(b[mid], a[mid], {x, y}) > 0) r = mid; else l = mid + 1; } return r; } int main() { bool is_first = true; while (scanf("%d", &n), n) { LL x1, y1, x2, y2; scanf("%d%lld%lld%lld%lld", &m, &x1, &y1, &x2, &y2); for (int i = 0; i < n; i ++ ) { LL u, l; scanf("%lld%lld", &u, &l); a[i] = {u, y1}, b[i] = {l, y2}; } a[n] = {x2, y1}, b[n] = {x2, y2}; if (is_first) is_first = false; else puts(""); memset(ans, 0, sizeof ans); while (m -- ) { LL x, y; scanf("%lld%lld", &x, &y); ans[find(x, y)] ++ ; } for (int i = 0; i <= n; i ++ ) printf("%d: %d\n", i, ans[i]); } return 0; }
XingLing的整体模板
#include<algorithm>
#include<cstdio>
#include<cmath>
/*一:【准备工作】*/
#define LD double
#define LL long long
#define Re register int
#define Vector Point
using namespace std;
const int N=262144+3;
const LD eps=1e-8,Pi=acos(-1.0);
inline int dcmp(LD a){return a<-eps?-1:(a>eps?1:0);}//处理精度
inline LD Abs(LD a){return a*dcmp(a);}//取绝对值
struct Point{
LD x,y;Point(LD X=0,LD Y=0){x=X,y=Y;}
inline void in(){scanf("%lf%lf",&x,&y);}
inline void out(){printf("%.2lf %.2lf\n",x,y);}
};
/*二:【向量】*/
inline LD Dot(Vector a,Vector b){return a.x*b.x+a.y*b.y;}//【点积】
inline LD Cro(Vector a,Vector b){return a.x*b.y-a.y*b.x;}//【叉积】
inline LD Len(Vector a){return sqrt(Dot(a,a));}//【模长】
inline LD Angle(Vector a,Vector b){return acos(Dot(a,b)/Len(a)/Len(b));}//【两向量夹角】
inline Vector Normal(Vector a){return Vector(-a.y,a.x);}//【法向量】
inline Vector operator+(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
inline Vector operator-(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
inline Vector operator*(Vector a,LD b){return Vector(a.x*b,a.y*b);}
inline bool operator==(Point a,Point b){return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);}//两点坐标重合则相等
/*三:【点、向量的位置变换】*/
/*1.【点、向量的旋转】*/
inline Point turn_P(Point a,LD theta){//【点A\向量A顺时针旋转theta(弧度)】
LD x=a.x*cos(theta)+a.y*sin(theta);
LD y=-a.x*sin(theta)+a.y*cos(theta);
return Point(x,y);
}
inline Point turn_PP(Point a,Point b,LD theta){//【将点A绕点B顺时针旋转theta(弧度)】
LD x=(a.x-b.x)*cos(theta)+(a.y-b.y)*sin(theta)+b.x;
LD y=-(a.x-b.x)*sin(theta)+(a.y-b.y)*cos(theta)+b.y;
return Point(x,y);
}
/*四:【图形与图形之间的关系】*/
/*1.【点与线段】*/
inline int pan_PL(Point p,Point a,Point b){//【判断点P是否在线段AB上】
return !dcmp(Cro(p-a,b-a))&&dcmp(Dot(p-a,p-b))<=0;//做法一
// return !dcmp(Cro(p-a,b-a))&&dcmp(min(a.x,b.x)-p.x)<=0&&dcmp(p.x-max(a.x,b.x))<=0&&dcmp(min(a.y,b.y)-p.y)<=0&&dcmp(p.y-max(a.y,b.y))<=0;//做法二
//PA,AB共线且P在AB之间(其实也可以用len(p-a)+len(p-b)==len(a-b)判断,但是精度损失较大)
}
inline LD dis_PL(Point p,Point a,Point b){//【点P到线段AB距离】
if(a==b)return Len(p-a);//AB重合
Vector x=p-a,y=p-b,z=b-a;
if(dcmp(Dot(x,z))<0)return Len(x);//P距离A更近
if(dcmp(Dot(y,z))>0)return Len(y);//P距离B更近
return Abs(Cro(x,z)/Len(z));//面积除以底边长
}
/*2.【点与直线】*/
inline int pan_PL_(Point p,Point a,Point b){//【判断点P是否在直线AB上】
return !dcmp(Cro(p-a,b-a));//PA,AB共线
}
inline Point FootPoint(Point p,Point a,Point b){//【点P到直线AB的垂足】
Vector x=p-a,y=p-b,z=b-a;
LD len1=Dot(x,z)/Len(z),len2=-1.0*Dot(y,z)/Len(z);//分别计算AP,BP在AB,BA上的投影
return a+z*(len1/(len1+len2));//点A加上向量AF
}
inline Point Symmetry_PL(Point p,Point a,Point b){//【点P关于直线AB的对称点】
return p+(FootPoint(p,a,b)-p)*2;//将PF延长一倍即可
}
/*3.【线与线】*/
inline Point cross_LL(Point a,Point b,Point c,Point d){//【两直线AB,CD的交点】
Vector x=b-a,y=d-c,z=a-c;
return a+x*(Cro(y,z)/Cro(x,y));//点A加上向量AF
}
inline int pan_cross_L_L(Point a,Point b,Point c,Point d){//【判断直线AB与线段CD是否相交】
return pan_PL(cross_LL(a,b,c,d),c,d);//直线AB与直线CD的交点在线段CD上
}
inline int pan_cross_LL(Point a,Point b,Point c,Point d){//【判断两线段AB,CD是否相交】
LD c1=Cro(b-a,c-a),c2=Cro(b-a,d-a);
LD d1=Cro(d-c,a-c),d2=Cro(d-c,b-c);
return dcmp(c1)*dcmp(c2)<0&&dcmp(d1)*dcmp(d2)<0;//分别在两侧
}
/*4.【点与多边形】*/
inline int PIP(Point *P,Re n,Point a){//【射线法】判断点A是否在任意多边形Poly以内
Re cnt=0;LD tmp;
for(Re i=1;i<=n;++i){
Re j=i<n?i+1:1;
if(pan_PL(a,P[i],P[j]))return 2;//点在多边形上
if(a.y>=min(P[i].y,P[j].y)&&a.y<max(P[i].y,P[j].y))//纵坐标在该线段两端点之间
tmp=P[i].x+(a.y-P[i].y)/(P[j].y-P[i].y)*(P[j].x-P[i].x),cnt+=dcmp(tmp-a.x)>0;//交点在A右方
}
return cnt&1;//穿过奇数次则在多边形以内
}
inline int judge(Point a,Point L,Point R){//判断AL是否在AR右边
return dcmp(Cro(L-a,R-a))>0;//必须严格以内
}
inline int PIP_(Point *P,Re n,Point a){//【二分法】判断点A是否在凸多边形Poly以内
//点按逆时针给出
if(judge(P[1],a,P[2])||judge(P[1],P[n],a))return 0;//在P[1_2]或P[1_n]外
if(pan_PL(a,P[1],P[2])||pan_PL(a,P[1],P[n]))return 2;//在P[1_2]或P[1_n]上
Re l=2,r=n-1;
while(l<r){//二分找到一个位置pos使得P[1]_A在P[1_pos],P[1_(pos+1)]之间
Re mid=l+r+1>>1;
if(judge(P[1],P[mid],a))l=mid;
else r=mid-1;
}
if(judge(P[l],a,P[l+1]))return 0;//在P[pos_(pos+1)]外
if(pan_PL(a,P[l],P[l+1]))return 2;//在P[pos_(pos+1)]上
return 1;
}
/*5.【线与多边形】*/
/*6.【多边形与多边形】*/
inline int judge_PP(Point *A,Re n,Point *B,Re m){//【判断多边形A与多边形B是否相离】
for(Re i1=1;i1<=n;++i1){
Re j1=i1<n?i1+1:1;
for(Re i2=1;i2<=m;++i2){
Re j2=i2<m?i2+1:1;
if(pan_cross_LL(A[i1],A[j1],B[i2],B[j2]))return 0;//两线段相交
if(PIP(B,m,A[i1])||PIP(A,n,B[i2]))return 0;//点包含在内
}
}
return 1;
}
/*五:【图形面积】*/
/*1.【任意多边形面积】*/
inline LD PolyArea(Point *P,Re n){//【任意多边形P的面积】
LD S=0;
for(Re i=1;i<=n;++i)S+=Cro(P[i],P[i<n?i+1:1]);
return S/2.0;
}
/*2.【圆的面积并】*/
/*3.【三角形面积并】*/
/*六:【凸包】*/
/*1.【求凸包】*/
inline bool cmp1(Vector a,Vector b){return a.x==b.x?a.y<b.y:a.x<b.x;};//按坐标排序
inline int ConvexHull(Point *P,Re n,Point *cp){//【Graham扫描法】求凸包
sort(P+1,P+n+1,cmp1);
Re t=0;
for(Re i=1;i<=n;++i){//下凸包
while(t>1&&dcmp(Cro(cp[t]-cp[t-1],P[i]-cp[t-1]))<=0)--t;
cp[++t]=P[i];
}
Re St=t;
for(Re i=n-1;i>=1;--i){//上凸包
while(t>St&&dcmp(Cro(cp[t]-cp[t-1],P[i]-cp[t-1]))<=0)--t;
cp[++t]=P[i];
}
return --t;//要减一
}
/*2.【旋转卡壳】*/
/*3.【半平面交】*/
struct Line{
Point a,b;LD k;Line(Point A=Point(0,0),Point B=Point(0,0)){a=A,b=B,k=atan2(b.y-a.y,b.x-a.x);}
inline bool operator<(const Line &O)const{return dcmp(k-O.k)?dcmp(k-O.k)<0:judge(O.a,O.b,a);}//如果角度相等则取左边的
}L[N],Q[N];
inline Point cross(Line L1,Line L2){return cross_LL(L1.a,L1.b,L2.a,L2.b);}//获取直线L1,L2的交点
inline int judge(Line L,Point a){return dcmp(Cro(a-L.a,L.b-L.a))>0;}//判断点a是否在直线L的右边
inline int halfcut(Line *L,Re n,Point *P){//【半平面交】
sort(L+1,L+n+1);Re m=n;n=0;
for(Re i=1;i<=m;++i)if(i==1||dcmp(L[i].k-L[i-1].k))L[++n]=L[i];
Re h=1,t=0;
for(Re i=1;i<=n;++i){
while(h<t&&judge(L[i],cross(Q[t],Q[t-1])))--t;//当队尾两个直线交点不是在直线L[i]上或者左边时就出队
while(h<t&&judge(L[i],cross(Q[h],Q[h+1])))++h;//当队头两个直线交点不是在直线L[i]上或者左边时就出队
Q[++t]=L[i];
}
while(h<t&&judge(Q[h],cross(Q[t],Q[t-1])))--t;
while(h<t&&judge(Q[t],cross(Q[h],Q[h+1])))++h;
n=0;
for(Re i=h;i<=t;++i)P[++n]=cross(Q[i],Q[i<t?i+1:h]);
return n;
}
/*4.【闵可夫斯基和】*/
Vector V1[N],V2[N];
inline int Mincowski(Point *P1,Re n,Point *P2,Re m,Vector *V){//【闵可夫斯基和】求两个凸包{P1},{P2}的向量集合{V}={P1+P2}构成的凸包
for(Re i=1;i<=n;++i)V1[i]=P1[i<n?i+1:1]-P1[i];
for(Re i=1;i<=m;++i)V2[i]=P2[i<m?i+1:1]-P2[i];
Re t=0,i=1,j=1;V[++t]=P1[1]+P2[1];
while(i<=n&&j<=m)++t,V[t]=V[t-1]+(dcmp(Cro(V1[i],V2[j]))>0?V1[i++]:V2[j++]);
while(i<=n)++t,V[t]=V[t-1]+V1[i++];
while(j<=m)++t,V[t]=V[t-1]+V2[j++];
return t;
}
/*5.【动态凸包】*/
/*七:【圆】*/
/*1.【三点确定一圆】*/
#define S(a) ((a)*(a))
struct Circle{Point O;LD r;Circle(Point P,LD R=0){O=P,r=R;}};
inline Circle getCircle(Point A,Point B,Point C){//【三点确定一圆】暴力解方程
LD x1=A.x,y1=A.y,x2=B.x,y2=B.y,x3=C.x,y3=C.y;
LD D=((S(x2)+S(y2)-S(x3)-S(y3))*(y1-y2)-(S(x1)+S(y1)-S(x2)-S(y2))*(y2-y3))/((x1-x2)*(y2-y3)-(x2-x3)*(y1-y2));
LD E=(S(x1)+S(y1)-S(x2)-S(y2)+D*(x1-x2))/(y2-y1);
LD F=-(S(x1)+S(y1)+D*x1+E*y1);
return Circle(Point(-D/2.0,-E/2.0),sqrt((S(D)+S(E)-4.0*F)/4.0));
}
inline Circle getcircle(Point A,Point B,Point C){//【三点确定一圆】向量垂心法
Point P1=(A+B)*0.5,P2=(A+C)*0.5;
Point O=cross_LL(P1,P1+Normal(B-A),P2,P2+Normal(C-A));
return Circle(O,Len(A-O));
}
/*2.【最小覆盖圆】*/
inline int PIC(Circle C,Point a){return dcmp(Len(a-C.O)-C.r)<=0;}//判断点A是否在圆C内
inline void Random(Point *P,Re n){for(Re i=1;i<=n;++i)swap(P[i],P[rand()%n+1]);}//随机一个排列
inline Circle Min_Circle(Point *P,Re n){//【求点集P的最小覆盖圆】
// random_shuffle(P+1,P+n+1);
Random(P,n);Circle C=Circle(P[1],0);
for(Re i=2;i<=n;++i)if(!PIC(C,P[i])){
C=Circle(P[i],0);
for(Re j=1;j<i;++j)if(!PIC(C,P[j])){
C.O=(P[i]+P[j])*0.5,C.r=Len(P[j]-C.O);
for(Re k=1;k<j;++k)if(!PIC(C,P[k]))C=getcircle(P[i],P[j],P[k]);
}
}
return C;
}
/*3.【三角剖分】*/
inline LD calc(Point A,Point B,Point O,LD R){//【三角剖分】
if(A==O||B==O)return 0;
Re op=dcmp(Cro(A-O,B-O))>0?1:-1;LD ans=0;
Vector x=A-O,y=B-O;
Re flag1=dcmp(Len(x)-R)>0,flag2=dcmp(Len(y)-R)>0;
if(!flag1&&!flag2)ans=Abs(Cro(A-O,B-O))/2.0;//两个点都在里面
else if(flag1&&flag2){//两个点都在外面
if(dcmp(dis_PL(O,A,B)-R)>=0)ans=R*R*Angle(x,y)/2.0;//完全包含了圆弧
else{//分三段处理 △+圆弧+△
if(dcmp(Cro(A-O,B-O))>0)swap(A,B);//把A换到左边
Point F=FootPoint(O,A,B);LD lenx=Len(F-O),len=sqrt(R*R-lenx*lenx);
Vector z=turn_P(F-O,Pi/2.0)*(len/lenx);Point B_=F+z,A_=F-z;
ans=R*R*(Angle(A-O,A_-O)+Angle(B-O,B_-O))/2.0+Cro(B_-O,A_-O)/2.0;
}
}
else{//一个点在里面,一个点在外面
if(flag1)swap(A,B);//使A为里面的点,B为外面的点
Point F=FootPoint(O,A,B);LD lenx=Len(F-O),len=sqrt(R*R-lenx*lenx);
Vector z=turn_P(F-O,Pi/2.0)*(len/lenx);Point C=dcmp(Cro(A-O,B-O))>0?F-z:F+z;
ans=Abs(Cro(A-O,C-O))/2.0+R*R*Angle(C-O,B-O)/2.0;
}
return ans*op;
}
int main(){}