数论六之计算几何干货——计算几何模板解释全集 及 模板检验训练场

点和向量及运算

点：用坐标表示，$P(x,y)$

向量：有大小，方向，与二维点坐标一一对应，$P=(x,y)$

向量的运算

定义向量$A=(x_1,y_1),B=(x_2,y_2)$

• 加减：$A\pm B=(x_1\pm x_2,y1\pm y_2)$

• 数乘：$kA=(k{x}_1,k{y}_1)$

• 模长：$|A|=\sqrt{x_1^2+y_1^2}$

• 角度：$\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}2(y_1,x_1)$（返回的是弧度制）

  $\tan \theta =\frac{y}{x}\to \mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}2(y,x)=\theta$

• 点积/内积($\cdot$)：$A\cdot B=x_1{x}_2+y_1{y}_2$

  几何意义：$A$乘以$B$在$A$上的投影，$A\cdot B=|A||B|\cos \theta$

  若$A\cdot B=0$则两向量垂直；$>0$夹角为锐角；$<0$为钝角

  

• 叉积/外积($\times$)：$A\times B=\left|\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\end{array}\right|=x_1{y}_2-x_2{y}_1=-B\times A$

  与行列式的几何意义联系：都表示$xoy$二维平面上以$A,B$为邻边的平行四边形的有向面积，$A\times B=|A||B|\sin \theta$

  若$A\times B=0$，则向量$A,B$共线；$>0$说明从$B$转到$A$是顺时针；$<0$则从$B$转到$A$是逆时针

  

• 旋转：把$A$逆时针旋转$\alpha$弧度得到新向量$A^{\prime }$
  $$A=(x_1,y_1)=r(\cos \theta ,\sin \theta )\Rightarrow A^{\prime }=r(\cos (\theta +\alpha ),\sin (\theta +\alpha ))$$$$=r(\cos \theta \cos \alpha -\sin \theta \sin \alpha ,\cos \theta \sin \alpha +\sin \theta \cos \alpha )$$$$=(r\cos \theta \cos \alpha -r\sin \theta \sin \alpha ,r\cos \theta \sin \alpha +r\sin \theta \cos \alpha )$$$$=(x_1\cos \alpha -y_1\sin \alpha ,x_1\sin \alpha +y_1\cos \alpha )$$

$\mathrm{acos}$接受$[-1,1]$的浮点数，输出的范围是$[0,\pi ]$
$\mathrm{asin}$接受$[-1,1]$的浮点数，输出的范围是$[-\frac{\pi }{2},\frac{\pi }{2}]$
$\mathrm{atan}$接受浮点数，输出的范围是$[-\frac{\pi }{2},\frac{\pi }{2}]$
$\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}2$接受两个浮点数$y,x$（坐标），输出的范围是$[-\pi ,\pi ]$
点在第一象限$[0,\frac{\pi }{2}]$，第二象限$[\frac{\pi }{2},\pi ]$，第三象限$[-\pi ,-\frac{\pi }{2}]$，第四象限$[-\frac{\pi }{2},0]$

struct vec {
	double x, y;
	vec(){}
	vec( double X, double Y ) { x = X, y = Y; }
	vec operator + ( vec t ) { return vec( x + t.x, y + t.y ); }
	vec operator - ( vec t ) { return vec( x - t.x, y - t.y ); }
	vec operator * ( double t ) { return vec( x * t, y * t ); }
	double dot( vec t ) { return x * t.x + y * t.y; }
	friend double dot( vec s, vec t ) { return s.x * t.x + s.y * t.y; }
	double cross( vec t ) { return x * t.y - y * t.x; }
	friend double cross( vec s, vec t ) { return s.x * t.y - s.y * t.x; }
	double angle() { return atan2( y, x ); }
	friend double angle( vec t ) { return atan2( t.y, t.x ); }
	vec rotate( double rad ) { return vec( x * cos( rad ) - y * sin( rad ), x * sin( rad ) + y * cos( rad ) ); }
	friend vec rotate( vec t, double rad ) { return vec( t.x * cos( rad ) - t.y * sin( rad ), t.x * sin( rad ) + t.y * cos( rad ) ); }
	double len() { return sqrt( dot( *this ) ); }
};

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	point operator + ( vec p ) { return point( x + p.x, y + p.y ); }
	vec operator - ( point v ) { return vec( x - v.x, y - v.y ); }
};

直线和线段

直线的表示方式

• 斜率式：$kx+b$

  优点：可以控制$x$的取值来得到线段或者射线

  缺点：不能表示$x=c$的直线，计算斜率$k$需要除法

• 一般式：$ax+by+c=0$

  优点：可以表示所有的直线，$a,b,c$不是固定的，可以同时扩大缩小

  尽量避免除法以提高精度，方便用解析几何的方法算东西

  缺点：参数会多一个，线段或者射线不好处理

• 两点式：$A,B$

  优点：简单

  缺点：没有考虑直线上的点的特点

• 一个点和直线的方向向量
  设向量$v=B-A$,则直线$AB$上的点可以表示为$A+vt$（$t$任取）

  优点：可以表示所有的直线，并且通过限制$t$的取值，可以得到线段或者射线

  如取$t$在$0-1$之间，得到线段，取$t>0$，得到射线

求解点到直线的距离/点在直线上

• 用解析几何的方法设直线$l:ax+by+c=0$，点$p(x_0,y_0)$到$l$的距离为$\frac{|a{x}_0+b{y}_0+c|}{\sqrt{a^2+b^2}}$

• 用叉积的方法

  用叉积求出以$p_2{p}_1$和$p_2{p}_0$为邻边的平行四边形的面积，然后除以底$p_1{p}_2$的长度

  

  同样通过叉积判断点$p$是否在直线上，如果在那么构成的平行四边形就是一条线面积为$0$

struct line {//直线由一个端点p加上表示方向的向量v表示 
	point p; vec v;
	line(){}
	line( point P, vec V ) { p = P, v = V; };
	friend double dis( line l, point p ) {//点p到直线l的距离
		vec v = p - l.p;
		double s = cross( l.v, v );//叉积求出平行四边形的面积 
		return fabs( s ) / l.v.len();//面积除以底边的长度就是高(距离) 
	}
	friend bool online( line l, point p ) {//判断点p是否在直线l上 => 共线叉积面积为0 
		vec v = p - l.p;
		double s = cross( l.v, v );
		return dcmp( s ) == 0;
	}
};

求解点到线段的距离/点在线段上

• 做垂线，垂足落在线段上，那么点到线段的距离就是点到线段所在直线的距离

• 垂足不落在线段上，此时点到线段的距离就是点到线段的两个端点的距离中较小的那个

$\angle p_0{p}_2{p}_1$是钝角且$\angle p_0{p}_1{p}_2$是锐角，则$p_0$到线段$p_{1}p+2$的距离是线段$p_0{p}_2$的长度

用点积可以求得夹角是否为钝角

$\angle p_0{p}_1{p}_2$是钝角且$\angle p_0{p}_2{p}_1$是锐角，则$p_0$到线段$p_1{p}_2$的距离是线段$p_0{p}_1$的长度

所以必须两个夹角都是锐角才能说明垂足在线段上

double dis( point p1, point p2, point p ) {//点p到线段(p1,p2)的距离 
	double d1 = dot( p - p1, p2 - p1 ), d2 = dot( p - p2, p1 - p2 );
	if( dcmp( d1 ) >= 0 && dcmp( d2 ) >= 0 ) {
		line l( p1, p2 - p1 );
		return dis( l, p );
	}
	else if( dcmp( d1 ) < 0 ) return len( p - p1 );
	else return len( p - p2 );
}

判断点$p$是否在线段上，首先要保证$p$在线段所在直线上，然后出现在两个端点中间

用点积判断，在线段中间则与两端点形成的向量应该是反向的，点积结果应为负

bool onseg( point p1, point p2, point p ) {//判断点p是否在由点p1和点p2构成的线段上 
	line l( p1, p2 - p1 );
	if( ! online( l, p ) ) return 0;//首先要先能在ab直线上
	//在线段ab中间->与两端点形成的向量应该是反向的->点积为负 
	if( dcmp( dot( p - p1, p - p2 ) ) <= 0 ) return 1;
	else return 0; 
}

求解两条线段是否相交

跨立实验

判断$p_1,p_2$是否在直线$p_3{p}_4$的两侧；再判断$p3,p4$是否在直线$p_1{p}_2$的两侧

跨立实验通过叉积求其正负判断出左边还是右边

两点分别引出与直线端点的两向量的叉积应为非负数

可能出现某条线段端点在另一条线段上，这种情况下的叉积是$0$

但是也要排除两条线段共线不交的叉积为$0$的情况

bool seg_intersect( point p1, point p2, point p3, point p4 ) {
//判断线段(p1,p2)和(p3,p4)是否有交->跨立实验->一线段两端点分居另一线段两方
	double s1 = cross( p2 - p1, p3 - p1 ), s2 = cross( p2 - p1, p4 - p1 );
	double s3 = cross( p4 - p3, p1 - p3 ), s4 = cross( p4 - p3, p2 - p3 );
	int ret1 = dcmp( s1 ) * dcmp( s2 ), ret2 = dcmp( s3 ) * dcmp( s4 );
	if( ret1 > 0 || ret2 > 0 ) return 0;
	if( ret1 == 0 && ret2 == 0 ) {
		if( onseg( p1, p2, p3 ) || onseg( p1, p2, p4 ) ) return 1;
        //可能交点是两线段的某个端点 
		else return 0;//共线不交
	}
	else return 1;	
}

求解两直线的交点

定义两直线$l:p_1+t(p_2-p_1)r:p_3+s(p_4-p_3)$

则交点可列等式$p_1+t(p_2-p_1)=p_3+s(p_4-p_3)$
$$\begin{gathered} \Rightarrow t(p_2-p_1)-s(p_4-p_3)=p_3-p_1 \\ \Rightarrow \{\begin{array}{ll}x:& t(x_2-x_1)-s(x_4-x_3)=x_3-x_1\\ y:& t(y_2-y_1)-s(y_4-y_3)=y_3-y_1\end{array} \end{gathered}$$

$$D_0=\left|\begin{array}{cc}{x}_2-x_1& -(x_4-x_3)\\ y_2-y_1& -(y_4-y_3)\end{array}\right|=-\left|\begin{array}{cc}{x}_2-x_1& x_4-x_3\\ y_2-y_1& y_4-y_3\end{array}\right|$$

$$D_t=\left|\begin{array}{cc}{x}_3-x_1& -(x_4-x_3)\\ y_3-y_1& -(y_4-y_3)\end{array}\right|=\left|\begin{array}{cc}{x}_1-x_3& x_4-x_3\\ y_1-y_3& y_4-y_3\end{array}\right|$$

$$D_s=\left|\begin{array}{cc}{x}_2-x_1& x_3-x_1\\ y_2-y_1& y_3-y_1\end{array}\right|=-\left|\begin{array}{cc}{x}_2-x_1& x_1-x_3\\ y_2-y_1& y_1-y_3\end{array}\right|$$

$$\begin{gathered} \vec{u}=p_2-p_1=(x_2-x_1,y_2-y_1) \\ \vec{v}=p_4-p_3=(x_4-x_3,y_4-y_3) \\ \vec{w}=p_1-p_3=(x_1-x_3,y_1-y_3) \end{gathered}$$

$$\{\begin{array}{l}{D}_0=-\vec{u}\times \vec{v}\\ D_s=-\vec{u}\times \vec{w}\\ D_t=\vec{w}\times \vec{v}\end{array}$$

$$\begin{gathered} t=\frac{D_t}{D_0}=\frac{\vec{v}\times \vec{w}}{\vec{u}\times \vec{v}} \\ s=\frac{D_s}{D_0}=\frac{\vec{u}\times \vec{w}}{\vec{u}\times \vec{v}} \end{gathered}$$

注意当且仅当$D_0=0$时，方程可能有多解（直线重合）/ 可能无解（直线平行）

否则可以根据$t$和$s$的取值判断交点是否在线段上

point line_intersect( line l, line r ) {//求交点
	vec v = l.p - r.p;//\vec{w}
	double s = cross( l.v, r.v );//l.v->\vec{u};r.v->\vec{v}
	if( dcmp( s ) != 0 ) {//两直线不重合
	 	double t = cross( r.v, v ) / s;
	 	return l.p + l.v * t;
	}
	else return l.p;
}

多边形

用多边形的顶点集合表示（顶点顺序不同，对应的多边形可能会不同

也可以按极角的规律顺序

求解多边形面积

• 凸多边形

  

  对凸多边形做三角剖分，然后求每个三角形的面积

  三角形面积：两个向量的叉积表示夹成的平行四边形的面积，除以$2$就是三角形的面

• 凹多边形

  • 计算出$△ABC$的面积

    

  • 计算出凹四边形$ABCD$的面积（实际上是$△ACD-△ABC$）

    

  • 计算出复杂五边形$ABCDE$的面积（实际上是四边形$ABCD-△ADE$）

    

    $PS$：这个五边形不是简单多边形，简单多边形是由边不交的线段组成的多边形，这里$AE$和$BC$相交了，不妨设交点是$G$

    那么这里算的面积实际上是凹四边形$GCDE$的面积减去$△AGB$的面积

  • 计算出凹六边形$ABCDEF$的面积（实际上是五边形$ABCDE+△AEF$）

    

  实际上，这种剖分方法已经不再适用于凹多边形了

其实只需要搞清楚什么时候加三角形面积什么时候减三角形面积

（发现$AC$在$AB$右边所以减，$AD$在$AC$左边所以加$...$）

用叉积算出有向面积求和即可，简单多边形的面积就是多个三角形的有向面积和

甚至可以不用多边形上的顶点，随便选一个点进行三角剖分都可以

double area( vector < point > polygon ) {//计算多边形面积 多边形由所有端点集合表示
	double s = 0; int n = polygon.size();
	for( int i = 0;i < n;i ++ )
		s += ( polygon[i] - point( 0, 0 ) ).cross( polygon[( i + 1 ) % n] - point( 0, 0 ) );
	return fabs( s / 2 ); 
}

求解多边形重心

重心公式$(\frac{x_1+x_2+x_3}{3},\frac{y_1+y_2+y_3}{3})$

多边形的重心并不是顶点坐标平均数，当然如果所有质量都在顶点上那么这个公式就是正确的

但是显然多边形是质量分布均匀的

同样的，把多边形进行三角剖分，把每个三角形的质量放到三角形重心处（大小为三角形有向面积），求加权平均数

point center( vector < point > polygon ) {
	double s = 0; int n = polygon.size(); point p( 0, 0 );
	for( int i = 0;i < n;i ++ ) {
		double t = ( polygon[i] - point( 0, 0 ) ).cross( polygon[( i + 1 ) % n] - point( 0, 0 ) );
		s += t;
		p.x += ( polygon[i].x + polygon[( i + 1 ) % n].x ) * t;
		p.y += ( polygon[i].y + polygon[( i + 1 ) % n].y ) * t;
	}
	p.x /= 3, p.x /= s, p.y /= 3, p.y /= s;
	return p;
}

求解判断定点与多边形的位置关系

点恰好在多边形边上，可以通过枚举多边形的边特判掉；剩下来只需要考虑位置在内外

• 射线法

  从定点随机射出一条引向无穷远点的射线，根据射线与多边形的交点个数判断位置

  如果交点为偶数，则在多边形外；否则在多边形内

  证明显然：如果是在外面那么射入多边形就必定伴随着射出多边形（多边形是个封闭图形）产生偶数个交点

  

bool inpolygon_line( point p, vector < point > polygon ) {//射线法判断点p是否在多边形内 
  	int n = polygon.size();
  	double r = ( rand() / double( RAND_MAX ) - 0.5 ) * 2 * M_PI;//随机一条射线 
  	vec v( cos( r ), sin( r ) );
  	bool ans = 0;
  	for( int i = 0;i < n;i ++ ) {
  		if( onseg( polygon[i], polygon[( i + 1 ) % n], p ) ) return 1;//特判点在多边形边上
  		if( inangle( polygon[i] - p, polygon[( i + 1 ) % n] - p, v ) ) ans ^= 1;
  	}
  	return ans; 
}

• 转角法

  沿多边形走一圈，累计绕点旋转了多少角度

  $0$：点在多边形外

  $\pm \pi$：点在多边形上

  $\pm 2\pi$：点在多边形内

  

  

  缺点：角度计算要用反三角函数，精度不高，速度较慢

  然则在实现中，不可能花费时间精度去真的求角度判定

  非常巧妙得利用多边形定点编号自带顺序分辨定点的位置特点（结合代码观图解）

  

  • 点$I$，点$L$

    所有点都在其下方/上方，根据$y$永远不可能触发d<=0&&d'>0的条件，程序判定正确

  • 点$K$，点$J$

    因为多边形自带顶点顺序，满足$c$条件时，$d$条件恰好与$if$语句相反，程序判定正确

  • 点$M$

    $\overrightarrow{CD}$，满足$i+1$在$i$左边且$i+1$在点下方，$i$在点上方，cnt--

    $\overrightarrow{FG}$，满足$i+1$在$i$右边且$i+1$在点上方，$i$在点下方，cnt++

    程序判定正确

  关键点在于，因为多边形自带顶点顺序且封闭，所以当点在多边形外的时候，如果点(竖直位置)在一条边的中间，那么必定会有另外一条反向边含有点

  有可能多边形是锯齿图/波浪图，多次触发，也一定是成对的

  bool inpolygon_rad( point p, vector < point > polygon ) {
  	int n = polygon.size(), cnt = 0;
  	for( int i = 0;i < n;i ++ ) {
  		if( onseg( polygon[i], polygon[( i + 1 ) % n], p ) ) return 1;
  		double c = cross( polygon[i] - p, polygon[( i + 1 ) % n] - p );
  		double d1 = polygon[i].y - p.y, d2 = polygon[( i + 1 ) % n].y - p.y;
  		if( c < 0 && d1 <= 0 && d2 > 0 ) cnt ++;
  		if( c > 0 && d2 <= 0 && d1 > 0 ) cnt --;
          //c>0:p与i+1构成的线在p与i构成的线左边(需要p与i的线逆时针旋转)
  	}
  	return cnt != 0;
  }
  

• 改进的转角法

  把被测点看成坐标原点，对多边形的一条边$P_i-P_{i+1}$，有四种情况

  • $P_{i+1}$在$P_i$的相同象限，弧度不变

  • $P_{i+1}$​在$P_i$​的下一象限，弧度加$\frac{\pi }{2}$

  • $P_{i+1}$在$P_i$的上一象限，弧度减$\frac{\pi }{2}$

  • $P_{i+1}$​​在$P_i$​​的相对象限，计算向量$\overrightarrow{O{P}_i}$​​和$\overrightarrow{O{P}_{i+1}}$​的叉积$MS$​

    • $MS=0$​，点在多边形上
    • $MS>0$​​，弧度加$\pi$
    • $MS<0$​，弧度减$\pi$

  实现时，需特判点在边上的情况，判断象限的时候$0$​​​​看成正数

  原始转交法的几何原理就是此时的跨象限

  bool inpolygon_Rad( point p, vector < point > polygon ) {
  	int n = polygon.size(), cnt = 0;
  	vector < vec > v;
  	for( int i = 0;i < n;i ++ )
  		v.push_back( polygon[i] - p );
  	for( int i = 0;i < n;i ++ ) {
  		vec v1 = v[i], v2 = v[( i + 1 ) % n ];
  		double c = cross( v1, v2 );
  		if( c == 0 && dot( v1, v2 ) <= 0 ) return 1;
  		int d1 = quadrant( v1 ), d2 = quadrant( v2 );
  		if( d2 == ( d1 + 1 ) % 4 ) cnt ++;
  		if( d2 == ( d1 + 3 ) % 4 ) cnt --;
  		if( d2 == ( d1 + 2 ) % 4 ) {
  			if( c > 0 ) cnt += 2;
  			else if( c < 0 ) cnt -= 2;
  		}
  	}
  	return cnt;
  }
  

• 对于凸多边形

  还可以利用叉积判断定点逆时针绕多边形转一圈都在边的左边，则定点在多边形内

凸包

定义：将最外层的点连接起来构成的凸多边形，能包含点集中所有的点

形象上理解就是，用橡皮筋撑到无穷大，然后放手，最后橡皮筋蹦到挂住的点围成的多边形

性质

• 凸包内任意两点构成的线段都被包含于凸包内
• 凸包是凸多边形
• 凸包的交仍是凸包
• $n$个点形成的凸包大小，在完全随机的情况下期望不大

graham扫描法

类似于斜率优化，维护凸壳

• 先找出算法开始点$p_0$——左下点（先找最下面的点，若有多个最下面点则选择最左的一个点）

• 然后对每个点进行极角排序（极角相同的距离点近的在前面

  极角：点与极轴的夹角

  极轴：过原点的水平线（$x$​轴）

  

• 先把$p_0,p_1,p_2$入栈

  

  把点$C,D,I$入栈

• 按极角排序枚举点

  • 先枚举$G$​

    

    看$D-I-G$发现$I$左拐到$G$，把$G$入栈

    （如果$I$右拐或者不变方向到$G$就把$I$弹栈）

    在左边还是右边的判断，通过叉积可得

  • 然后枚举$F$

    

    $I-G-F$不转向，将$G$弹栈；$D-I-F$左转，不弹$I$，停止；将$F$入栈

  • 接着枚举$H$​

    $I-F-H$右转，将$F$弹栈；$D-I-H$左转，保留停止；将$H$入栈

    

  • 再枚举$E$

    $I-H-E$左转，将$E$入栈

    

  • 考虑$A$

    $H-E-A$右转，弹$E$；$I-H-A$左转，保留停止；入栈$A$

    

  • 最后是$B$

    $H-A-B$左转，入栈，整个程序结束

    

时间复杂度$O(n\log n)$​，瓶颈在排序

point Po[maxn], sta[maxn], p[maxn];
int graham() {
	for( int i = 0;i < n;i ++ ) p[i] = Po[i];
	int o = 0;
	for( int i = 1;i < n;i ++ )
		if( p[i].y < p[o].y || ( p[i].y == p[o].y && p[i].x < p[o].x ) ) o = i;
	swap( p[o], p[0] );
    //p[0]是起点 不参与排序 所以从p+1开始
	sort( p + 1, p + n, []( point s, point t ) {
		vec v1 = s - p[0], v2 = t - p[0];
		double c = cross( v1, v2 );
		return dcmp( c ) > 0 || ( dcmp( c ) == 0 && dcmp( v1.dot() - v2.dot() ) < 0 );
	} );
	sta[0] = p[0], sta[1] = p[1]; int top = 1;
	for( int i = 2;i < n;i ++ ) {
		while( top && dcmp( cross( sta[top] - sta[top - 1], p[i] - sta[top - 1] ) ) <= 0 ) top --;
		sta[++ top] = p[i];
	}
	return top + 1;
}

缺点：如果要保留凸包边上的点，则刚才的极角排序就不适用了

如果只有在向右拐的时候才弹栈，直走的时候不弹栈，也不能保留凸包边上所有的点

e.g.

• 先$A,B,C$​入栈
• 然后$D$​入栈
• 此时遇到$AF,AE$​极角相同，按距离先排$F$​
• $C-D-F$​左转，$F$​入栈
• $D-F-E$​右转，结果最后$E$​​​被弹出栈了

如果规定极角相同离远点小，则$F$​活$B$​​亡，总归是顾此失彼

graham扫描法加强版

• 把对所有点进行的极角排序改造（按$x$排序）

• 前两个点$p_0,p_1$入栈

• 枚举后面每个点

  如果栈内元素个数$>1$，设栈顶$top$，栈顶下一个$lst$

  考虑$lst-top-p_i$是否向左拐

  右拐或者直走都弹栈顶，直到向左拐或只剩下一个栈顶，将$p_i$​入栈

• 操作后得到点集的下凸壳

• 把极角排序后的点倒过来做一遍得到点集的上凸壳

e.g.

• 所有点排序

• 先把点$1,2$入栈

  

• $2$弹栈，$3$进栈

  

• $4$进栈

  

• $4$出栈，$5$进栈

  

• $5$出栈，$6$进栈

  

• $6$出栈，$7$进栈

  

• $8$进栈

• $8$出栈，$9$​进栈

• $9$​​出栈，$10$​​​进栈

• $11$进栈

• $11$​​出栈，$12$​​​进栈

• $12$出栈，$13$进栈

• $14$进栈

• 得到下凸壳

  

• 同理倒序枚举扫描找出上凸壳，凸包刻画成功

  

实际上以$y$​作为第一关键字排序，先求右半凸包，再求左半凸包也一样

如果叉积为$0$，即不改变方向时不弹栈，就可以保留凸包边上的点

point Po[maxn], p[maxn], sta[maxn];
int GraHam( int n ) {
	for( int i = 0;i < n;i ++ ) p[i] = Po[i];
	sort( p, p + n, []( point s, point t ) { return s.x == t.x ? s.y < t.y : s.x < t.x; } );
	int top = 0;
	for( int i = 0;i < n;i ++ ) {
		while( top > 1 && cross( sta[top - 1] - sta[top - 2], p[i] - sta[top - 2] ) <= 0 )
			top --;
		sta[top ++] = p[i];
	}
	int t = top;
	for( int i = n - 2;~ i;i -- ) {
		while( top > t && cross( sta[top - 1] - sta[top - 2], p[i] - sta[top - 2] ) <= 0 )
			top --;
		sta[top ++] = p[i];
	}
	if( n > 1 ) return top - 1;//倒着求凸包的时候多放了个1点进去
	else return top;
}

更重要的区别在于，加强版是根据坐标排序，而极角排序显然会带来更大的误差

所以还是选择维护上下凸壳的加强版更优（亲身体验：UVA-10256 the great divide）

圆

圆的表示：圆心$C$和半径$r$

圆的一般方程形式：$(x-x_0{)}^2+(y-y_0{)}^2=r^2$​

已知圆心角表示圆上一点：$(x_0+r\cos \theta ,y_0+r\sin \theta )$

求解圆与直线的交点

判定相交

算出圆心到直线的距离$d$

• $d=r$，一个交点
• $d>r$，不相交
• $d<r$，两个交点

求解交点

算出垂足$D$​（已知直线$AB$可得$AD$直线方向向量，从$A$点走方向向量长度$d$

$|CD|=|BD|=\sqrt{r^2-d^2}$

利用$BC$直线的方向向量同理找到所有交点坐标

求解圆与圆的交点

判定相交：圆心距$d\le r_1+r_2$

求解交点

• 联立方程求解

• 几何法

  已知$|AC|,|AB|,|BC|$余弦定理可求$\angle BCA=\alpha$

  $CA$​直线已知，可求得$CA$​直线的单位方向向量$\overrightarrow{ca}$​

  将$\overrightarrow{ca}$逆时针/顺时针旋转$\alpha$，从$C$开始走长度$|BC|$的向量

求解定点与定圆的切线

切线数量

• 点在圆内，无
• 点在圆上，$1$条
• 点在圆外，$2$条

求解切线

• 数学解析法

  点$(x_0,y_0)$​​​在圆$(x-a{)}^2-(y-b{)}^2=r^2$​​​的切线解析式为$(x_0-a{)}^2+(y_0-b{)}^2=r^2$​​

• 几何法

  

  已知$|AC|=d,|AB|=r$，勾股定理求得$|BC|$长度

  余弦定理可解$\angle BCA\alpha$，又知直线$CA$方向向量

  逆时针/顺时针旋转$\alpha$得$CB$直线方向向量，从$C$出发走长度$|BC|$

求解圆公切线

公切线数量判定

• 重合：$d=0$，无数条
• 内含：$0<d<r_2-r_1$，$0$条
• 内切：$d=r_2-r_1$，$1$条
• 相交：$r_2-r_1<d<r_2+r_1$，$2$条
• 外切：$d=r_1+r_2$，$1$条
• 外离：$r_1+r_2<d$，$2$条外公切线，$2$条内公切线

求解公切线

• 重合，内含不用求

• 内切/外切的一条内/外公切线对应的是过唯一交点的切线

• 两条公切线

  

  已知$|AC|=d,|AD|=r_1,|CE|=r_2$，垂足为$D,E$

  相似求解$\frac{FA}{FC}=\frac{AD}{CE}=\frac{FD}{FE}$，$|AF|,|FD|,|DE|$长度可得

  余弦定理得$\angle AFD$

  $CA$直线已知，从$A$出发走$\overrightarrow{ca}$长度$|AF|$得到两外公切线交点$F$

  逆时针/顺时针旋转相同角度$\alpha$，走一定长度得到与圆的交点，两点确定一条直线

  （内切线同理）

pick定理

点阵中顶点在格点上的多边形面积公式$2S=2a+b-2$

其中$a$表示多边形内部的点数，$b$表示多边形边界上的点数，$S$表示多边形的面积

一般是反过来使用求$a$

$S$多边形面积可以用叉积解决

边界上的点数$b$用$\gcd$​解决：求和枚举多边形的边，一边两点的横纵坐标差的最大公约数$+1$

设线段$(x_1,y_1)-(x_2,y_2)\Rightarrow cnt=\gcd (|x_2-x_1|,|y_2-y_1|)+1$

pick定理扩展：一条线段下方有多少顶点

等价于求${\sum }_{x=0}^n\lfloor \frac{ax+b}{c}\rfloor$

求解圆相交面积

#include <bits/stdc++.h>
using namespace std;
#define eps 1e-10

double pi = acos( -1.0 );

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	friend double dis( point p1, point p2 ) { return sqrt( ( p1.x - p2.x ) * ( p1.x - p2.x ) + ( p1.y - p2.y ) * ( p1.y - p2.y ) ); }
};

struct circle {
	point o; double r;
	circle(){}
	circle( point O, double R ) { o = O, r = R; }
};

double area_circle( circle c1, circle c2 ) {
	double d = dis( c1.o, c2.o );
	if( dcmp( d - c1.r - c2.r ) >= 0 ) return 0; //外切或外离
	if( dcmp( d - fabs( c1.r - c2.r ) ) <= 0 ) { //内切或内含 
		double r = min( c1.r, c2.r );
		return pi * r * r;
	}
	//相交
	double rad1 = acos( ( c1.r * c1.r + d * d - c2.r * c2.r ) / ( 2 * c1.r * d ) ); 
	double rad2 = acos( ( c2.r * c2.r + d * d - c1.r * c1.r ) / ( 2 * c2.r * d ) );
	return rad1 * c1.r * c1.r + rad2 * c2.r * c2.r - c1.r * d * sin( rad1 );
}

欧拉定理(图论)

设$G$为任意联通的平面图，$n$表示顶点数，$m$表示边数，$r$表示面数（去除封闭面后剩下的“背景”也算完整的一面）

有$n-m+r=2$​恒成立

求解最小圆覆盖

平面上有一些点，用半径最小的圆覆盖这些点，用随机增量法

定义$C_i$​为覆盖前$i$​​个点的最小包围圆

• $p_i\in C_i-1$​，则$C_i=C_{i-1}$​

• $p_i\notin C_{i-1}$​时的处理，需要对$C_i$​​另外更新，要求$C_i$​​覆盖点集 { p 1 , p 2 , … , p i } \{p_1,p_2,…,p_i\} {p1​,p2​,…,pi​}​

  先让 C i = { p 1 , p i } C_i=\{p_1,p_i\} Ci​={p1​,pi​}，逐个判断点 { p 2 , p 3 , … , p i − 1 } \{p_2,p_3,…,p_{i-1}\} {p2​,p3​,…,pi−1​}

  如果存在$p_j\notin C_i$，则 C i = { p j , p i } C_i=\{p_j,p_i\} Ci​={pj​,pi​}

  同时，再依次判断点 { p 1 , p 2 , … , p j − 1 } \{p_1,p_2,…,p_{j-1}\} {p1​,p2​,…,pj−1​}​，如果存在$p_k\notin C_i$ ，则 C i = { p k , p j , p i } C_i=\{p_k,p_j,p_i\} Ci​={pk​,pj​,pi​}

  因为三点唯一确定一个圆，所以只需在此基础上判断其他的点是否位于圆内，不停地更新$p_k$​

  当$k$​​​​​的循环完成时，退出循环，更新$p_j$​​​​

  当$j$​​​的循环结束时，退出循环，更新$p_i$​​

  当$i=n$时，算法完成

从算法直观上来看，强烈意识感觉这是个假算法

可能会有这样的问题：$ABC$​​​​三点构成一个圆，现在加入$D$​​​点，假如$D$​​不在$◯ABC$​内

那么按照算法，$D$要在新构造的圆上

• 先枚举$A$​​​​，$AD$​​​构成圆的直径，假如$B$​​不在$◯AD$​内
• 让$BD$​​​​构成圆的直径
• 假如$A$​​​不在$◯BD$​​内，让$ABD$​构成一个圆
• 假如$C$​​​不在$◯ABD$​​内，让$CD$​构成圆的直径
• 假如$A$​​​不在$◯CD$​​内，让$ACD$​构成一个圆
• 假如$B$​​​不在$◯ACD$​​内，让$BCD$​构成一个圆

这里出现了直觉上的问题：有没有$A$​不在$◯BCD$内的可能

事实上是没有的

根据以上描述，平面上四个点$ABCD$，任意三个点作圆，另外一个点都不在作出的圆内

• 设$ABCD$​是凹四边形

  那么对这个凹四边形求凸包，得到三角形

  一定存在一个点在另外三个点构成的三角形内部（或者边上）

  点在三角形内，那么点自然在三角形的外接圆内

• 设$ABCD$​​是凸四边形

  假设四点共圆，那么对角和为$180^{\circ }$​

  如果某个点在此圆外

  由四边形的内角和为$360^{\circ }$，一组对角和小于$180^{\circ }$，另外一组对角和大于$180^{\circ }$​

  真正的圆必定是会被枚举到大于此圆的，此时又将某个点包含在圆内了

  所以极限边界就是四点共圆

必然存在一个点在另外三个点形成的圆内

因此算法流程走完后必定求得最小圆覆盖

时间复杂度

前$i$个点在包围$◯C_i$上也就$3$个点，所以点$i$在$◯C_i$上（不在$◯C_{i-1}$内）的概率是很低的，只有$\frac{3}{i}$

所以对于$i$，需要枚举$j$的概率是$\frac{3}{i}$

对于$j$，需要枚举k的概率也是$\frac{3}{j}$​

总的复杂度是线性的，当然这是建立在点随机排序的情况下，最坏还是$O(n^3)$​的

#define eps 1e-12

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	friend double len( point s, point t ) { return sqrt( ( s.x - t.x ) * ( s.x - t.x ) + ( s.y - t.y ) * ( s.y - t.y ) ); }
	point operator + ( point p ) { return point( x + p.x, y + p.y ); }
	point operator / ( double k ) { return point( x / k, y / k ); }
	double operator - ( point p ) { return x * x + y * y - p.x * p.x - p.y * p.y; }
}p[maxn];

double X( point p1, point p2 ) {
	return p1.x - p2.x;
}

double Y( point p1, point p2 ) {
	return p1.y - p2.y;
}

point get_c( point p1, point p2, point p3 ) {
	double t = 2 * ( Y( p2, p1 ) * X( p3, p1 ) - Y( p3, p1 ) * X( p2, p1 ) );
	return point( ( Y( p2, p1 ) * ( p3 - p1 ) - Y( p3, p1 ) * ( p2 - p1 ) ) / t, ( X( p3, p1 ) * ( p2 - p1 ) - X( p2, p1 ) * ( p3 - p1 ) ) / t );
}

int n; point c; double r;
void circle_cover() {
	random_shuffle( p + 1, p + n + 1 );
	c = p[1], r = 0;
	for( int i = 2;i <= n;i ++ )
		if( dcmp( len( c, p[i] ) - r ) > 0 ) {
			c = p[i], r = 0;
			for( int j = 1;j < i;j ++ )
				if( dcmp( len( c, p[j] ) - r ) > 0 ) {
					c = ( p[i] + p[j] ) / 2, r = len( c, p[i] );
					for( int k = 1;k < j;k ++ )
						if( dcmp( len( c, p[k] ) - r ) > 0 )
							c = get_c( p[i], p[j], p[k] ), r = len( c, p[i] );
				}
		}
}

求解最近点对

平面上有n个点，求最近点对

分治

• 对点按$x$坐标排序，作一条竖线$x=mid$​​​​​​​​​​，把点集分成两半，$L$​​​​​​​​​和$R$​​

  假设现在已经求出$L$​​​​​的最近点对距离$d_l$​​​​，$R$​​​的最近点对距离$d_r$​​

  最后的答案要么是$d=\min (d_l,d_r)$​，要么是一个$L$​中的点和一个$R$​中的点的最短距离

• 考虑$L$中的点和$R$中的点的最短距离如何求

  首先$L$​中的点，$x+d<mid$​的点肯定不行，把可行的点记作$L^{\prime }$

  同理$R$中的点，$x>mid+d$的点肯定不行，把可行的点记作$R^{\prime }$

  设$P=L^{\prime }\bigcup R^{\prime }$

  • 简单想法就暴力枚举$P$​中的点对，更新答案，但可能$P$​中的点数也很多，最坏还是$O(n^2)$​

  • 把$P$中的点按照$y$排序

    取定$P$​​中的点$p_i$​​，枚举到$p_j$​，使得$y_{p_{j+1}}-y_{p_i}\ge d$​，就不枚举了

    看上去貌似仍然是$O(n^2)$，事实上并不是，只需要枚举$O(n)$次

    因为$L$​​和$R$​中的任意两点距离均$\ge d$

    所以在重复点去掉的情况下， 对于任意$p_i$，满足$y_{p_j}-y_{p_i}<d$​​​​​的点只有常数个

    

闵可夫斯基和(Minkowski Sum)

给定点集$A,B$，它们的闵可夫斯基和$A+B$被定义为$\{a+b|a\in A,b\in B\}$​

可以理解为八其中一个图形沿着另外一个图形的边做平移得到的新图形

考虑两个凸包的闵可夫斯基和

由于结果也是一个凸包，所以结果集合的这些边也有极角序

并且结果集合上的边都来自于$A,B$两个凸包上的边

相当于对$A,B$两个凸包上的边根据极角序做归并排序，就是闽可夫斯基和

先对$A,B$两个凸包按照极角序排序

一般来说$A,B$两个凸包的点都是逆时针顺序的

所以只需要循环位移$A,B$两个凸包上的点，线性的把$y$坐标最小的点移到最前面就可以了

然后在$A,B$两个凸包有序的情况下归并得到$A+B$

• 假设现在做到$A_i$和$B_j$两个点，把$A_i+B_j$放入集合
• 用叉积比较$A_{i+1}-A_i,B_{j+1}-B_j$​​的极角序，极角序较小的指针$+1$​，如果相同则都$+1$

可能会有三点共线的情况，只需要丢进去来一遍求凸包即可

vec v1[maxn], v2[maxn];
point PoA[maxn], PoB[maxn];

int minkowski( int n, int m, point *Po ) {
	for( int i = 0;i < n;i ++ )
		v1[i] = PoA[( i + 1 ) % n] - PoA[i];
	for( int i = 0;i < m;i ++ )
		v2[i] = PoB[( i + 1 ) % m] - PoB[i];
	Po[0] = PoA[0] + PoB[0];
	int i = 0, j = 0, cnt = 1;
	while( i < n && j < m ) {
		if( dcmp( cross( v1[i], v2[j] ) ) >= 0 )
			Po[cnt] = Po[cnt - 1] + v1[i ++], cnt ++;
		else
			Po[cnt] = Po[cnt - 1] + v2[j ++], cnt ++;
	}
	while( i < n ) Po[cnt] = Po[cnt - 1] + v1[i ++], cnt ++;
	while( j < m ) Po[cnt] = Po[cnt - 1] + v2[j ++], cnt ++;
	return cnt;
} 

闵可夫斯基和也可以用来求解凸包间的最短距离

• 在两个凸包内分别任意取两个点$a,b$
• 这两个点作差得到$\overrightarrow{d}=a-b$
  那么由$a,b$​​的任意性，对于所有$a,b$，对$\overrightarrow{d}$​取最小值就是答案
• 把$B$​​凸包的点全部取反(横纵坐标取相反数)得到$B^{\prime }$​​
• 计算闽可夫斯基和$A+B^{\prime }$​，实际上计算的是$A-B$
• 然后计算点$(0,0)$​到凸包$A-B$​的最短距离即可

求定点到凸包的最短距离

• 枚举凸包的边，求点到线段的距离
• 参考动态凸包的方法，以凸包内一点为原点进行极角排序，二分找到目标线段求距离

欢迎来到——各大模板训练场！！！

Segments——判断线段相交

problem

如果所有线段在某条直线上的投影有交点，那么从这个交点引出一条垂直于直线的线必定交于所有线段

把垂直的线稍微掰一掰，直到恰好穿过两个端点

所以枚举任意两个端点引出一条直线，判断这直线是否与所有线段相交即可

两个端点可能原来就隶属于统一线段，但一定不能是同一个点

#include <cstdio>
#include <cmath>
using namespace std;
#define eps 1e-8
#define maxn 105

struct point {
	double x, y;
	point() {}
	point( double X, double Y ) { x = X, y = Y; }
	point operator + ( point t ) { return point( x + t.x, y + t.y ); }
	point operator - ( point t ) { return point( x - t.x, y - t.y ); }
	friend double cross( point s, point t ) { return s.x * t.y - s.y * t.x; }
	double dot( point t ) { return x * t.x + y * t.y; }
	double len() { return sqrt( dot( *this ) ); }
};

struct line {
	point s, t;
	line(){}
	line( point S, point T ) { s = S, t = T; }
}Line[maxn];

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

bool intersect( line l, line r ) {
	return dcmp( cross( l.s - r.s, l.t - r.s ) ) * dcmp( cross( l.s - r.t, l.t - r.t ) ) <= 0;
}

int T, n;

bool check( line l ) {
	if( ! dcmp( ( l.t - l.s ).len() ) ) return 0;
	for( int i = 1;i <= n;i ++ )
		if( ! intersect( l, Line[i] ) ) return 0;
	return 1;
}

int main() {
	scanf( "%d", &T );
	again :
	while( T -- ) {
		scanf( "%d", &n );
		for( int i = 1;i <= n;i ++ ) {
			double lx, ly, rx, ry;
			scanf( "%lf %lf %lf %lf", &lx, &ly, &rx, &ry );
			Line[i] = line( point( lx, ly ), point( rx, ry ) );
		}
		for( int i = 1;i <= n;i ++ )
			for( int j = i;j <= n;j ++ )
				if( check( line( Line[i].s, Line[j].s ) ) || check( line( Line[i].s, Line[j].t ) ) || check( line( Line[i].t, Line[j].s ) ) || check( line( Line[i].t, Line[j].t ) ) ) {
					printf( "Yes!\n" );
					goto again;
				}
		printf( "No!\n" );
	}
	return 0;
}

Points Within——点与多边形位置

problem

模板判断点是否在多边形内部

#include <cmath>
#include <ctime>
#include <cstdio>
#include <vector>
#include <iostream>
#include <algorithm>
using namespace std;
#define eps 1e-8

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

struct vec {
	double x, y;
	vec(){}
	vec( double X, double Y ) { x = X, y = Y; }
	friend double cross( vec s, vec t ) { return s.x * t.y - s.y * t.x; }
	friend double dot( vec s, vec t ) { return s.x * t.x + s.y * t.y; }
};

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	vec operator - ( point t ) { return vec( x - t.x, y - t.y ); }
};

struct line {
	point p; vec v;
	line(){}
	line( point P, vec V ) { p = P, v = V; }
	friend bool online( line l, point p ) { return dcmp( cross( l.v, p - l.p ) ) == 0; }
};

vector < point > polygon;

bool onseg( point p1, point p2, point p ) {
	line l( p1, p2 - p1 );
	if( ! online( l, p ) ) return 0;
	else if( dcmp( dot( p - p1, p - p2 ) ) <= 0 ) return 1;
	else return 0;
}

bool inangle( vec v1, vec v2, vec v ) {
	if( dcmp( cross( v1, v2 ) ) < 0 ) swap( v1, v2 );
	return cross( v1, v ) >= 0 && cross( v2, v ) <= 0;
}

bool inpolygon( point p ) {
	int n = polygon.size(), ans = 0;
	double r = ( rand() / double( RAND_MAX ) - 0.5 ) * 2 / M_PI;
	vec v( cos( r ), sin( r ) ); 
	for( int i = 0;i < n;i ++ ) {
		if( onseg( polygon[i], polygon[( i + 1 ) % n], p ) ) return 1;
		if( inangle( polygon[i] - p, polygon[( i + 1 ) % n] - p, v ) ) ans ^= 1;
	}
	return ans;
}
int main() {
	int n, m, T = 1;
	while( scanf( "%d", &n ) && n ) {
		scanf( "%d", &m );
		polygon.clear();
		for( int i = 1;i <= n;i ++ ) {
			double x, y;
			scanf( "%lf %lf", &x, &y );
			polygon.push_back( point( x, y ) );
		}
		if( T != 1 ) printf( "\n" );
		printf( "Problem %d:\n", T ), T ++;
		while( m -- ) {
			double x, y;
			scanf( "%lf %lf", &x, &y );
			point p( x, y );
			if( inpolygon( p ) ) printf( "Within\n" );
			else printf( "Outside\n" );
		}
	}
	return 0;
}

Rotational Painting——找重心、求凸包、判断点在线段上

problem

先找到重心，然后求凸包，枚举凸包上的线段，判断重心垂足是否落在线段上

#include <cmath>
#include <cstdio>
#include <vector>
#include <algorithm>
using namespace std;
#define maxn 50005
#define eps 1e-10

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

struct vec {
	double x, y;
	vec(){}
	vec( double X, double Y ) { x = X, y = Y; }
	vec operator + ( vec v ) { return vec( x + v.x, y + v.y ); }
	vec operator - ( vec v ) { return vec( x - v.x, y - v.y ); }
	vec operator * ( double k ) { return vec( x * k, y * k ); }
	double dot() { return x * x + y * y; }
	friend double dot( vec u, vec v ) { return u.x * v.x + u.y * v.y; }
	friend double cross( vec u, vec v ) { return u.x * v.y - u.y * v.x; }
};

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	vec operator - ( point p ) { return vec( x - p.x, y - p.y ); };
	point operator + ( vec v ) { return point( x + v.x, y + v.y ); }
	friend double cross( point s, point t ) { return s.x * t.y - s.y * t.x; }
};

struct line {
	point p; vec v;
	line(){}
	line( point P, vec V ) { p = P, v = V; }
};

vector < point > polygon;
int T, n;

point center() {
	double s = 0; point p( 0, 0 );
	for( int i = 0;i < n;i ++ ) {
		double t = cross( polygon[i], polygon[( i + 1 ) % n] );
		s += t;
		p.x += ( polygon[i].x + polygon[( i + 1 ) % n].x ) * t;
		p.y += ( polygon[i].y + polygon[( i + 1 ) % n].y ) * t;
	}
	p.x /= 3, p.x /= s, p.y /= 3, p.y /= s;
	return p;
}

bool check( point p1, point p2, point p ) {
	return dcmp( dot( p - p1, p2 - p1 ) ) > 0 && dcmp( dot( p - p2, p1 - p2 ) ) > 0;
}

point Po[maxn], sta[maxn], p[maxn];

int graham() {
	for( int i = 0;i < n;i ++ ) p[i] = Po[i];
	int o = 0;
	for( int i = 1;i < n;i ++ )
		if( p[i].y < p[o].y || ( p[i].y == p[o].y && p[i].x < p[o].x ) ) o = i;
	swap( p[o], p[0] );
	sort( p + 1, p + n, []( point s, point t ) {
		vec v1 = s - p[0], v2 = t - p[0];
		double c = cross( v1, v2 );
		return dcmp( c ) > 0 || ( dcmp( c ) == 0 && dcmp( v1.dot() - v2.dot() ) < 0 );
	} );
	sta[0] = p[0], sta[1] = p[1]; int top = 1;
	for( int i = 2;i < n;i ++ ) {
		while( top && dcmp( cross( sta[top] - sta[top - 1], p[i] - sta[top - 1] ) ) <= 0 ) top --;
		sta[++ top] = p[i];
	}
	return top + 1;
}

int main() {
	scanf( "%d", &T );
	while( T -- ) {
		scanf( "%d", &n );
		polygon.clear();
		for( int i = 0;i < n;i ++ ) {
			double x, y;
			scanf( "%lf %lf", &x, &y );
			Po[i] = point( x, y );
			polygon.push_back( Po[i] );
		}
		point c = center();
		int m = graham(), ans = 0;
		for( int i = 0;i < m;i ++ )
			if( check( sta[i], sta[( i + 1 ) % m], c ) ) ans ++;
		printf( "%d\n", ans );
	}
	return 0;
}

Professor’s task——数据结构(set)维护动态凸包

problem

#include <set>
#include <cmath>
#include <cstdio>
using namespace std;
#define int long long

struct point {
	int x, y, dis;
	double rad;
	point(){}
	point( int X, int Y ) { x = X, y = Y; }
	point operator + ( point p ) { return point( x + p.x, y + p.y ); }
	point operator - ( point p ) { return point( x - p.x, y - p.y ); }
	friend int cross( point s, point t ) { return s.x * t.y - s.y * t.x; }
}Po[5], o;

bool operator < ( const point s, const point t ) { 
	return s.rad < t.rad || ( s.rad == t.rad && s.dis < t.dis );
}

set < point > s;
#define it set < point > :: iterator

it left( it i ) {
	if( i == s.begin() ) i = s.end();
	i --;
	return i;
}

it right( it i ) {
	i ++;
	if( i == s.end() ) i = s.begin();
	return i;
}

int dis( point p ) {
	return ( p.x - o.x ) * ( p.x - o.x ) + ( p.y - o.y ) * ( p.y - o.y );
}

int Q;

signed main() {
	scanf( "%lld", &Q );
	for( int i = 1, x;i <= 3;i ++ ) {
		scanf( "%lld %lld %lld", &x, &Po[i].x, &Po[i].y );
		o = o + Po[i], Po[i].x *= 3, Po[i].y *= 3;
	}
	for( int i = 1;i <= 3;i ++ ) {
		Po[i].dis = dis( Po[i] );
		Po[i].rad = atan2( Po[i].y - o.y, Po[i].x - o.x );
		s.insert( Po[i] );
	}
	Q -= 3;
	int opt, x, y;
	while( Q -- ) {
		scanf ( "%lld %lld %lld", &opt, &x, &y );
		x *= 3, y *= 3;
		point p( x, y );
		p.dis = dis( p );
		p.rad = atan2( y - o.y, x - o.x );
		it r = s.lower_bound( p );
		if( r == s.end() ) r = s.begin();
		it l = left( r );
		if( opt & 1 ) {
			if( cross( p - point(*l), point(*r) - point(*l) ) <= 0 ) continue;
			s.insert( p );
			it now = s.find( p );
			l = left( now ); it lst = left( l );
			while( cross( point(*l) - point(*lst), point(*now) - point(*lst) ) <= 0 )
				s.erase( l ), l = lst, lst = left( lst );
			r = right( now ); it nxt = right( r );
			while( cross( point(*r) - point(*now), point(*nxt) - point(*now) ) <= 0 )
				s.erase( r ), r = nxt, nxt = right( nxt );
		}
		else {
			if( cross( p - point(*l), point(*r) - point(*l) ) <= 0 ) printf( "YES\n" );
			else printf( "NO\n" );
		}
	}
	return 0;
}

The Great Divide——判断凸包是否相交

problem

 #include <bits/stdc++.h>
using namespace std;
#define eps 1e-10
#define maxn 505

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

struct vec {
	double x, y;
	vec(){}
	vec( double X, double Y ) { x = X, y = Y; }
	vec operator + ( vec v ) { return vec( x + v.x, y + v.y ); }
	vec operator - ( vec v ) { return vec( x - v.x, y - v.y ); }
	vec operator * ( double k ) { return vec( x * k, y * k ); }
	double dot() { return x * x + y * y; }
	friend double dot( vec u, vec v ) { return u.x * v.x + u.y * v.y; }
	friend double cross( vec u, vec v ) { return u.x * v.y - u.y * v.x; }
};

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	point operator + ( vec v ) { return point( x + v.x, y + v.y ); }
	vec operator - ( point p ) { return vec( x - p.x, y - p.y ); }
}R[maxn], B[maxn], hallR[maxn], hallB[maxn], p[maxn];

struct line {
	point p; vec v;
	line(){}
	line( point P, vec V ) { p = P, v = V; }
	friend bool online( line l, point p ) { return dcmp( cross( l.v, p - l.p ) ) == 0; }
};

bool onseg( point p1, point p2, point p ) {
	line l( p1, p2 - p1 );
	if( ! online( l, p ) ) return 0;
	if( dcmp( dot( p - p1, p - p2 ) ) <= 0 ) return 1;
	else return 0;
}

bool seg_intersect( point p1, point p2, point p3, point p4 ) {
	double s1 = cross( p2 - p1, p3 - p1 ), s2 = cross( p2 - p1, p4 - p1 );
	double s3 = cross( p4 - p3, p1 - p3 ), s4 = cross( p4 - p3, p2 - p3 );
	int ret1 = dcmp( s1 ) * dcmp( s2 ), ret2 = dcmp( s3 ) * dcmp( s4 );
	if( ret1 > 0 || ret2 > 0 ) return 0;
	else if( ret1 == 0 && ret2 == 0 ) {
		if( onseg( p1, p2, p3 ) || onseg( p1, p2, p4 ) ) return 1;
		else return 0;
	} else return 1;
}

int graham( int n, point *Po, point *sta ) {
	for( int i = 0;i < n;i ++ ) p[i] = Po[i];
	sort( p, p + n, []( point s, point t ) { return s.x == t.x ? s.y < t.y : s.x < t.x; } );
	int top = 0;
	for( int i = 0;i < n;i ++ ) {
		while( top > 1 && cross( sta[top - 1] - sta[top - 2], p[i] - sta[top - 1] ) <= 0 )
			top --;
		sta[top ++] = p[i];
	}
	int t = top;
	for( int i = n - 2;~ i;i -- ) {
		while( top > t && cross( sta[top - 1] - sta[top - 2], p[i] - sta[top - 1] ) <= 0 )
			top --;
		sta[top ++] = p[i];
	}
	if( n > 1 ) top --;
	return top;
}

bool inpolygon_rad( point p, point *polygon, int n ) {
	int cnt = 0;
	for( int i = 0;i < n;i ++ ) {
		if( onseg( polygon[i], polygon[( i + 1 ) % n], p ) ) return 1;
		double c = cross( polygon[i] - p, polygon[( i + 1 ) % n] - p );
		double d1 = polygon[i].y - p.y, d2 = polygon[( i + 1 ) % n].y - p.y;
		if( c < 0 && d1 <= 0 && d2 > 0 ) cnt ++;
		if( c > 0 && d2 <= 0 && d1 > 0 ) cnt --;
	}
	return cnt != 0;
}

int main() {
	int n, m;
	next : 
	while( scanf( "%d %d", &n, &m ) ) {
		if( ! n && ! m ) return 0;
		for( int i = 0;i < n;i ++ )
			scanf( "%lf %lf", &R[i].x, &R[i].y );
		for( int i = 0;i < m;i ++ )
			scanf( "%lf %lf", &B[i].x, &B[i].y );
		n = graham( n, R, hallR );
		m = graham( m, B, hallB );
		for( int i = 0;i < n;i ++ )
			if( inpolygon_rad( R[i], hallB, m ) ) { printf( "No\n" ); goto next; }
		for( int i = 0;i < m;i ++ )
			if( inpolygon_rad( B[i], hallR, n ) ) { printf( "No\n" ); goto next; }
		for( int i = 0;i < n;i ++ )
			for( int j = 0;j < m;j ++ )
				if( seg_intersect( hallR[i], hallR[( i + 1 ) % n], hallB[j], hallB[( j + 1 ) % m] ) ) {
					printf( "No\n" );
					goto next;
				}
		printf( "Yes\n" );	
	}
	return 0;
}

Area——求凸包，求多变形面积，pick定理

problem

从$(0,0)$开始，题目给的相当于是$n$​个向量帮忙搭建凸包

套用pick定理

#include <cmath>
#include <cstdio>
#include <vector>
using namespace std;

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	point operator + ( point p ) { return point( x + p.x, y + p.y ); }
	friend double cross( point s, point t ) { return s.x * t.y - s.y * t.x; }
};

vector < point > polygon;
int T, n;

double area() {
	double s = 0;
	for( int i = 0;i < n;i ++ )
		s += cross( polygon[i], polygon[( i + 1 ) % n] );
	return fabs( s );
}

int gcd( int x, int y ) {
	if( ! y ) return x;
	else return gcd( y, x % y );
}

int main() {
	scanf( "%d", &T );
	for( int Case = 1;Case <= T;Case ++ ) {
		polygon.clear();
		scanf( "%d", &n );
		point p( 0, 0 );
		polygon.push_back( p );
		for( int i = 1;i <= n;i ++ ) {
			double x, y;
			scanf( "%lf %lf", &x, &y );
			p = p + point( x, y );
			polygon.push_back( p );
		}
		if( Case > 1 ) printf( "\n" );
		printf( "Scenario #%d:\n", Case );
		double s = area(); int v = 0;
		for( int i = 0;i < n;i ++ ) {
			point p1 = polygon[i], p2 = polygon[( i + 1 ) % n];
			v += gcd( fabs( p1.x - p2.x ), fabs( p1.y - p2.y ) );
		}
		printf( "%d %d %.1f\n", ( int( s ) + 2 - v ) / 2 , v, s / 2 );
	}
	return 0;
}

That Nice Euler Circuit——欧拉定理

套用欧拉定理，注意重合点要删掉

#include <cmath>
#include <cstdio>
#include <algorithm>
using namespace std;
#define eps 1e-9
#define maxn 305

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

struct vec {
	double x, y;
	vec(){}
	vec( double X, double Y ) { x = X, y = Y; }
	vec operator * ( double k ) { return vec( x * k, y * k ); }
	friend double dot( vec u, vec v ) { return u.x * v.x + u.y * v.y; }
	friend double cross( vec u, vec v ) { return u.x * v.y - u.y * v.x; }
};

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	point operator + ( vec v ) { return point( x + v.x, y + v.y ); }
	vec operator - ( point p ) { return vec( x - p.x, y - p.y ); }
}Po[maxn], V[maxn * maxn];

struct line {
	point p; vec v;
	line(){}
	line( point P, vec V ) { p = P, v = V; }
};

bool strict_intersect( point p1, point p2, point p3, point p4 ) {
	double s1 = cross( p2 - p1, p3 - p1 ), s2 = cross( p2 - p1, p4 - p1 );
	double s3 = cross( p4 - p3, p1 - p3 ), s4 = cross( p4 - p3, p2 - p3 );
	return dcmp( s1 ) * dcmp( s2 ) < 0 && dcmp( s3 ) * dcmp( s4 ) < 0;
}

point intersection( line l, line r ) {
	double s = cross( l.v, r.v );
	if( dcmp( s ) != 0 ) {
		double t = cross( r.v, l.p - r.p ) / s;
		return l.p + l.v * t;
	}
}

bool operator < ( const point &p1, const point &p2 ) {
	return p1.x == p2.x ? p1.y < p2.y : p1.x < p2.x;
}

bool operator == ( const point &p1, const point &p2 ) {
	return dcmp( p1.x - p2.x ) == 0 && dcmp( p1.y - p2.y ) == 0;
}

bool onseg( point l, point r, point p ) {
	return dcmp( cross( l - p, r - p ) ) == 0 && dcmp( dot( l - p, r - p ) ) < 0;
}

int main() {
	int n, T = 0;
	while( scanf( "%d", &n ) && n ) {
		for( int i = 0;i < n;i ++ )
			scanf( "%lf %lf", &Po[i].x, &Po[i].y ), V[i] = Po[i];
		n --; int cnt = n, E = n;
		for( int i = 0;i < n;i ++ )
			for( int j = i + 1;j < n;j ++ )
				if( strict_intersect( Po[i], Po[i + 1], Po[j], Po[j + 1] ) )
					V[cnt ++] = intersection( line( Po[i], Po[i + 1] - Po[i] ), line( Po[j], Po[j + 1] - Po[j] ) );
		sort( V, V + cnt );
		cnt = unique( V, V + cnt ) - V;
		for( int i = 0;i < cnt;i ++ )
			for( int j = 0;j < n;j ++ )
				if( onseg( Po[j], Po[j + 1], V[i] ) ) E ++;
		printf( "Case %d: There are %d pieces.\n", ++ T, 2 - cnt + E );
	}
	return 0;
}

Buried memory——求三角形的外心、最小圆覆盖

problem

#include <bits/stdc++.h>
using namespace std;
#define maxn 100005
#define eps 1e-12

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	friend double len( point s, point t ) { return sqrt( ( s.x - t.x ) * ( s.x - t.x ) + ( s.y - t.y ) * ( s.y - t.y ) ); }
	point operator + ( point p ) { return point( x + p.x, y + p.y ); }
	point operator / ( double k ) { return point( x / k, y / k ); }
	double operator - ( point p ) { return x * x + y * y - p.x * p.x - p.y * p.y; }
}p[maxn];

double X( point p1, point p2 ) {
	return p1.x - p2.x;
}

double Y( point p1, point p2 ) {
	return p1.y - p2.y;
}

point get_c( point p1, point p2, point p3 ) {
	double t = 2 * ( Y( p2, p1 ) * X( p3, p1 ) - Y( p3, p1 ) * X( p2, p1 ) );
	return point( ( Y( p2, p1 ) * ( p3 - p1 ) - Y( p3, p1 ) * ( p2 - p1 ) ) / t, ( X( p3, p1 ) * ( p2 - p1 ) - X( p2, p1 ) * ( p3 - p1 ) ) / t );
}

int n; point c; double r;
void circle_cover() {
	random_shuffle( p + 1, p + n + 1 );
	c = p[1], r = 0;
	for( int i = 2;i <= n;i ++ )
		if( dcmp( len( c, p[i] ) - r ) > 0 ) {
			c = p[i], r = 0;
			for( int j = 1;j < i;j ++ )
				if( dcmp( len( c, p[j] ) - r ) > 0 ) {
					c = ( p[i] + p[j] ) / 2, r = len( c, p[i] );
					for( int k = 1;k < j;k ++ )
						if( dcmp( len( c, p[k] ) - r ) > 0 )
							c = get_c( p[i], p[j], p[k] ), r = len( c, p[i] );
				}
		}
}

int main() {
	while( scanf( "%d", &n ) && n ) {
		for( int i = 1;i <= n;i ++ )
			scanf( "%lf %lf", &p[i].x, &p[i].y );
		circle_cover();
		printf( "%.2f %.2f %.2f\n", c.x, c.y, r );
	}
	return 0;
}

「ICPC World Finals 2017」机场构建 Airport Construction——求凹凸多边形内的最长线段

source

#include <cmath>
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
#define maxn 205
#define double long double
const double eps = 1e-7;
int n;
double ans;

int dcmp( double x ) {
	return fabs( x ) < eps ? 0 : ( x > 0 ? 1 : -1 );
}

struct vec {
	double x, y;
	vec(){}
	vec( double X, double Y ) { x = X, y = Y; }
	vec operator + ( vec v ) { return vec( x + v.x, y + v.y ); }
	vec operator - ( vec v ) { return vec( x - v.x, y - v.y ); }
	vec operator * ( double k ) { return vec( x * k, y * k ); }
	friend double len( vec v ) { return v.x * v.x + v.y * v.y; }
	friend double dis( vec v ) { return sqrt( len( v ) ); }
	friend double dot( vec u, vec v ) { return u.x * v.x + u.y * v.y; }
	friend double cross( vec u, vec v ) { return u.x * v.y - u.y * v.x; }
	friend bool order( vec v ) { return dcmp( v.y ) < 0 || ( dcmp( v.y ) == 0 && dcmp( v.x ) < 0 ); }//order规定一种指向 
	friend bool sameline( vec u, vec v ) { return dcmp( cross( u, v ) ) == 0; }//共线判断:两向量同向或反向 叉积面积为0 一维直线在二维的面积为0 
	friend bool samedirection( vec u, vec v ) { return sameline( u, v ) && dcmp( dot( u, v ) ) >= 0; }//同向判断:首先共线然后点积为非负->锐角 
};

struct point {
	double x, y;
	point(){}
	point( double X, double Y ) { x = X, y = Y; }
	point operator + ( point p ) { return point( x + p.x, y + p.y ); }
	point operator + ( vec v ) { return point( x + v.x, y + v.y ); }
	point operator / ( double k ) { return point( x / k, y / k ); }
	vec operator - ( point p ) { return vec( x - p.x, y - p.y ); }
}polygon[maxn]; 

struct line {
	point p; vec v;
	line(){}
	line( point P, vec V ) { p = P, v = V; }
}seg[maxn];

bool onseg( point l, point r, point p ) {
	return dcmp( cross( l - p, r - p ) ) == 0 && dcmp( dot( l - p, r - p ) ) <= 0;
}

bool onseg( line MS, point p ) {
	point l = MS.p, r = MS.p + MS.v;
	return onseg( l, r, p );
}

bool inpolygon( point p ) {
	int cnt = 0;
	for( int i = 0;i < n;i ++ ) {
		if( onseg( polygon[i], polygon[( i + 1 ) % n], p ) ) return 1;
		double c = cross( polygon[i] - p, polygon[( i + 1 ) % n] - p );
		double d1 = polygon[i].y - p.y, d2 = polygon[( i + 1 ) % n].y - p.y;
		if( c < 0 && d1 <= 0 && d2 > 0 ) cnt ++;
		if( c > 0 && d2 <= 0 && d1 > 0 ) cnt --;
	}
	return cnt != 0;
}

point seg_intersect( line l, point p1, point p2 ) {
	line r( p1, p2 - p1 );
	return l.p + l.v * ( cross( r.v, l.p - r.p ) / cross( l.v, r.v ) );
}

point p[maxn];

void solve( line MS ) {
	int cnt = 0;
	for( int i = 0;i < n;i ++ )
		if( sameline( MS.v, polygon[( i + 1 ) % n] - polygon[i] ) ) continue;//与int main里跳过线段同理 
		else {
			point intersection = seg_intersect( MS, polygon[i], polygon[( i + 1 ) % n] );
			if( onseg( polygon[i], polygon[( i + 1 ) % n], intersection ) ) p[++ cnt] = intersection;
			//如果线段与多边形的某些线段有交 记录下来
			//多边形可能是凹多边形 所以这种情况是可能出现的 
		}
	sort( p + 1, p + cnt + 1, []( point s, point t ) { return s.x == t.x ? s.y < t.y : s.x < t.x; } );
	int l = 1;
	while( l < cnt && ! onseg( MS, p[l]) ) l ++;
	int r = l;
	while( r < cnt && onseg( MS, p[r + 1] ) ) r ++;
	//两个while确定MS线段的首尾端点都被包含在多边形内 
	for( int i = l;i < r;i ++ )
		if( ! inpolygon( ( p[i] + p[i + 1] ) / 2 ) ) return;//然后新的整MS线段要被完全包含在多边形内才行 看图
	while( l > 1 && inpolygon( ( p[l] + p[l - 1] ) / 2 ) ) l --;//然后尽可能的往外拓展 p点是拍了序的 所以可以保证 
	while( r < cnt && inpolygon( ( p[r] + p[r + 1] ) / 2 ) ) r ++;
	ans = max( ans, dis( p[r] - p[l] ) );
}

int main() {
	scanf( "%d", &n );
	for( int i = 0, x, y;i < n;i ++ ) {
		scanf( "%d %d", &x, &y );
		polygon[i] = point( x, y ); 
	}
	for( int i = 0;i < n;i ++ ) {
		int ip = 0;
		for( int j = 0;j < n;j ++ )
			if( i == j || order( polygon[j] - polygon[i] ) ) continue;
			else seg[++ ip] = line( polygon[i], polygon[j] - polygon[i] );
		sort( seg + 1, seg + ip + 1, []( line l, line r ) {//极角排序建构多边形边集 
			int c = dcmp( cross( l.v, r.v ) );
			return c > 0 || ( c == 0 && dcmp( len( l.v ) - len( r.v ) ) < 0 );
		} );
		for( int j = 1;j <= ip;j ++ )
			if( j > 1 && sameline( seg[j - 1].v, seg[j].v ) ) continue;
			/*
			solve做的操作是把一条线段尽可能地边长撑到多边形边界为止 所以一个直线上的线段做一次就可以了 
			用共线sameline或者同向samedirection判断都是一样的
			因为已经极角排序过了->逆时针 
			*/ 
			else solve( seg[j] ); 
	}
	printf( "%.9Lf\n", ans );
	return 0;
}