Logo root的博客

博客

计算几何初步

2023-11-23 17:50:32 By root

出自:https://www.luogu.com.cn/blog/wjyyy/geometry1

引言

计算几何对于很多刚接触OI的同学们是不会涉及的,当逐步进入省选范畴之后,计算几何的应用就非常广了。

你需要知道?

在学习计算几何(至少是这篇文章)之前,你需要了解以下这些东西:
平面直角坐标系
向量及其运算
三角剖分求面积

平面直角坐标系

平面直角坐标系就是笛卡尔坐标系,这一部分会在人教版七年级数学中系统地学习。
在这篇文章中,我们会经常用到坐标表示,包括下面的向量。同时,我们的点和向量都是通过坐标来保存的。

向量及其运算

向量

向量是既有大小又有方向的量,可以用一个带有箭头的线段表示。
如果一个向量的起点是点A,终点是点B,那么这个向量可以被表示为$\vec{AB}$ ,或者用一个小写英文字母表示为$\vec{a}$⃗或$\mathbf{a}$(教科书上加粗时也是斜体)。
向量具有平移不变性,对于任意一个向量,都可以把它的起点平移到原点。因此我们可以在坐标系中用一对有序实数对来表示向量,也就是这个向量起点为原点时终点的坐标(x,y)。

向量的模

模就是模长,向量的模就是向量的长度。表示为$|\vec{a}|$,对于 $\vec{a}$ =(x,y),$|\vec{a}|=\sqrt{x^2+y^2}$

相反向量

如果两个向量大小相等,方向相反,那么这两个向量互为相反向量,$\vec{a}$的相反向量表示为$-\vec{a}$,$|\vec{a}|=|-\vec{a}|$ .对于 $\vec{a}$ =(x,y),$|-\vec{a}|=(-x,-y)$

垂直向量

如果两个向量所在的直线互相垂直,那么这两个向量互称垂直向量,表示为$\vec{a} \perp \vec{b}$ 。

共线向量

如果两个向量所在的直线互相平行,那么这两个向量互为共线向量,表示为$\vec{a} \parallel \vec{b}$ 。因为两向量平行,所以平移其中一个向量,总能和另一个向量共线。

向量的加减

前面提到了向量具有平移不变性,那么$\vec{a} + \vec{b}$就相当于把向量$\vec{a}, \vec{b}$ 首尾相接(假设$\vec{a}$的首接在$\vec{b}$的尾),此时$\vec{a}$的尾指向$\vec{b}$ 的首就是相加后的向量。对于$\vec{a}=(x_1,y_1)$,$\vec{b}=(x_2,y_2)$,$\vec{a}+\vec{b} =(x_1+x_2,y_1+y_2)$.
向量的减法实际上是一个向量加上另一个向量的相反向量,也就是$\vec{a}-\vec{b}$=$\vec{a} + (-\vec{b})$,可以先把$\vec{b}$向量变成$-\vec{b}$ ,然后再进行加法计算。 $\vec{a}=(x_1,y_1)$,$\vec{b}=(x_2,y_2)$,$\vec{a}-\vec{b} =(x_1-x_2,y_1-y_2)$.
向量的加减法具有交换律。

零向量

零向量的模长为0,方向任意,它与任何一个向量都共线,但在高中数学中不认为它与任何一个向量都垂直。坐标为(0,0)。

向量的数乘

向量的数乘实际上是对向量进行放缩。数乘的形式 是⃗ $ \lambda\vec{a}$,其中常数$\lambda \in R$是对向量放缩的程度,可正可负可为零。数乘不改变向量的方向,因此 $\vec{a}\parallel \lambda\vec{a}$ 。对于$\vec{a}=(x,y)$,$\lambda\vec{a}=(\lambda x$,$\lambda y)$。

向量的数量积

向量的数量积,也叫点积、内积,几何意义为一个向量在另一个向量上的投影再乘上第二个向量的模长。$\vec{a}$ 与$\vec{b}$的数量积表示为$\vec{a}\cdot\vec{b}$ ,是一个实数。计算式为 $\vec{a}\cdot\vec{b}=|\vec{a}||\vec{b}|\cos\theta$ ($\theta$为夹角)。引入坐标后,通过三角恒等变换的推导,对于$\vec{a}=(x_1,y_1), \vec{b}=(x_2,y_2),\vec{a}\cdot\vec{b}=x_1x_2+y_1y_2$。

如果两个向量同向(共线),那么它们的数量积为他们的模长之积。
如果两个向量夹角$\lt 90^\circ$,那么它们的数量积为正。
如果两个向量夹角$=90^\circ$,那么他们的数量积为0,因为$\cos 90^\circ=0$。
如果两个向量夹角$\gt 90^\circ$,那么它们的数量积为负。
如果两个向量反向(共线),那么它们的数量积为他们的模长之积的相反数。
由向量的数量积可以判断它们的夹角和投影。向量的数量积具有交换律。

向量的外积

向量的外积,也叫叉积,几何意义是两向量由平行四边形法则围成的面积。外积是一个向量,垂直于原来两个向量所在的平面。对于$\vec{a}=(x_1,y_1),\vec{b}=(x_2,y_2),\vec{a}\times\vec{b}$的方向就要用右手从$\vec{a}$ 沿着不大于平角的方向旋转,大拇指的方向就是外积的方向。如果外积方向朝外,那么值为正,否则值为负。计算式为$\vec{a}\times\vec{b}=|\vec{a}||\vec{b}|\sin\theta=x_1y_2-x_2y_1$,因为有谁转向谁的区别,所以外积没有交换律。 $\vec{a}\times\vec{b}$的正负也可以理解为,把$\vec{a}$逆时针方向转到$\vec{b}$的方向,夹角为$\theta$。当$0\leq\theta\lt\pi$时值为正;当$\pi\leq\theta\lt 2\pi$时值为负。

如果$\vec{a}\parallel\vec{b}$ ,那么由$\sin 0^\circ=\sin 180^\circ=0,\vec{a}\times\vec{b}=\vec{0}$,得到的数字结果也是0。
如果$\vec{b}$的终点在$\vec{a}$的左侧(假设已经共起点),那么$\vec{a}\times\vec{b}\gt 0,\vec{a}$握向$\vec{b}$转过的角逆时针不超过$180^\circ$,结果为正。
如果$\vec{b}$的终点在$\vec{a}$的右侧,那么$\vec{a}\times\vec{b}\lt 0,\vec{a}$握向$\vec{b}$转过的角逆时针超过了$180^\circ$ ,或者说顺时针不超过$180^\circ$,大拇指朝内,结果为负。
由向量外积可以判断两向量的旋转关系,同时可以方便地求出点到直线的距离。分别会用于下面的凸包和旋转卡壳。

一般情况下我们储存点就会把它放在向量中,用向量$\vec{OA}$来表示A 这个点。因此下文的A×B 就表示向量$\vec{OA}$与向量$\vec{OB}$的叉积。

线段交点

upd on 2019.3.12

根据面积我们可以求出两线段交点。

我们对于两个线段AB,CD,可以通过面积和相似来求出交点O的坐标。因为A的坐标已知,实际上是通过求出$\vec{AO}$来确定O的坐标。
我们根据上面的叉积公式求面积,得到$S_{\triangle ACD}=\frac{\vec{AD}\times\vec{AC}}{2}$,$S_{\triangle BCD}=\frac{\vec{BC}\times\vec{AC}}{2}$。 再作垂线$AE\bot CD,BF \bot CD,S_{\triangle ACD}=\frac{|CD|\times|AE|}{2},S_{\triangle BCD}=\frac{|CD|\times|BF|}{2}$。则$\frac{∣AE∣}{BF}=\frac{S_{\triangle ACD}}{S_{\triangle BCD}} $此时△AEO∽△BFO,得到$\frac{|AE|}{|BF|}=\frac{|AO|}{|BO|}$。

又因为$\vec{AO}=\frac{|AO|}{|AB|} \cdot \vec{AB}=\frac{|AO|}{|AO|+|BO|}\cdot\vec{AB}$,令原点为O',则$\vec{O'O}=\vec{O'A}+\vec{AO},\vec{O'O}$的坐标即为O 的坐标。

向量旋转

upd on 2019.3.1

这个知识点虽然在这篇文章里没有应用,但是是计算几何的基本操作技巧。
设$\vec{a}=(x,y)$,倾角为$\theta$,长度为l=$\sqrt{x^2+y^2}$则$x=l\sin\theta,y=l\sin\theta$。令其顺时针旋转α度角,得到向量$b=(l\cos(\theta+\alpha),l\sin(\theta+\alpha))$。
由三角恒等变换得,$\vec{b}=(l(\cos\theta \cos\alpha-\sin\theta\sin\alpha),l(\sin\theta\cos\alpha+\cos\theta\sin\alpha))$
化简,去括号
$\vec{b}$ =(lcosθcosα−lsinθsinα,lsinθcosα+lcosθsinα)
把上面的x,y 代回来得
$\vec{b}$ =(xcosα−ysinα,ycosα+xsinα)
即使不知道三角恒等变换,这个式子也很容易记下来。

三角剖分求面积

三角剖分可以用来求任意多边形的面积。

我们只需要知道一个多边形的各个连续顶点分别在哪里即可。

上面提到了向量的叉积是表示的是有向面积,方向用正负表示,此时正负就派上用场了。

我们如果把相邻每两个顶点与原点构成的向量的叉积的数值的一半依次累加起来,就能得到多边形的面积。如下图:红色面积表示负的,蓝色面积表示正的。 $S_{ABCDEF}=\frac{\vec{OA}\times\vec{OB}+\vec{OB}\times\vec{OC}+\cdots+\vec{OF}\times\vec{OA}}{2}$
首先计算$\vec{OA}\times\vec{OB}$ ,得到负的面积$-S_{\triangle OAB}$ ,在图上表示为红色部分;然后计算$S_{\triangle OBC}$ ,此时在图上表示为抵消$-S_{\triangle OBP}$ ,并额外加上正的$S_{\triangle PBC}$ 。

5.png
三角剖分1

然后依次进行下去
6.png
三角剖分2

可以看出,多边形外部出现过的红色总会被消掉,因为多边形是闭合的。当出现了一条相对原点向右的边时,则之后必然会出现一条相对原点向左的边来封闭这个多边形。

如果原点在多边形内部就更好理解了。多边形凹的部分会被”拐回来的部分“给减掉,剩下的部分都可以直接相加。

最终对于n 边形 $A_1A_2 …A_n$ ,面积计算公式就是$S=x_ny_1-y_nx_1+\sum_{i=1}^{n-1}(x_iy_{i+1}-y_ix_{i+1})$

当用向量存点时,计算可以用相邻向量叉积之和,$S=A_n\times A_1+\sum_{i=1}^{n-1}A_i\times A_{i+1}$。

凸包

定义

凸包就是求一个周长最小的,并且能包围所有给定点的多边形。当多边形表面存在凹陷时,根据三角不等式 $$\begin{cases} a+b>c \\[2ex] b+c>a \\[2ex] a+c>b\end{cases}$$ ,一定没有直接把最短那条边连起来优。

性质

如果按逆时针方向看,凸包上每两条相邻的边都是向左拐的。比如说,与边$\vec{a_i}$顺时针方向相邻的是边$\vec{a_{i+1}$ ,那么对于任意$i\in [1,n)$(n 为凸包上点的个数)有$\vec{a_i}\times\vec{a_{i+1}}>0,\vec{a_n}\times\vec{a_1}>0$。因为用右手从$\vec{a_i}$的方向顺小于平角的一边握向$\vec{a_{i+1}}$ ,大拇指总会指向外边。

求法 首先把所有点以横坐标为第一关键字,纵坐标为第二关键字排序。

显然排序后最小的元素和最大的元素一定在凸包上。而且因为是凸多边形,我们如果从一个点出发逆时针走,轨迹总是 “左拐” 的,一旦出现右拐,就说明这一段不在凸包上。因此我们可以用一个单调栈来维护上下凸壳。

因为从左向右看,上下凸壳所旋转的方向不同,为了让单调栈起作用,我们首先升序枚举求出下凸壳,然后降序求出上凸壳。

求凸壳时,一旦发现即将进栈的点(P)和栈顶的两个点($S_1,S_2$ ,其中$S_1$为栈顶 )行进的方向向右旋转,即叉积小于$0:\vec{S_2S_1}\times\vec{S_1P}\lt 0$,则弹出栈顶,回到上一步,继续检测,直到$\vec{S_2S_1}\times\vec{S_1P}\geq 0$ 或者栈内仅剩一个元素为止。

通常情况下不需要保留位于凸包边上的点,因此上面一段中$\vec{S_2S_1}\times\vec{S_1P}\lt 0$这个条件中的 “<” 可以视情况改为≤,同时后面一个条件应改为>。

最后不要忘了把最小的元素与栈顶进行比较,以保证最后一段也是凸壳。

时间复杂度 O(nlogn)。

代码

std::sort(p+1,p+1+n);
stk[++tp]=1;
for(int i=2;i<=n;++i)
{
    while(tp>1&&(p[stk[tp]]-p[stk[tp-1]])*(p[i]-p[stk[tp]])<=0)
        used[stk[tp--]]=0;
    used[i]=1;
    stk[++tp]=i;
}
int qaq=tp;
for(int i=n-1;i>0;--i)
    if(!used[i])
    {
        while(tp>qaq&&(p[stk[tp]]-p[stk[tp-1]])*(p[i]-p[stk[tp]])<=0)
            used[stk[tp--]]=0;
        used[i]=1;
        stk[++tp]=i;
    }
for(int i=1;i<=tp;++i)
    h[i]=p[stk[i]];

最后h数组中存的就是凸包了。(下标1和tp存的是同一个点)

练习 UVA 11626 Convex Hull

洛谷 P2742 圈奶牛

洛谷P3829 [SHOI2012] 信用卡凸包

旋转卡壳

前言
Warning:前言部分仅供引入,严禁参考。

刚学旋转卡壳时,是看到了一句话:“同时旋转两条平行线直到其中一条和多边形的一条边重合”。那不就是旋转到夹角较小的那一条边上去吗?

$∠1 \le ∠2$,那么我们把平行线逆时针旋转过∠1 的大小,就可以和多边形的一边重合了。

于是就出现了这样的代码:

/***vct结构体内部***/
friend db operator ^(vct a,vct b)//求夹角
{
    db tmp=a.x*b.x+a.y*b.y;
    tmp/=a.dis();
    tmp/=b.dis();
    return fabs(acos(tmp));
}

/***主函数中***/
while(t2<tp)
{
    db d1=t^(h[t1+1]-h[t1]);//夹角1
    db d2=t_^(h[t2+1]-h[t2]);//夹角2
    if(d1<d2)
    {
        ++t1;
        t.rot(d1);
        t_.rot(d1);
    }
    else
    {
        ++t2;
        t.rot(d2);
        t_.rot(d2);
    }
    db tmp=(h[t1]-h[t2]).dis();
    ans=ans>tmp?ans:tmp;
}

然而这样会产生较大的精度误差,显然不是计算几何想要的。

凸多边形的切线

如果一条直线与凸多边形有交点,并且整个凸多边形都在这条直线的一侧,那么这条直线就是该凸多边形的一条切线。

对踵点

如果过凸多边形上两点作一对平行线,使得整个多边形都在这两条线之间,那么这两个点被称为一对对踵点。

凸多边形的直径

即凸多边形上任意两个点之间距离的最大值。直径一定会在对踵点中产生,如果两个点不是对踵点,那么两个点中一定可以让一个点向另一个点的对踵点方向移动使得距离更大。并且点与点之间的距离可以体现为线与线之间的距离,在非对踵点之间构造平行线,一定没有在对踵点构造平行线优,这一点可以通过平移看出。

旋转卡壳

这是一张比较标致的动图,生动地体现了旋转卡壳的过程,NOIWC2017中的《膜你抄》中也用了这张图。

为了减小误差(上面那个傻乎乎的方法误差就比较大,并方便计算,我们求出每一条边的对踵点。通过上面的图可以看出,边的对踵点同时是这两个点的对踵点,因此并没有遗漏。

同时,通过边来求,可以很好地运用叉积的计算,减小了误差。

对于最下面这条边来说,对于上面四个顶点,他们分别与这条边构成了四个同底不同高的三角形,而三角形面积可以用叉积算出来,三角形的面积此时也就反映了点到直线的距离了。由于面积是单峰的,那么对单独一条边的对踵点,我们可以用三分来求,对于所有边,我们可以用two-pointer来O(n) 求出。因为随着边在逆时针枚举,它的对踵点T也在逆时针移动。

补充:对于三角形面积

因为 $$ \begin{cases} \vec{AB} \times \vec{AC} =|AB||AC| \sin \theta \\[2ex] S_{\triangle ABC}=\frac{1}{2}|AB||AC|\sin\theta,\theta=\angle BAC \\[2ex] S_{\triangle ABC}=\frac{1}{2}|BC|h, \end{cases} $$ 三式联立得ℎ=$\frac{\vec{AB} \times \vec{AC}}{|BC|}$ 。
因此可以用旋转卡壳来求凸多边形直径。

步骤

首先求出凸包,目的是把所有点按逆时针排序。然后枚举逆时针边,定义T 为当前的对踵点,若 T 逆时针变动能增大到当前边的距离则改变T,直到继续改变会导致距离减小为止。每次求出边(设为AB)的对踵点$T_i$后,分别用$∣AT_i∣,∣BT_i∣$ 来更新答案,即凸多边形的直径。 时间复杂度O(n)。

代码

db dis(vct a,vct b,vct x)//求出x到ab的距离
{
    return fabs((a-x)*(b-x)/(a-b).dis());
}
int t=1;
db ans=0.0;
for(int i=1;i<tp;++i)
{
    while(dis(h[i],h[i+1],h[t])< = dis(h[i],h[i+1],h[t+1]))
        t=t%(tp-1)+1;
    ans=std::max(ans,std::max((h[i]-h[t]).dis(),(h[i+1]-h[t]).dis()));//这里的dis是向量的模长
}

练习

洛谷 P1452 Beauty Contest

UVA10173 Smallest Bounding Rectangle $O(n^2)$求最小覆盖矩形

洛谷 P3187 [HNOI2007] 最小矩形覆盖 $O(nlogn)$求最小覆盖矩形并输出矩形顶点(bzoj有spj)

半平面交

定义

半平面

一条直线和直线的一侧。半平面是一个点集,因此是一条直线和直线的一侧构成的点集。当包含直线时,称为闭半平面;当不包含直线时,称为开半平面。

解析式一般为 Ax+By+C≥0。

在计算几何中用向量表示,整个题统一以向量的左侧或右侧为半平面。

半平面

半平面交

半平面交是指多个半平面的交集。因为半平面是点集,所以点集的交集仍然是点集。在平面直角坐标系围成一个区域。

这就很像普通的线性规划问题了,得到的半平面交就是线性规划中的可行域。一般情况下半平面交是有限的,经常考察面积等问题的解决。

它可以理解为向量集中每一个向量的右侧的交,或者是下面方程组的解。 $$ \begin{cases} A_1x+B_1y+C \geq 0 \\[2ex] A_2x+B_2y+C \geq 0 \\[2ex] … \end{cases} $$

多边形的核

如果一个点集中的点与多边形上任意一点的连线与多边形没有其他交点,那么这个点集被称为多边形的核。

把多边形的每条边看成是首尾相连的向量,那么这些向量在多边形内部方向的半平面交就是多边形的核。

解法 - S & I算法 极角排序 以前一直不会极角排序太丢人了甚至不会写Graham求凸包 C语言有一个库函数叫做atan2(double y,double x),可以返回θ∈[−π,π](网上很多地方说左边是开的,不过对排序没有影响),$\theta=arctan \frac{y}{x}$。

直接以向量为自变量,调用这个函数,以返回值为关键字排序,得到新的边(向量)集。

排序时,如果遇到共线向量(且方向相同),则取靠近可行域的一个。比如两个向量的极角相同,而我们要的是向量的左侧半平面,那么我们只需要保留左侧的向量。判断方法是取其中一个向量的起点或终点与另一个比较,检查是在左边还是在右边。

维护单调队列

因为半平面交是一个凸多边形,所以需要维护一个凸壳。因为后来加入的只可能会影响最开始加入的或最后加入的边(此时凸壳连通),只需要删除队首和队尾的元素,所以需要用单调队列。

我们遍历排好序了的向量,并维护另一个交点数组。当单队中元素超过2个时,他们之间就会产生交点。

对于当前向量,如果上一个交点在这条向量表示的半平面交的异侧,那么上一条边就没有意义了。

单调队列

如上图,假设取向量左侧半平面。极角排序后,遍历顺序应该是$\vec{a}\rightarrow \vec{b} \rightarrow \vec{c}$。当$\vec{a}$和$\vec{b}$入队时,在交点数组里会产生一个点D(交点数组保存队列中相同下标的向量与前一向量的交点)。
接下来枚举到$\vec{c}$时,发现D 在$\vec{c}$的右侧。而因为产生D 的向量的极角一定比$\vec{⃗c}$ 要小,所以产生D的向量(指$\vec{b}$)就对半平面交没有影响了。

还有一种可能的情况是快结束的时候,新加入的向量会从队首开始造成影响。

队首影响

仍然假设取向量左侧半平面。加入向量$\vec{f}$之后,第一个交点G就在⃗$\vec{f}$的右侧,我们把上面的判断标准逆过来看,就知道此时应该删除向量$\vec{a}$,也即队首的向量。

最后用队首的向量排除一下队尾多余的向量。因为队首的向量会被后面的约束,而队尾的向量不会。此时它们围成了一个环,因此队首的向量就可以约束队尾的向量。

得到半平面交 如果半平面交是一个凸 n 边形,最后在交点数组里会得到n 个点。我们再把它们首尾相连,就是一个统一方向(顺或逆时针)的 n 多边形。

此时就可以用三角剖分求面积了。(求面积是最基础的考法)

偶尔会出现半平面交不存在或面积为0的情况,注意考虑边界。

时间复杂度O(nlogn)。

注意事项 upd on 2019.3.1

昨天晚上 @WG6101N 同学跟我说两个循环写反了会挂得很惨…稍微研究了一下原因。

当出现一个可以把队列里的点全部弹出去的向量(即所有队列里的点都在该向量的右侧),则我们必须先处理队尾,再处理队首。因此在循环中,我们先枚举--r;的部分,再枚举++l;的部分,才不会错。原因如下。

一般情况下,我们在队列 (队列顺序为{$\vec{u},\vec{v}$}) 后面加一条边(向量$\vec{w}$),会产生一个交点N ,缩小$\vec{v}$ 后面的范围。

但是毕竟每次操作都是一般的,因此可能会有把M 点“挤出去”的情况。

如果此时出现了向量 $\vec{a}$ ,使得M 在$\vec{a}$的右侧,那么M 就要出队了。此时如果从队首枚举++l,显然是扩大了范围。实际上M 点是由$\vec{u}$和$\vec{v}$共同构成的,因此需要考虑影响到现有进程的是⃗$\vec{u}$还是$\vec{v}$。而因为我们在极角排序后,向量是逆时针顺序,所以$\vec{v}$的影响要更大一些。 就如上图,如果M 确认在$\vec{a}$的右侧,那么此时$\vec{v}$ 的影响一定不会对半平面交的答案作出任何贡献。

而我们排除队首的原因是当前向量的限制比队首向量要大,这个条件的前提是队列里有不止两个线段(向量),不然就会出现上面的情况。

所以一定要先排除队尾再排除队首。

代码
比较部分

friend bool operator <(seg x,seg y)
{
    db t1=atan2((x.b-x.a).y,(x.b-x.a).x);
    db t2=atan2((y.b-y.a).y,(y.b-y.a).x);//求极角
    if(fabs(t1-t2)>eps)//如果极角不等
        return t1<t2;
    return (y.a-x.a)*(y.b-x.a)>eps;//判断向量x在y的哪边,令最靠左的排在最左边
}

增量部分

//pnt its(seg a,seg b)表示求线段a,b的交点
//s[]是极角排序后的向量
//q[]是向量队列
//t[i]是s[i-1]与s[i]的交点
//【码风】队列的范围是(l,r]
//求的是向量左侧的半平面
int l=0,r=0;
for(int i=1;i<=n;++i)
    if((s[i]==s[i-1])==false)
    {
        //务必先检查队尾再检查队首
        while(r-l>1&&(s[i].b-t[r])*(s[i].a-t[r])>eps)//如果上一个交点在向量右侧则弹出队尾
            --r;
        while(r-l>1&&(s[i].b-t[l+2])*(s[i].a-t[l+2])>eps)//如果第一个交点在向量右侧则弹出队首
            ++l;
        q[++r]=s[i];
        if(r-l>1)
            t[r]=its(q[r],q[r-1]);//求新交点
    }
while(r-l>1&&(q[l+1].b-t[r])*(q[l+1].a-t[r])>eps)//注意删除多余元素
    --r;
t[r+1]=its(q[l+1],q[r]);//再求出新的交点
++r;
//这里不能在t里面++r需要注意一下……

练习

POJ 2451 Uyuw's Concert 注意边界

POJ 1279 Art Gallery 求多边形的核

洛谷 P4196 [CQOI2006]凸多边形

虽然上面三个都是板子,但是需要注意的细节不同。

我的计算几何做题记录

线段树总结

2023-09-07 16:28:33 By root

普通线段树数组

基于将区间拆分为$O(logn)$区间的并这一思想来实现。其中附加信息的维护方式变化很灵活。函数

区间加。利用懒标记直接维护。 区间加&区间乘。两个懒标记同时维护,其中乘法更新要带动加法更新。 区间染色。一个懒标记直接覆盖便可。 区间gcd。 树上。有点相似树剖,直接利用DFS序转化为线性,而后普通维护就行了。 维护最小值的位置。query函数的类型做为结构体,每次返回最小值及其位置。 利用取模以后缩小一半的性质,一个数$a$最多被模$loga$次。所以暴力作便可。还要须要维护最大值判断不用取模的状况。 Codeforces 292 E. Copying Data优化

给出两个序列$a,b$。动态对$b$进行修改,修改操做是指定$a$的一段区间$[x,x+k]$,粘贴到$b$的$[y,y+k]$上。单点询问$b$。spa

若是咱们知道某一个位置最后一次是被哪一次操做覆盖的,那么就能够知道这个位置对应的值。排序

由此问题转化为求解某一个位置最后一次被哪一个操做覆盖。这显然是一个区间赋值(染色)问题。利用线段树维护——只须要用懒标记强行覆盖便可。继承

Codeforces 914 D. Bash and a Tough Math Puzzle递归

给出一个序列。动态单点更新。每次询问一段区间,并给出一个数字$x$。问能不能经过修改区间内的一个数字,使区间的$gcd=x$。io

这道题告诉我gcd也是知足区间性质的。假如一段区间内全部数都是$x$的倍数,那么我能够将其中一个修改成$x$使得全部数的$gcd$为$x$;假如一段区间内只有一个数不是$x$的倍数,那么我能够将那个数修改成$x$使得全部数的$gcd$为$x$;假如一个区间内有超过一个数不是$x$倍数则无解。由此解决这个问题只须要统计区间内有多少个数不是$x$的倍数便可。扩展

咱们能够像维护区间和同样来维护区间$gcd$。若是一个区间的$gcd$是$x$的倍数,意味着这里面全是$x$的倍数;不然,继续深究,若是到达叶节点那么直接判断并返回。注意,这样的作法违背了线段树$O(logn)$查询的原则,因此须要作一些优化——一旦总数超过1就跳出。复杂度玄学。gc

Codeforces 877 E. Danil and a Part-time Job

给出一棵树,点权为bool。动态询问子树和,同时会对点权进行修改。修改的规则是:选择一棵子树$v$,子树中全部节点的点权都异或1。

用DFS序将问题转化为线性。而后考虑如何完成取反操做。很像Splay里面的那个区间翻转操做,所以咱们的标记的更新时异或1便可。而后将总和与总数做差便可。

Codeforces 383 C. Propagating tree

给出一棵树,有点权,动态查询子树和。其中会对点权进行修改,修改的规则是:选择一颗子树$v$,给$v$这个点+k,$v$的直接儿子-k,$v$的孙子+k……以此类推。

首先利用树的DFS序转化为线性问题。关键在于分层有正负的问题。虽然哪一层正哪一层负不肯定,可是相对正负关系必定是肯定的。咱们能够假设奇数层为正,偶数层为负。假设偶数层要更新为正,则更新一个负数。最后统计的时候依然按照奇数正偶数负的方式来统计就能获得正确答案。这种思想感受仍是很神奇很重要的。

树的深度统计是要放在递归以前作的,这里我写挂了。

Codeforces 438 D. The Child and Sequence

给出一个序列$a_i$($n \leq 10^5, a_i \leq 1e9$),动态询问区间和,单点修改,修改操做是区间取模。

区间取模貌似不能直接维护。事实上咱们能够暴力,加上一点小小的优化。假设咱们要让一段区间对$P$取模。那么若是区间内的最大值小于$P$,那么可忽略此操做。另外咱们发现,一个数取模之后至少缩小一半。由此$a_i$最多被模$loga_i$次。由此,复杂度近似是$O(nlognlog1e9)$。

取模以后缩小一半的性质证实:设$a%b=c$,其中$a>b>c$。可表示为$a=kb+c$,其中$k\geq 1$。由于$b>c,k \geq 1$,因此$kb>c$。因此$a>2c$

主席树

单点更新仅产生$O(logn)$个点。故复杂度为$O(mlogn)$

(一)维护区间第$k$小。经过维护区间内每一个数的出现次数,而后在线段树上二分实现。所谓区间,就是指两个历史版本的差。这里利用了前缀和的思想。

(二)维护树上路径上的第$k$小。转化为路径上每一个数的出现次数。利用树上差分的思想。

Codeforces 961 E. Tufurama

给出一个序列$a_i$($a_i \leq 1e9$),问有多少对$(x,y)$知足$a_x \geq y$且$a_y \geq x$。其中$1 \leq x

对于每个$y$,考虑$x∈[1,y-1]$。那么$x$要知足的条件有$x \leq a_y$,$a_x \geq y$。由前者结合前提可得$x∈[1,min\{y-1,a_y\}]$。

至此,问题转化为求解$x∈[1,min\{y-1,a_y\}]$里,$a_x \geq y$的$x$个数。这是一个求解区间大于$k$的个数问题,利用主席树解决便可。注意须要离散化,细节不少。

一类特殊的出现次数问题

Codeforces 813 E. Army Creation

给出一个序列$a_i$($n,a_i \leq 10^5$),动态询问区间内元素的最大个数,知足相同元素不超过$k$个。

不能维护每一个元素出现次数,这样的话很是难$log(n)$统计答案。这里有一个巧妙的转化:预处理出$a_i$以后的第$k$个元素的位置,记为$b_i$。若是$b_i \leq r$,那么就不选$a_i$;若是$b_i > r$,那么就选。这样很巧妙地限制了一段区间内咱们最多只取靠后的$k$个相同元素。所以咱们利用$b_i$创建主席树,问题转化为求解区间内大于$r$的元素个数,至关于求$(r,+\infty)$的区间和。

luogu 1972 HH的项链

给出一个序列$a_i$($n,a_i \leq 5·10^5$),动态询问区间不一样的元素个数。

容易发现这就是上面那题$k=1$的特殊状况。理论上直接主席树解决便可。

但咱们考虑一种线段树的离线作法。对于一个区间的右边界,若是咱们只统计每种元素最靠近右边界的那个位置,就不会重复了,这样区间和就是不一样元素的个数了。咱们能够这样解决:预处理出与$a_i$相等的元素上一次出现的位置$las_i$。将全部询问区间按照右端点排序,用线段树维护每一个位置,0表示不算,1表示算。右端点往右扩展时,碰到的这个数若是以前出现过,就在先前那个位置将1改成0。这样就保证了只统计最靠后的。

Codeforces 1000 F. One Occurrence

给出一个序列$a_i$($n,a_i \leq 5·10^5$),动态询问一个区间内只出现一次的数(多解可任意)。

考虑线段树离线(或主席树)。咱们一样只考虑最靠近右边界的,对于这样的$a_i$,若是上一次出现的位置$

Codeforces 833 B. Bakery

给出一个序列$a_i$($n,a_i \leq 35000$),要求把序列划分为$k$段($k \leq 50$),每段的价值为其中不一样的数的个数。求最大总价值。

明显是一个划分DP的模型。有$dp_{i,j}=max\{dp_{k,j-1}+distinct(k+1,i)\}$。直接转移须要枚举$k$,复杂度就是$O(kn^2)$级别的。所以须要考虑别的方法。

考虑已知$dp_{k,j-1}$,要考虑它对转移的贡献就是要求出全部的$distinct(k+1,i)$。考虑对于区间$(k,i]$内的一个位置$p$能不能对它产生贡献——一种数值的数最多只能产生一个贡献,不如规定最左侧的数产生贡献。这就让咱们联想到去维护一个$pre_p$,咱们只考虑$pre_p \leq k$的位置。针对这种问题,咱们让每个$p$去覆盖区间$[pre_p,p)$。因而位置$k$被覆盖了几回就表明了$distinct(k+1,i)$。用线段树维护区间修改和最大值查询便可。

总结

这类要维护区间内出现了多少个不一样的数的问题,咱们的处理方法都是:对于重复的元素,只看其中一个。或是最左侧的,或是最右侧的。判断最左最右的方法就是维护一个$pre$或者$nxt$数组。线段树维护的通常是位置,附加值表示这个位置算不算或者总共被后面覆盖几回之类的。

可持久化数组

可经过主席树来维护。数组就是主席树的叶节点,上面那些节点实际上是不存放附加信息的。

线段树合并

在树上要求解各子树的答案时,而且用线段树维护时,要合并子树的信息来获得整个树的信息,这时就要用到线段树合并。

线段树合并的思想很简单,就是合并。考虑到空间的限制,要求动态开点(时空复杂度$O(n \log n)$)

当合并根节点为$a,b$的两个线段树时,有两种代码实现的方法。一种是直接将合并后的结果放在$a$上,一种是每次合并新开一个节点做为新的根。对于题目仅仅要咱们输出一次每一个节点的答案时,可使用第一种方法,由于咱们不须要历史信息。而若是题目是在线的屡次询问,那么此时显然须要利用到历史信息的,此时就要选择方法二了。由于第一种方法咱们会在合并的时候直接继承某一个子树,这样在modify操做的时候就会直接影响到历史的线段树致使答案错误了。

Codeforces 600 E. Lomsat gelral

给出一颗树,每一个节点有一个颜色。询问每一个子树内颜色数量最多的那些颜色的总和(颜色的值)

用线段树维护最大值,再附加维护全部最大值的颜色总和。线段树合并便可。

luogu 3605 Promotion Counting

给出一颗树,每一个点有一个权值。询问每一个节点,其全部子树内的节点中权值比它大的节点个数。

离散化后,权值线段树维护和,线段树合并便可。

luogu 3899 谈笑风生

给出一棵树,每次询问一个$a$节点,问有多少个有序数对$(a,b,c)$,知足$a,b$都是$c$的祖先,且$a,b$之间距离不超过$k$

显然要用线段树来维护某一个深度的节点,关键是维护什么?咱们发现应该维护每一个深度的$size$,这样才能够直接计算答案。由于题目是在线询问的,因此每一次合并的时候要开点。

uoj配置说明

2023-02-14 15:04:11 By root

题目 只有超级管理员有新建题目的权限。如果你是超级管理员,你可以点击题库页面右下方的 “新建题目” 按钮来新建一个题目。然后,你就可以在题目后台对题目进行管理。后台分为三个选项卡:

编辑:题面编辑页面 管理者:题目管理员列表管理页面 数据:题目数据管理页面 如果你不是超级管理员,请联系超级管理员帮忙添加题目,然后帮你加上管理该题目的权限。

下面我们将逐一介绍这三个选项卡(当然还有最后个 “返回” 选项卡退出后台,就不介绍了),然后我们将详细介绍题目数据需要满足的格式要求。

一、选项卡:编辑 该页面可以对题面本身进行编辑,还可以管理题目的标签。

▶ 编辑题面 题面用 Markdown 编写。

理论上题面是可以自由编写的,但还是有一些推荐的格式。可参照 UOJ 上对应类型的题面(传统题、交互题、提交答案题等)

一些推荐的规则:

中文与英文、数字之间加一个空格隔开。 输入输出样例用

 标签包围,并且以一个空行结尾。(方便大家复制样例到终端后不用再打回车)
题面中最高级标题为三级标题。这是因为题目名称是二级标题。
名称超过一个字符的数学符号要用 mathrm。例如 \mathrm{sgn}, \mathrm{lca}。
注意 \max, \min, \gcd 都是有现成的。
注意 xor 这种名称超过一个字符的二元运算符请用 \mathbin{\mathrm{xor}}。
一行内的计算机输入输出和常量字符串表达式请用  标签,例如 请输出 “NO SOLUTION”, aaa 中有三个 a。
一行内的计算机代码请用 ` 括起来,就像上面的规则那样。
可参考下面这个例子:

读入一个整数 $n$,表示题目中提到的 $n$ 位大爷 AC 的总题数。 请分别输出他们 AC 的总题数。

如果你不能正确读入,那么你将获得 $0$ 分。

前 $3$ 个测试点你正确读入即可获得 $6$ 分,第 4 个测试点你正确读入只能获得 $3$ 分。

如果你不会做这道题,请直接输出 “WOBUHUI”。

下面是一个样例:

233

▶ 编辑标签 只需要用英文逗号隔开标签填入文本框就行。

推荐你使用如下几条规则填写标签:

标签的目的是标出题目类型,方便用户检索题目。一般来说,标签顺序基本为从小范围到大范围。 最前面的几个标签是这题所需要的前置技能,这里假定 “二分查找” 之类过于基础的技能选手已经掌握。 接下来是这道题的大方法,比如 “贪心”、“DP”、“乱搞”、“构造”、“分治”…… 接下来,如果这道题是非传统题,用一个标签注明非传统题类型,比如 “提交答案”、“交互式”、“通讯”。 接下来,如果这道题是模板题,用一个标签注明 “模板题”。 接下来,如果这道题是不用脑子想就能做出的题,例如 NOIP 第一题难度,用一个标签注明 “水题”。 最后,如果这题的来源比较重要,用一个标签注明。比如 “UOJ Round”、“NOI”、“WC”。 前置技能中,“数学” 太过宽泛不能作为标签,但 “数论” 可以作为前置技能。 如果有多个解法,每个解法的前置技能和大方法都不太一样,那么尽可能都标上去。 “乱搞” 标签不宜滥用。 二、标签页:管理者 可以通过这个页面增加或删除题目管理员。+mike 表示把 mike 加入题目管理员列表,-mike 表示把 mike 从题目管理员列表中移除。

题目管理员和超级管理员都能看到题目后台,但只有题目管理员才能通过 SVN 管理题目数据。所以,必要时可以将超级管理员也添加为题目管理员。当然,除了通过 SVN 管理题目数据之外,直接通过网页浏览器管理也是可以的。

三、标签页:数据 UOJ 使用 SVN 管理题目数据,SVN 仓库地址可以在本页面上找到。配置题目一般分两步,第一步是上传题目数据,第二步是让 UOJ 根据上传来的数据,生成一份供测评使用的数据包。下面我们将依次进行介绍。

▶ 上传原始数据 上传数据有三种方式。第一种方法是点击 “浏览 SVN 仓库” 进行上传。点击后会进入到一个文件浏览器界面,然后再点击右上角上传,就可以把本地文件直接拖拽过来上传了。

第二种方法是点击 “直接上传一个压缩包” 进行上传。点击后会弹出一个对话框,然后你就可以直接上传一个包含题目数据的 zip 压缩包了。上传后你可以点击 “浏览 SVN 仓库” 检查上传内容是否正确。

第三种方式是通过本地 SVN 客户端上传。这个有点复杂,你需要了解一点 SVN 的基础知识。不过你如果已经学会了使用 git 的话,SVN 是不难掌握的。

题外话:在我们开发 UOJ 核心代码的 2014 年,SVN 和 git 似乎还难分伯仲。当时考虑到 Codeforces 附属的造题网站 Polygon 使用的是 SVN,所以我们最后让 UOJ 使用的也是 SVN 而非大家可能更熟知的 git。

想要使用 SVN,你得先获取你 UOJ 用户所对应的 SVN 密码。当你点击了按钮 “我会用 SVN 客户端,请发我 SVN 账号密码” 之后,UOJ 就会把你的 SVN 账号密码发送到你的邮箱。这个 SVN 密码是不随题目变化而变化的,所以你只用获取一次,然后好好保存。

但是,初始时 UOJ 是没有配置过用于给大家发邮件的邮箱的。所以如果你还没有配置过,请直接在数据库的 user_info 表中找到 svn_password 一栏存储的值;如果你愿意折腾一下,可以在 /var/www/uoj/app/.config.php 里配置一下邮箱,详情见安装教程咯~

下面是使用 SVN 上传的一点指南:

下载 SVN,即 subversion。如果是 win 可以用 TortoiseSVN,如果是 ubuntu 可以直接 apt install,也可以使用 smartsvn。 学习怎么使用 SVN,学到大约会用 checkout、add、remove、commit 的程度就可以了。 在题目数据管理页面,获取 svn 密码。 根据“数据”标签页上的 SVN 仓库地址,checkout 这道题。 在你的 working copy 下,建立名为 1 的文件夹,然后在 1 这个文件夹下放置题目数据。这一点非常重要,如果不放在该文件夹下,UOJ 将无法识别。 搞完之后 commit。 在网页上点“与SVN仓库同步”,就可以啦! 将 UOJ 设计为数据都放在文件夹 1 下初衷是为了方便大家存同一题目的不同版本(其中主版本叫 1),但这样确实给传题的管理员们造成了很多困惑。所以如果网速可观,推荐大家点击 “浏览 SVN 仓库” 上传数据,这样你就不用考虑文件夹 1 的问题了,直接将数据放在根目录即可。

▶ 生成测评数据包 UOJ 这一部分的交互界面做得有点晦涩。。

上传了原始数据之后,你可以点击右侧按钮 “与svn仓库同步” 来让 UOJ 对你上传的原始数据进行格式检查,并最终生成测评时使用的数据包。一般而言,这一步骤主要完成的是如下操作:

检查原始数据中的输入输出文件里的换行符。如果是 \r\n,会被自动替换为 \n。多余的行末空格也会被替换掉。 编译测评所需的文件,比如答案检查器(chk)、标准程序(std)、数据检验器(val)。 将 download 文件夹下的文件打包成压缩包,作为该题的附件供选手下载。 总之,你需要在上传数据之后点击 “与svn仓库同步”,检查是否有报错。顺利生成测评数据包之后,这道题的测评功能才会正常工作。

另外还要注意 Hack 功能的开关。初始时,题目是允许 Hack 的。这意味着你只有在上传 std 和 val 后,“与svn仓库同步” 才会停止报错。但如果你不想写,或者暂时不打算写,你可以点击 “禁止使用 Hack” 的按钮来生成测评数据包。

四、题目数据格式 下面我们先从传统题的配置讲起,然后再依次展开介绍各类非传统题的配置方式。

▶ 传统题的数据格式 1. 数据配置文件 假设你要上传一道传统题,输入文件是 www1.in, www2.in, ... 这个样子,输出文件 www1.out, www2.out, ...。你想设时间限制为 1s,空间限制为 256MB,那么你需要在上传这些输入输出文件的同时,编写一个类似下面这样的数据配置文件 problem.conf:

use_builtin_judger on use_builtin_checker ncmp n_tests 10 n_ex_tests 5 n_sample_tests 1 input_pre www input_suf in output_pre www output_suf out time_limit 1 memory_limit 256 output_limit 64 当然你也可以通过点击 “修改数据配置” 按钮,来在线修改该数据配置文件的内容。

测评类型:use_builtin_judger on 这一行指定了该题使用 UOJ 内置的测评器进行测评,绝大多数情况下都只需要使用内置的测评器。后面我们将介绍如何使用自定义的测评器。

输入输出文件:UOJ 会根据 printf("%s%d.%s", input_pre, i, input_suf) 这种方式来匹配第 i 个输入文件,你可以将 input_pre 和 input_suf 改成自己想要的值。输出文件同理。

额外数据:extra test 是指额外数据,在 AC 的情况下会测额外数据,如果某个额外数据通不过会被倒扣3分。额外数据中,输入文件命名形如 ex_www1.in, ex_www2.in, ... 这个样子,一般来说就是 printf("ex_%s%d.%s", input_pre, i, input_suf),输出文件同理。

样例数据:样例必须被设置为前若干个额外数据。所以数据配置文件里有一个 n_sample_tests 参数,表示的是前多少组是样例。你需要把题面中的样例和大样例统统放进额外数据里。

答案检查器(chk):数据配置文件中,use_builtin_checker ncmp 的意思是将本题的 checker 设置为 ncmp。checker 是指判断选手输出是否正确的答案检查器。一般来说,如果输出结果为整数序列,那么用 ncmp 就够了。ncmp 会比较标准答案的整数序列和选手输出的整数序列。如果是忽略所有空白字符,进行字符串序列的比较,可以用 wcmp。如果你想按行比较(不忽略行末空格,但忽略文末回车),可以使用 fcmp。

如果想使用自定义的 checker,请使用 Codeforces 的 testlib 库编写 checker,并命名为 chk.cpp,然后在 problem.conf 中去掉 “use_builtin_checker” 这一行。

时空限制:time_limit 控制的是一个测试点的时间限制,单位为秒,可以是小数(至多三位小数,即最高精确到毫秒)。

memory_limit 控制的是一个测试点的空间限制,单位为 MB。output_limit 控制的是程序的输出长度限制,单位也为 MB。注意这两个限制都不能为小数。

  1. 配置 Hack 功能 如果想让你的题目支持 Hack 功能,那么你还要额外上传点程序。

数据检验器(val):数据检验器的功能是检验你的题目的输入数据是否满足题目要求。比如你输入保证是个连通图,validator 就得检验这个图是不是连通的。但比如 “保证数据随机” “保证数据均为人手工输入” 等就无法进行检验。请使用 Codeforces 的 testlib 库编写,写好后命名为 val.cpp 即可。

标准程序(std):标准程序的功能是给定输入,生成标准答案。默认时空限制与选手程序相同。Hack 时,chk 将会根据 std 的输出,判断选手的输出是否正确。

  1. 辅助测评程序 我们暂且将 chk、val、std 这种称为辅助测评程序。上面我们已经介绍了这些辅助测评程序各自的用途,下面我们介绍两个有用的配置。

配置语言:默认情况下,类似 std.cpp 的文件名将会被识别为 C++ 03(即编写 UOJ 核心代码时 g++ 的默认 C++ 语言标准)。如果想使用 C++ 11、C++ 14、C++ 17、C++ 20,请将文件名改为 std11.cpp、std14.cpp、std17.cpp、std20.cpp。类似的也可以有 chk14.cpp、val20.cpp。下面是文件后缀名到语言的映射表(详见 app/models/UOJLang.php):

<?php $suffix_map = [ '20.cpp' => 'C++20', '17.cpp' => 'C++17', '14.cpp' => 'C++14', '11.cpp' => 'C++11', '.cpp' => 'C++', // C++ 03 '.c' => 'C', '.pas' => 'Pascal', '2.7.py' => 'Python2.7', '.py' => 'Python3', '7.java' => 'Java7', '8.java' => 'Java8', '11.java' => 'Java11', '14.java' => 'Java14', ]; 配置时空限制:默认情况下,chk 和 val 的时空限制为 5s + 256MB,std 的时空限制与选手程序相同。如果你想让这些辅助测评程序拥有更宽松的时空限制,你可以在 problem.conf 中加入

_time_limit 100

_memory_limit 1024 来将其时间限制设置为 100s,空间限制设置为 1024MB。其中,chk 对应的 为 checker,val 对应的 为 validator,std 对应的 为 standard。 4. 测试点分值 默认情况下,如果测试点全部通过,那么直接给 100 分。否则,如果你有 n 个测试点,那么每个测试点的分值为 floor(100/n)。

如果你想改变其中某个测试点的分值,可以在 problem.conf 里加入类似 point_score_5 20 这样的配置(意为把第 5 个测试点的分值改为 20)。

UOJ 支持出题人设置测试点的分值的数值类型。共有两种:整数和实数。在整数模式中,每个测试点的总分值必须为整数,实际得分也为整数;在实数模式中,测试点的分值可以有小数,最后实际得分会用浮点数进行运算,四舍五入到小数点后若干位。默认模式为整数模式,也可手动设置:score_type int;实数模式设置方法形如 score_type real-3,意为实数模式 + 保留小数点后 3 位。可以通过改变后面的这个 “3” 来控制保留到小数点后几位(可以设为 0 到 10 之间的任意整数)。

有些题目(特别是提答题)会有给单个测试点打部分分的需求。此时请使用自定义答案检查器 chk,并使用 quitp 函数。在整数模式下,假设测试点满分为 ,则

quitp(x, "haha"); 将会给 floor(a round(100 x) / 100) 分。即,先将 四舍五入到小数点后两位,得到 ,然后给分 。这样设计是因为 UOJ 期望你输出的 是一个两位小数。当你这个测试点想给 分( 为整数)的时候,你输出的两位小数 需要满足

因此,最稳妥的办法是:

quitp(ceil(100.0 * p / a) / 100, "haha"); 虽然这个上取整有点反直觉。 为什么整数模式下,得分是下取整的?因为下取整才可以保证答案有误的时候得分不为满分。如果想要四舍五入,可以使用实数模式。在实数模式下,若设置的精度为 ,那么 UOJ 期望你精确地输出 ,而最后得分会是将 四舍五入到小数点后第 位所得到的结果。因此,当你想给 分( 可以是小数),那么直接像下面这样写就好了:

quitp((double)p / a, "haha"); 如果你的题目是一道以子任务形式评分的题目,那么你可以用如下语法划分子任务和设定分值:

n_subtasks 6 subtask_end_1 5 subtask_score_1 10 subtask_end_2 10 subtask_score_2 10 subtask_end_3 15 subtask_score_3 10 subtask_end_4 20 subtask_score_4 20 subtask_end_5 25 subtask_score_5 20 subtask_end_6 40 subtask_score_6 30 可以用 subtask_type_5 将第 5 个子任务的评分类型设置为 packed 或 min。其中 packed 表示错一个就零分,min 表示测评子任务内所有测试点并取得分的最小值。默认值为 packed。

可以用 subtask_used_time_type_5 将第 5 个子任务统计程序用时的方式设置为 sum 或 max。其中 sum 表示将子任务内测试点的用时加起来,min 表示取子任务内测试点的用时最大值。默认值为 sum。

可以用 subtask_dependence_5 3 来将第 3 个子任务添加为第 5 个子任务的依赖任务,相当于测评时第 5 个子任务会拥有第 3 个子任务的所有测试点。所以如果第 3 个子任务中有一个错了,并且第 5 个子任务的评分类型为 packed,那么第 5 个子任务也会自动零分。如果你想给一个子任务添加多个依赖任务,比如同时依赖 3 和 4,那么你可以用如下语法:

subtask_dependence_5 many subtask_dependence_5_1 3 subtask_dependence_5_2 4 5. 目录结构 下面是存储一道传统题数据的文件夹所可能具有的结构:

www1.in, www2.in, ... :输入文件 www1.out, www2.out, ... :输出文件 ex_www1.in, ex_www2.in, ... :额外输入文件 ex_www1.out, ex_www2.out, ... :额外输出文件 problem.conf:数据配置文件 chk.cpp:答案检查器 val.cpp:数据检验器 std.cpp:标准程序 download/:一个文件夹,该文件夹下的文件会被打包成压缩包,作为该题的附件供选手下载。 require/:一个文件夹,该文件夹下的文件会被复制到选手程序运行时所在的工作目录。 注意:如果你是手动通过 SVN 客户端操作,那么你需要把这些文件都放在名为 1/ 的文件夹下。

▶ 提交答案题的配置 范例:

use_builtin_judger on submit_answer on n_tests 10 input_pre www input_suf in output_pre www output_suf out So easy!

不过还没完,因为一般来说你还要写个 chk。当然如果你是道 最小割计数 就不用写了,直接用 use_builtin_checker ncmp 就行。

假设你已经写好 chk 了,你会发现你还需要发给选手一个本地可以运行的版本,这怎么办呢?如果只是传参进来一个测试点编号,检查测试点是否正确,那么只要善用 registerTestlib 就行了。如果你想让 checker 与用户交互,请使用 registerInteraction。特别地,如果你想让 checker 通过终端与用户进行交互,请使用如下代码(有点丑啊)

char *targv[] = {argv[0], "inf_file.out", (char*)"stdout"};

ifdef __EMSCRIPTEN__

setvbuf(stdin, NULL, _IONBF, 0);

endif

registerInteraction(3, targv);

那么怎么下发文件呢?你可以把编译好的 chk 放在文件夹 download 下,该文件夹下的所有文件都会被打包发给选手。

交叉编译小技巧:(本段写于 2016 年)由于你的本地 checker 可能需要发布到各个平台,所以你也许需要交叉编译。交叉编译的意思是在一个平台上编译出其他平台上的可执行文件。好!下面是 Ubuntu 上交叉编译小教程时间。首先不管三七二十一装一些库:(请注意 libc 和 libc++ 可能有更新的版本)

sudo apt-get install libc6-dev-i386 sudo apt-get install lib32stdc++6 lib32stdc++-4.8-dev sudo apt-get install mingw32 然后就能编译啦:

g++ localchk.cpp -o checker_linux64 -O2 g++ localchk.cpp -m32 -o checker_linux32 -O2 i586-mingw32msvc-g++ localchk.cpp -o checker_win32.exe -O2 问:那其他系统怎么办?有个叫 emscripten 的东西,可以把 C++ 代码编译成 javascript。请务必使用 UOJ 上的 testlib,而不是 CF 上的那个版本。由于 testlib 有些黑科技,UOJ 上的 testlib 对 emscripten 才是兼容的。编译就用

em++ localchk.cpp -o checker.js -O2 ▶ 交互题的配置 UOJ 内部并没有显式地支持交互式,而是提供了 require、implementer 和 token 这几个东西。

  1. require 文件夹 除了前文所述的 download 这个文件夹,还有一个叫 require 的神奇文件夹。在测评时,这个文件夹里的所有文件均会被移动到与选手源程序同一目录下。在这个文件夹下,你可以放置交互库,供编译时使用。

另外,require 文件夹的内容显然不会被下发……如果你想下发一个样例交互库,可以自己放到 download 文件夹里。

  1. implementer 代码 再来说 implementer 的用途。如果你在 problem.conf 里设置 with_implementer on 的话,各个语言的编译命令将变为:

C++: g++ code implementer.cpp code.cpp -lm -O2 -DONLINE_JUDGE C: gcc code implementer.c code.c -lm -O2 -DONLINE_JUDGE Pascal: fpc implementer.pas -o code -O2 这样,一个交互题就大致做好了。你需要写 implementer.cpp、implementer.c 和 implementer.pas。当然你不想支持某个语言的话,就可以不写对应的交互库。

如果是 C/C++,正确姿势是搞一个统一的头文件放置接口,让 implementer 和选手程序都 include 这个头文件。把主函数写在 implementer 里,连起来编译时程序就是从 implementer 开始执行的了。

如果是 Pascal,正确姿势是让选手写一个 pascal unit 上来,在 implementer.pas 里 uses 一下。除此之外也要搞另外一个 pascal unit,里面存交互库的各种接口,让 implementer 和选手程序都 uses 一下。

  1. token 机制 下面我们来解释什么是 token。考虑一道典型的交互题,交互库跟选手函数交流得很愉快,最后给了个满分。此时交互库输出了 “AC!!!” 的字样,也可能输出了选手一共调用了几次接口。但是聪明的选手发现,只要自己手动输出一下 “AC!!!” 然后果断退出似乎就能骗过测评系统了哈哈哈。

这是因为选手函数和交互库已经融为一体成为了一个统一的程序,测评系统分不出谁是谁。这时候,token 就出来拯救世界了。

在 problem.conf 里配置一个奇妙的 token,比如 token wahaha233,然后把这个 token 写进交互库的代码里(线上交互库,不是样例交互库)。在交互库准备输出任何东西之前,先输出一下这个 token。当测评系统要给选手判分时,首先判断文件的第一行是不是 token,如果不是,直接零分。这样就避免了奇怪的选手 hack 测评系统了。这里的 token 判定是在调用 checker 之前,测评系统会把 token 去掉之后的选手输出喂给 checker。(所以一道交互题的 token 要好好保护好)

不过这样的防御手段并不是绝对的。理论上只要选手能编写一个能理解交互库的机器码的程序,试图自动模拟交互库行为,那是很难防住的。

另外,记得在 C/C++ 的交互库里给全局变量开 static,意为仅本文件里的代码可访问(即仅交互库可访问)。

▶ 简单通信题的配置 这里我们考虑最简单的通信题。有一个交互器 interactor,选手程序的输出会变成 interactor 的输入,而 interactor 的输出会变成选手程序的输入。

(有时这种题目也叫交互题。。但为了不致混淆我们这里称之为通信题。)

如果你想配置这样的题目,只需要在 problem.conf 中加一行

interaction_mode on 然后编写程序 interactor.cpp,跟输入输出文件放在一起即可。interactor 也属于辅助测评程序,可以参照 “辅助测评程序” 这一小节的内容进行配置。比如,通过添加一行 interactor_time_limit 2 可将 interactor 的时间限制设为 2s。 UOJ 的交互器 interactor 与 Codeforces 的交互器有一点点不同。在 Codeforces 的测评逻辑里,交互器与选手程序交互后,会把一些信息输出到文件,然后再交由答案检查器 chk 作出最终评分。但是 UOJ 里的交互器相当于 interactor + chk,也就是说交互器需要直接给出最终评分。所以,当你在使用 testlib 库编写交互器时,需要按照编写答案检查器的方式,在 main 函数里第一句写上

registerTestlibCmd(argc, argv); 之后你就可以通过 ouf 获取选手程序输出,通过 cout 或者 printf 向选手程序输出信息(不要忘了刷新缓冲区),通过 quitf 或者 quitp 打分。注意,不要按照 Codeforces 的习惯调用 registerInteraction。 如果你想实现更复杂一点通信,比如多个程序相互之间通信,那么你可能需要参考下一节 “意想不到的非传统题” 的内容自己编写 judger 了。

▶ 意想不到的非传统题 假如你有一道题不属于上述任何一种题目类型 —— 那你就得自己写测评器了!

测评器(judger)的任务是:给你选手提交的文件,输出一个测评结果。

把 problem.conf 里的 use_builtin_judger on 这行去掉,放一个 judger 程序的源代码放在题目数据目录下,就好了。

但怎么编写 judger 呢。。。?首先你需要写个 Makefile,UOJ 会根据你的 Makefile 和 judger 源代码生成 judger。如果你要使用自定义答案检查器 chk,那么你也需要在 Makefile 里面指定 chk 的编译方式。比如下面这个例子:

export INCLUDE_PATH CXXFLAGS = -I$(INCLUDE_PATH) -O2

all: chk judger

% : %.cpp $(CXX) $(CXXFLAGS) $< -o $@ 什么?你问 judger 应该怎么写?

这部分文档我只能先【坑着】。不过这里有一点有助于你自己看懂代码的指南:

你可以先从 UOJ 内置的测评器读起,在 uoj_judger/builtin/judger/judger.cpp; 你会发现它引用了头文件 uoj_judger.h,大概读一读; uoj_judger.h 调用了沙箱 uoj_judger/run/run_program.cpp,大概理解一下沙箱的参数和输入输出方式; uoj_judger.h 还会用 uoj_judger/run/compile.cpp 来编译程序,用 uoj_judger/run/run_interaction.cpp 来处理程序间的通信,大概读一读; 回过头来仔细理解一下 uoj_judger.h 的代码逻辑。 不过说实话 uoj_judger.h 封装得并不是很好。所以有计划重写一个封装得更好的库 uoj_judger_v2.h,这个库拖拖拉拉写了一半,但暂且还是先放了出来,供大家参考。期待神犇来造一造轮子!

最后,如果要使用自定义 judger 的话,就只有超级管理员能点击 “与 SVN 仓库同步” 了。这是因为 judger 的行为几乎不受限制,只能交给值得信赖的人编写。不过无论你的 judger 如何折腾,测评也是有 10 分钟的时间限制的。

四、样题 在这个仓库里我们公开了一些来自 UOJ 官网的题目数据,供大家参考。

计算几何学习

2022-11-29 21:31:23 By root

画图工具
计算几何资料

一、引言
  计算机的出现使得很多原本十分繁琐的工作得以大幅度简化,但是也有一些在人们直观看来很容易的问题却需要拿出一套并不简单的通用解决方案,比如几何问题。作为计算机科学的一个分支,计算几何主要研究解决几何问题的算法。在现代工程和数学领域,计算几何在图形学、机器人技术、超大规模集成电路设计和统计等诸多领域有着十分重要的应用。在本文中,我们将对计算几何常用的基本算法做一个全面的介绍,希望对您了解并应用计算几何的知识解决问题起到帮助。

二、目录
本文整理的计算几何基本概念和常用算法包括如下内容:
1. 矢量的概念
2. 矢量加减法
3. 矢量叉积
4. 折线段的拐向判断
5. 判断点是否在线段上
6. 判断两线段是否相交
7. 判断线段和直线是否相交
8. 判断矩形是否包含点
9. 判断线段、折线、多边形是否在矩形中
10. 判断矩形是否在矩形中
11. 判断圆是否在矩形中
12. 判断点是否在多边形中
13. 判断线段是否在多边形内
14. 判断折线是否在多边形内
15. 判断多边形是否在多边形内
16. 判断矩形是否在多边形内
17. 判断圆是否在多边形内
18. 判断点是否在圆内
19. 判断线段、折线、矩形、多边形是否在圆内
20. 判断圆是否在圆内
21. 计算点到线段的最近点
22. 计算点到折线、矩形、多边形的最近点
23. 计算点到圆的最近距离及交点坐标
24. 计算两条共线的线段的交点
25. 计算线段或直线与线段的交点
26. 求线段或直线与折线、矩形、多边形的交点
27. 求线段或直线与圆的交点
28. 凸包的概念
29. 凸包的求法
30. 半平面交
31. 辛普森积分
32. 圆的反演
三、算法介绍
1.矢量的概念:

  如果一条线段的端点是有次序之分的,我们把这种线段成为有向线段(directed segment)。如果有向线段p1p2的起点p1在坐标原点,我们可以把它称为矢量(vector)p2。
2.矢量加减法:
   设二维矢量P = ( x1, y1 ),Q = ( x2 , y2 ),则矢量加法定义为: P + Q = ( x1 + x2 , y1 + y2 ),同样的,矢量减法定义为: P - Q = ( x1 - x2 , y1 - y2 )。显然有性质 P + Q = Q + P,P - Q = - ( Q - P )。
3.  矢量叉积:
  计算矢量叉积是与直线和线段相关算法的核心部分。设矢量P = ( x1, y1 ),Q = ( x2, y2 ),则矢量叉积定义为由(0,0)、p1、p2和p1+p2所组成的平行四边形的带符号的面积,即:P × Q = x1*y2 - x2*y1,其结果是一个标量。显然有性质 P × Q = - ( Q × P ) 和 P × ( - Q ) = - ( P × Q )。一般在不加说明的情况下,本文下述算法中所有的点都看作矢量,两点的加减法就是矢量相加减,而点的乘法则看作矢量叉积。
  叉积的一个非常重要性质是可以通过它的符号判断两矢量相互之间的顺逆时针关系:
  若 P × Q > 0 , 则P在Q的顺时针方向。
  若 P × Q < 0 , 则P在Q的逆时针方向。
  若 P × Q = 0 , 则P与Q共线,但可能同向也可能反向。
4.  折线段的拐向判断:
  折线段的拐向判断方法可以直接由矢量叉积的性质推出。对于有公共端点的线段p0p1和p1p2,通过计算(p2 - p0) × (p1 - p0)的符号便可以确定折线段的拐向:
  若(p2 - p0) × (p1 - p0) > 0,则p0p1在p1点拐向右侧后得到p1p2。
  若(p2 - p0) × (p1 - p0) < 0,则p0p1在p1点拐向左侧后得到p1p2。
  若(p2 - p0) × (p1 - p0) = 0,则p0、p1、p2三点共线。
具体情况可参考下图:
5.  判断点是否在线段上:
  设点为Q,线段为P1P2 ,判断点Q在该线段上的依据是:( Q - P1 ) × ( P2 - P1 ) = 0 且 Q 在以 P1,P2为对角顶点的矩形内。前者保证Q点在直线P1P2上,后者是保证Q点不在线段P1P2的延长线或反向延长线上,对于这一步骤的判断可以用以下过程实现:
  ON-SEGMENT(pi,pj,pk)
  if min(xi,xj) <= xk <= max(xi,xj) and min(yi,yj) <= yk <= max(yi,yj)
  then return true;
  else return false;
  特别要注意的是,由于需要考虑水平线段和垂直线段两种特殊情况,min(xi,xj)<=xk<=max(xi,xj)和min(yi,yj)<=yk<=max(yi,yj)两个条件必须同时满足才能返回真值。
6.  判断两线段是否相交:
  我们分两步确定两条线段是否相交:
  (1)快速排斥试验
    设以线段 P1P2 为对角线的矩形为R, 设以线段 Q1Q2 为对角线的矩形为T,如果R和T不相交,显然两线段不会相交。
  (2)跨立试验
    如果两线段相交,则两线段必然相互跨立对方。若P1P2跨立Q1Q2 ,则矢量 ( P1 - Q1 ) 和( P2 - Q1 )位于矢量( Q2 - Q1 ) 的两侧,即( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) < 0。上式可改写成( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) > 0。当 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 时,说明 ( P1 - Q1 ) 和 ( Q2 - Q1 )共线,但是因为已经通过快速排斥试验,所以 P1 一定在线段 Q1Q2上;同理,( Q2 - Q1 ) ×(P2 - Q1 ) = 0 说明 P2 一定在线段 Q1Q2上。所以判断P1P2跨立Q1Q2的依据是:( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。同理判断Q1Q2跨立P1P2的依据是:( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。具体情况如下图所示:
  在相同的原理下,对此算法的具体的实现细节可能会与此有所不同,除了这种过程外,大家也可以参考《算法导论》上的实现。
7.  判断线段和直线是否相交:
  有了上面的基础,这个算法就很容易了。如果线段P1P2和直线Q1Q2相交,则P1P2跨立Q1Q2,即:( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。
8.  判断矩形是否包含点:
  只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。
9.  判断线段、折线、多边形是否在矩形中:
  因为矩形是个凸集,所以只要判断所有端点是否都在矩形中就可以了。
10.  判断矩形是否在矩形中:
  只要比较左右边界和上下边界就可以了。
11.  判断圆是否在矩形中:
  很容易证明,圆在矩形中的充要条件是:圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。
12.  判断点是否在多边形中:
  判断点P是否在多边形中是计算几何中一个非常基本但是十分重要的算法。以点P为端点,向左方作射线L,由于多边形是有界的,所以射线L的左端一定在多边形外,考虑沿着L从无穷远处开始自左向右移动,遇到和多边形的第一个交点的时候,进入到了多边形的内部,遇到第二个交点的时候,离开了多边形,……所以很容易看出当L和多边形的交点数目C是奇数的时候,P在多边形内,是偶数的话P在多边形外。
  但是有些特殊情况要加以考虑。如图下图(a)(b)(c)(d)所示。在图(a)中,L和多边形的顶点相交,这时候交点只能计算一个;在图(b)中,L和多边形顶点的交点不应被计算;在图(c)和(d) 中,L和多边形的一条边重合,这条边应该被忽略不计。如果L和多边形的一条边重合,这条边应该被忽略不计。
  为了统一起见,我们在计算射线L和多边形的交点的时候,1。对于多边形的水平边不作考虑;2。对于多边形的顶点和L相交的情况,如果该顶点是其所属的边上纵坐标较大的顶点,则计数,否则忽略;3。对于P在多边形边上的情形,直接可判断P属于多边行。由此得出算法的伪代码如下:
count ← 0;
以P为端点,作从右向左的射线L;
for 多边形的每条边s
do if P在边s上
then return true;
if s不是水平的
then if s的一个端点在L上
if 该端点是s两端点中纵坐标较大的端点
then count ← count+1
else if s和L相交
then count ← count+1;
if count mod 2 = 1
then return true;
else return false;
  其中做射线L的方法是:设P'的纵坐标和P相同,横坐标为正无穷大(很大的一个正数),则P和P'就确定了射线L。
  判断点是否在多边形中的这个算法的时间复杂度为O(n)。
  另外还有一种算法是用带符号的三角形面积之和与多边形面积进行比较,这种算法由于使用浮点数运算所以会带来一定误差,不推荐大家使用。
13.  判断线段是否在多边形内:
  线段在多边形内的一个必要条件是线段的两个端点都在多边形内,但由于多边形可能为凹,所以这不能成为判断的充分条件。如果线段和多边形的某条边内交(两线段内交是指两线段相交且交点不在两线段的端点),因为多边形的边的左右两侧分属多边形内外不同部分,所以线段一定会有一部分在多边形外(见图a)。于是我们得到线段在多边形内的第二个必要条件:线段和多边形的所有边都不内交。
  线段和多边形交于线段的两端点并不会影响线段是否在多边形内;但是如果多边形的某个顶点和线段相交,还必须判断两相邻交点之间的线段是否包含于多边形内部(反例见图b)。
  因此我们可以先求出所有和线段相交的多边形的顶点,然后按照X-Y坐标排序(X坐标小的排在前面,对于X坐标相同的点,Y坐标小的排在前面,这种排序准则也是为了保证水平和垂直情况的判断正确),这样相邻的两个点就是在线段上相邻的两交点,如果任意相邻两点的中点也在多边形内,则该线段一定在多边形内。
  证明如下:
  命题1:
    如果线段和多边形的两相邻交点P1 ,P2的中点P' 也在多边形内,则P1, P2之间的所有点都在多边形内。
  证明:
    假设P1,P2之间含有不在多边形内的点,不妨设该点为Q,在P1, P'之间,因为多边形是闭合曲线,所以其内外部之间有界,而P1属于多边行内部,Q属于多边性外部,P'属于多边性内部,P1-Q-P'完全连续,所以P1Q和QP'一定跨越多边形的边界,因此在P1,P'之间至少还有两个该线段和多边形的交点,这和P1P2是相邻两交点矛盾,故命题成立。证毕。
  由命题1直接可得出推论:
  推论2:
  设多边形和线段PQ的交点依次为P1,P2,……Pn,其中Pi和Pi+1是相邻两交点,线段PQ在多边形内的充要条件是:P,Q在多边形内且对于i =1, 2,……, n-1,Pi ,Pi+1的中点也在多边形内。
  在实际编程中,没有必要计算所有的交点,首先应判断线段和多边形的边是否内交,倘若线段和多边形的某条边内交则线段一定在多边形外;如果线段和多边形的每一条边都不内交,则线段和多边形的交点一定是线段的端点或者多边形的顶点,只要判断点是否在线段上就可以了。
  至此我们得出算法如下:
if 线端PQ的端点不都在多边形内
then return false;
点集pointSet初始化为空;
for 多边形的每条边s
do if 线段的某个端点在s上
then 将该端点加入pointSet;
else if s的某个端点在线段PQ上
then 将该端点加入pointSet;
else if s和线段PQ相交 // 这时候已经可以肯定是内交了
then return false;
将pointSet中的点按照X-Y坐标排序;
for pointSet中每两个相邻点 pointSeti ] , pointSeti+1]
do if pointSet i] , pointSeti+1] 的中点不在多边形中
then return false;
return true;
  这个过程中的排序因为交点数目肯定远小于多边形的顶点数目n,所以最多是常数级的复杂度,几乎可以忽略不计。因此算法的时间复杂度也是O(n)。

14.  判断折线是否在多边形内:
  只要判断折线的每条线段是否都在多边形内即可。设折线有m条线段,多边形有n个顶点,则该算法的时间复杂度为O(m*n)。
15.  判断多边形是否在多边形内:
  只要判断多边形的每条边是否都在多边形内即可。判断一个有m个顶点的多边形是否在一个有n个顶点的多边形内复杂度为O(m*n)。
16.  判断矩形是否在多边形内:
  将矩形转化为多边形,然后再判断是否在多边形内。
17  判断圆是否在多边形内:
  只要计算圆心到多边形的每条边的最短距离,如果该距离大于等于圆半径则该圆在多边形内。计算圆心到多边形每条边最短距离的算法在后文阐述。
18.  判断点是否在圆内:
  计算圆心到该点的距离,如果小于等于半径则该点在圆内。
19.  判断线段、折线、矩形、多边形是否在圆内:
  因为圆是凸集,所以只要判断是否每个顶点都在圆内即可。
20.  判断圆是否在圆内:
  设两圆为O1,O2,半径分别为r1, r2,要判断O2是否在O1内。先比较r1,r2的大小,如果r1<r2则O2不可能在O1内;否则如果两圆心的距离大于r1 - r2 ,则O2不在O1内;否则O2在O1内。
21.  计算点到线段的最近点:
  如果该线段平行于X轴(Y轴),则过点point作该线段所在直线的垂线,垂足很容易求得,然后计算出垂足,如果垂足在线段上则返回垂足,否则返回离垂足近的端点;如果该线段不平行于X轴也不平行于Y轴,则斜率存在且不为0。设线段的两端点为pt1和pt2,斜率为:k = ( pt2.y - pt1. y ) / (pt2.x - pt1.x );该直线方程为:y = k* ( x - pt1.x) + pt1.y。其垂线的斜率为 - 1 / k,垂线方程为:y = (-1/k) * (x - point.x) + point.y 。
  联立两直线方程解得:x = ( k^2 * pt1.x + k * (point.y - pt1.y ) + point.x ) / ( k^2 + 1) ,y = k * ( x - pt1.x) + pt1.y;然后再判断垂足是否在线段上,如果在线段上则返回垂足;如果不在则计算两端点到垂足的距离,选择距离垂足较近的端点返回。
22.  计算点到折线、矩形、多边形的最近点:
  只要分别计算点到每条线段的最近点,记录最近距离,取其中最近距离最小的点即可。
23.  计算点到圆的最近距离及交点坐标:
  如果该点在圆心,因为圆心到圆周任一点的距离相等,返回UNDEFINED。
  连接点P和圆心O,如果PO平行于X轴,则根据P在O的左边还是右边计算出最近点的横坐标为centerPoint.x - radius 或 centerPoint.x + radius。如果PO平行于Y轴,则根据P在O的上边还是下边计算出最近点的纵坐标为 centerPoint.y -+radius或 centerPoint.y - radius。如果PO不平行于X轴和Y轴,则PO的斜率存在且不为0,这时直线PO斜率为k = ( P.y - O.y )/ ( P.x - O.x )。直线PO的方程为:y = k * ( x - P.x) + P.y。设圆方程为(x - O.x ) ^2 + ( y - O.y ) ^2 = r ^2,联立两方程组可以解出直线PO和圆的交点,取其中离P点较近的交点即可。
24.  计算两条共线的线段的交点:
  对于两条共线的线段,它们之间的位置关系有下图所示的几种情况。图(a)中两条线段没有交点;图 (b) 和 (d) 中两条线段有无穷焦点;图 (c) 中两条线段有一个交点。设line1是两条线段中较长的一条,line2是较短的一条,如果line1包含了line2的两个端点,则是图(d)的情况,两线段有无穷交点;如果line1只包含line2的一个端点,那么如果line1的某个端点等于被line1包含的line2的那个端点,则是图(c)的情况,这时两线段只有一个交点,否则就是图(b)的情况,两线段也是有无穷的交点;如果line1不包含line2的任何端点,则是图(a)的情况,这时两线段没有交点。
25.  计算线段或直线与线段的交点:
  设一条线段为L0 = P1P2,另一条线段或直线为L1 = Q1Q2 ,要计算的就是L0和L1的交点。
 1. 首先判断L0和L1是否相交(方法已在前文讨论过),如果不相交则没有交点,否则说明L0和L1一定有交点,下面就将L0和L1都看作直线来考虑。
 2. 如果P1和P2横坐标相同,即L0平行于Y轴
  a) 若L1也平行于Y轴,
    i. 若P1的纵坐标和Q1的纵坐标相同,说明L0和L1共线,假如L1是直线的话他们有无穷的交点,假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点(该方法在前文已讨论过);
    ii. 否则说明L0和L1平行,他们没有交点;
  b) 若L1不平行于Y轴,则交点横坐标为P1的横坐标,代入到L1的直线方程中可以计算出交点纵坐标;
 3. 如果P1和P2横坐标不同,但是Q1和Q2横坐标相同,即L1平行于Y轴,则交点横坐标为Q1的横坐标,代入到L0的直线方程中可以计算出交点纵坐标;
 4. 如果P1和P2纵坐标相同,即L0平行于X轴
  a) 若L1也平行于X轴,
    i. 若P1的横坐标和Q1的横坐标相同,说明L0和L1共线,假如L1是直线的话他们有无穷的交点,假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点(该方法在前文已讨论过);
    ii. 否则说明L0和L1平行,他们没有交点;
  b) 若L1不平行于X轴,则交点纵坐标为P1的纵坐标,代入到L1的直线方程中可以计算出交点横坐标;
 5. 如果P1和P2纵坐标不同,但是Q1和Q2纵坐标相同,即L1平行于X轴,则交点纵坐标为Q1的纵坐标,代入到L0的直线方程中可以计算出交点横坐标;
 6. 剩下的情况就是L1和L0的斜率均存在且不为0的情况
  a) 计算出L0的斜率K0,L1的斜率K1 ;
  b) 如果K1 = K2
    i. 如果Q1在L0上,则说明L0和L1共线,假如L1是直线的话有无穷交点,假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点(该方法在前文已讨论过);
    ii. 如果Q1不在L0上,则说明L0和L1平行,他们没有交点。
  c) 联立两直线的方程组可以解出交点来
  这个算法并不复杂,但是要分情况讨论清楚,尤其是当两条线段共线的情况需要单独考虑,所以在前文将求两条共线线段的算法单独写出来。另外,一开始就先利用矢量叉乘判断线段与线段(或直线)是否相交,如果结果是相交,那么在后面就可以将线段全部看作直线来考虑。需要注意的是,我们可以将直线或线段方程改写为ax+by+c=0的形式,这样一来上述过程的部分步骤可以合并,缩短了代码长度,但是由于先要求出参数,这种算法将花费更多的时间。
26.  求线段或直线与折线、矩形、多边形的交点:
  分别求与每条边的交点即可。
27.  求线段或直线与圆的交点:
  设圆心为O,圆半径为r,直线(或线段)L上的两点为P1,P2。
  1. 如果L是线段且P1,P2都包含在圆O内,则没有交点;否则进行下一步。
  2. 如果L平行于Y轴,
   a) 计算圆心到L的距离dis;
   b) 如果dis > r 则L和圆没有交点;
   c) 利用勾股定理,可以求出两交点坐标,但要注意考虑L和圆的相切情况。
  3. 如果L平行于X轴,做法与L平行于Y轴的情况类似;
  4. 如果L既不平行X轴也不平行Y轴,可以求出L的斜率K,然后列出L的点斜式方程,和圆方程联立即可求解出L和圆的两个交点;
  5. 如果L是线段,对于2,3,4中求出的交点还要分别判断是否属于该线段的范围内。
28.  凸包的概念:
  点集Q的凸包(convex hull)是指一个最小凸多边形,满足Q中的点或者在多边形边上或者在其内。下图中由红色线段表示的多边形就是点集Q={p0,p1,...p12}的凸包。
29.  凸包的求法:
  现在已经证明了凸包算法的时间复杂度下界是O(n*logn),但是当凸包的顶点数h也被考虑进去的话,Krikpatrick和Seidel的剪枝搜索算法可以达到O(n*logh),在渐进意义下达到最优。最常用的凸包算法是Graham扫描法和Jarvis步进法。本文只简单介绍一下Graham扫描法,其正确性的证明和Jarvis步进法的过程大家可以参考《算法导论》。
  对于一个有三个或以上点的点集Q,Graham扫描法的过程如下:
  令p0为Q中Y-X坐标排序下最小的点
  设<p1,p2,...pm>为对其余点按以p0为中心的极角逆时针排序所得的点集(如果有多个点有相同的极角,除了距p0最远的点外全部移除
  压p0进栈S
  压p1进栈S
  压p2进栈S
for i ← 3 to m
do while 由S的栈顶元素的下一个元素、S的栈顶元素以及pi构成的折线段不拐向左侧
对S弹栈
压pi进栈S
return S;
  此过程执行后,栈S由底至顶的元素就是Q的凸包顶点按逆时针排列的点序列。需要注意的是,我们对点按极角逆时针排序时,并不需要真正求出极角,只需要求出任意两点的次序就可以了。而这个步骤可以用前述的矢量叉积性质实现。
四、结语
  尽管人类对几何学的研究从古代起便没有中断过,但是具体到借助计算机来解决几何问题的研究,还只是停留在一个初级阶段,无论从应用领域还是发展前景来看,计算几何学都值得我们认真学习、加以运用,希望这篇文章能带你走进这个丰富多彩的世界。
目录
㈠ 点的基本运算
1. 平面上两点之间距离 1
2. 判断两点是否重合 1
3. 矢量叉乘 1
4. 矢量点乘 2
5. 判断点是否在线段上 2
6. 求一点饶某点旋转后的坐标 2
7. 求矢量夹角 2

㈡ 线段及直线的基本运算
1. 点与线段的关系 3
2. 求点到线段所在直线垂线的垂足 4
3. 点到线段的最近点 4
4. 点到线段所在直线的距离 4
5. 点到折线集的最近距离 4
6. 判断圆是否在多边形内 5
7. 求矢量夹角余弦 5
8. 求线段之间的夹角 5
9. 判断线段是否相交 6
10.判断线段是否相交但不交在端点处 6
11.求线段所在直线的方程 6
12.求直线的斜率 7
13.求直线的倾斜角 7
14.求点关于某直线的对称点 7
15.判断两条直线是否相交及求直线交点 7
16.判断线段是否相交,如果相交返回交点 7

㈢ 多边形常用算法模块
1. 判断多边形是否简单多边形 8
2. 检查多边形顶点的凸凹性 9
3. 判断多边形是否凸多边形 9
4. 求多边形面积 9
5. 判断多边形顶点的排列方向,方法一 10
6. 判断多边形顶点的排列方向,方法二 10
7. 射线法判断点是否在多边形内 10
8. 判断点是否在凸多边形内 11
9. 寻找点集的graham算法 12
10.寻找点集凸包的卷包裹法 13
11.判断线段是否在多边形内 14
12.求简单多边形的重心 15
13.求凸多边形的重心 17
14.求肯定在给定多边形内的一个点 17
15.求从多边形外一点出发到该多边形的切线 18
16.判断多边形的核是否存在 19

㈣ 圆的基本运算
1 .点是否在圆内 20
2 .求不共线的三点所确定的圆 21

㈤ 矩形的基本运算
1.已知矩形三点坐标,求第4点坐标 22

㈥ 常用算法的描述 22

㈦ 补充
1.两圆关系: 24
2.判断圆是否在矩形内: 24
3.点到平面的距离: 25
4.点是否在直线同侧: 25
5.镜面反射线: 25
6.矩形包含: 26
7.两圆交点: 27
8.两圆公共面积: 28
9. 圆和直线关系: 29
10. 内切圆: 30
11. 求切点: 31
12. 线段的左右旋: 31
13.公式: 32





/* 需要包含的头文件 */
#include <cmath >

/* 常用的常量定义 */
const double INF = 1E200
const double EP = 1E-10
const int MAXV = 300
const double PI = 3.14159265

/* 基本几何结构 */
struct POINT
{
double x;
double y; POINT(double a=0, double b=0) { x=a; y=b;} //constructor
};
struct LINESEG
{
POINT s;
POINT e; LINESEG(POINT a, POINT b) { s=a; e=b;}
LINESEG() { }
};
struct LINE // 直线的解析方程 a*x+b*y+c=0 为统一表示,约定 a >= 0
{
double a;
double b;
double c; LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;}
};

/********************\
* *
* 点的基本运算 *
* *
\********************/

double dist(POINT p1,POINT p2) // 返回两点之间欧氏距离
{
return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) );
}
bool equal_point(POINT p1,POINT p2) // 判断两个点是否重合
{
return ( (abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP) );
}

/******************************************************************************
r=multiply(sp,ep,op),得到(sp-op)*(ep-op)的叉积
r>0:ep在矢量opsp的逆时针方向;
r=0:opspep三点共线;
r<0:ep在矢量opsp的顺时针方向
*******************************************************************************/

double multiply(POINT sp,POINT ep,POINT op)
{
return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y));
}

/*******************************************************************************
r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积,如果两个矢量都非零矢量
r<0:两矢量夹角为锐角;r=0:两矢量夹角为直角;r>0:两矢量夹角为钝角
*******************************************************************************/
double dotmultiply(POINT p1,POINT p2,POINT p0)
{
return ((p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y));
}

/* 判断点p是否在线段l上,条件:(p在线段l所在的直线上)&& (点p在以线段l为对角线的矩形内) */
bool online(LINESEG l,POINT p)
{
return((multiply(l.e,p,l.s)==0)
&&( ( (p.x-l.s.x)*(p.x-l.e.x)<=0 )&&( (p.y-l.s.y)*(p.y-l.e.y)<=0 ) ) );
}

// 返回点p以点o为圆心逆时针旋转alpha(单位:弧度)后所在的位置
POINT rotate(POINT o,double alpha,POINT p)
{
POINT tp;
p.x-=o.x;
p.y-=o.y;
tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x;
tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y;
return tp;
}

/* 返回顶角在o点,起始边为os,终止边为oe的夹角(单位:弧度)
角度小于pi,返回正值
角度大于pi,返回负值
可以用于求线段之间的夹角
*/
double angle(POINT o,POINT s,POINT e)
{
double cosfi,fi,norm;
double dsx = s.x - o.x;
double dsy = s.y - o.y;
double dex = e.x - o.x;
double dey = e.y - o.y;

cosfi=dsx*dex+dsy*dey;
norm=(dsx*dsx+dey*dey)*(dex*dex+dey*dey);
cosfi /= sqrt( norm );

if (cosfi >= 1.0 ) return 0;
if (cosfi <= -1.0 ) return -3.1415926;

fi=acos(cosfi);
if (dsx*dey-dsy*dex>0) return fi; // 说明矢量os 在矢量 oe的顺时针方向
return -fi;
}





/*****************************\
* *
* 线段及直线的基本运算 *
* *
\*****************************/

/* 判断点与线段的关系,用途很广泛
本函数是根据下面的公式写的,P是点C到线段AB所在直线的垂足

AC dot AB
r = ---------
||AB||^2
(Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
= -------------------------------
L^2

r has the following meaning:

r=0 P = A
r=1 P = B
r<0 P is on the backward extension of AB
r>1 P is on the forward extension of AB
0<r<1 P is interior to AB
*/
double relation(POINT p,LINESEG l)
{
LINESEG tl;
tl.s=l.s;
tl.e=p;
return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e));
}

// 求点C到线段AB所在直线的垂足 P
POINT perpendicular(POINT p,LINESEG l)
{
double r=relation(p,l);
POINT tp;
tp.x=l.s.x+r*(l.e.x-l.s.x);
tp.y=l.s.y+r*(l.e.y-l.s.y);
return tp;
}
/* 求点p到线段l的最短距离,并返回线段上距该点最近的点np
注意:np是线段l上到点p最近的点,不一定是垂足 */
double ptolinesegdist(POINT p,LINESEG l,POINT &np)
{
double r=relation(p,l);
if(r<0)
{
np=l.s;
return dist(p,l.s);
}
if(r>1)
{
np=l.e;
return dist(p,l.e);
}
np=perpendicular(p,l);
return dist(p,np);
}

// 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别
double ptoldist(POINT p,LINESEG l)
{
return abs(multiply(p,l.e,l.s))/dist(l.s,l.e);
}

/* 计算点到折线集的最近距离,并返回最近点.
注意:调用的是ptolineseg()函数 */
double ptopointset(int vcount,POINT pointset,POINT p,POINT &q)
{
int i;
double cd=double(INF),td;
LINESEG l;
POINT tq,cq;

for(i=0;i<vcount-1;i++)
{
l.s=pointseti ];
l.e=pointset;
td=ptolinesegdist(p,l,tq);
if(td<cd)
{
cd=td;
cq=tq;
}
}
q=cq;
return cd;
}
/* 判断圆是否在多边形内.ptolineseg()函数的应用2 */
bool CircleInsidePolygon(int vcount,POINT center,double radius,POINT polygon)
{
POINT q;
double d;
q.x=0;
q.y=0;
d=ptopointset(vcount,polygon,center,q);
if(d<radius||fabs(d-radius)<EP)
return true;
else
return false;
}

/* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意:如果想从余弦求夹角的话,注意反余弦函数的定义域是从 0到pi */
double cosine(LINESEG l1,LINESEG l2)
{
return (((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x) +
(l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) );
}
// 返回线段l1与l2之间的夹角 单位:弧度 范围(-pi,pi)
double lsangle(LINESEG l1,LINESEG l2)
{
POINT o,s,e;
o.x=o.y=0;
s.x=l1.e.x-l1.s.x;
s.y=l1.e.y-l1.s.y;
e.x=l2.e.x-l2.s.x;
e.y=l2.e.y-l2.s.y;
return angle(o,s,e);
}
// 如果线段u和v相交(包括相交在端点处)时,返回true
bool intersect(LINESEG u,LINESEG v)
{
return( (max(u.s.x,u.e.x)>=min(v.s.x,v.e.x))&& //排斥实验
(max(v.s.x,v.e.x)>=min(u.s.x,u.e.x))&&
(max(u.s.y,u.e.y)>=min(v.s.y,v.e.y))&&
(max(v.s.y,v.e.y)>=min(u.s.y,u.e.y))&&
(multiply(v.s,u.e,u.s)*multiply(u.e,v.e,u.s)>=0)&& //跨立实验
(multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0));
}


// (线段u和v相交)&&(交点不是双方的端点) 时返回true
bool intersect_A(LINESEG u,LINESEG v)
{
return((intersect(u,v))&&
(!online(u,v.s))&&
(!online(u,v.e))&&
(!online(v,u.e))&&
(!online(v,u.s)));
}


// 线段v所在直线与线段u相交时返回true;方法:判断线段u是否跨立线段v
bool intersect_l(LINESEG u,LINESEG v)
{
return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0;
}


// 根据已知两点坐标,求过这两点的直线解析方程: a*x+b*y+c = 0 (a >= 0)
LINE makeline(POINT p1,POINT p2)
{
LINE tl;
int sign = 1;
tl.a=p2.y-p1.y;
if(tl.a<0)
{
sign = -1;
tl.a=sign*tl.a;
}
tl.b=sign*(p1.x-p2.x);
tl.c=sign*(p1.y*p2.x-p1.x*p2.y);
return tl;
}

// 根据直线解析方程返回直线的斜率k,水平线返回 0,竖直线返回 1e200
double slope(LINE l)
{
if(abs(l.a) < 1e-20)return 0;
if(abs(l.b) < 1e-20)return INF;
return -(l.a/l.b);
}

// 返回直线的倾斜角alpha ( 0 - pi)
double alpha(LINE l)
{
if(abs(l.a)< EP)return 0;
if(abs(l.b)< EP)return PI/2;
double k=slope(l);
if(k>0)
return atan(k);
else
return PI+atan(k);
}

// 求点p关于直线l的对称点
POINT symmetry(LINE l,POINT p)
{
POINT tp;
tp.x=((l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b);
tp.y=((l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b);
return tp;
}

// 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交,返回true,且返回交点p
bool lineintersect(LINE l1,LINE l2,POINT &p) // 是 L1,L2
{
double d=l1.a*l2.b-l2.a*l1.b;
if(abs(d)<EP) // 不相交
return false;
p.x = (l2.c*l1.b-l1.c*l2.b)/d;
p.y = (l2.a*l1.c-l1.a*l2.c)/d;
return true;
}

// 如果线段l1和l2相交,返回true且交点由(inter)返回,否则返回false
bool intersection(LINESEG l1,LINESEG l2,POINT &inter)
{
LINE ll1,ll2;
ll1=makeline(l1.s,l1.e);
ll2=makeline(l2.s,l2.e);
if(lineintersect(ll1,ll2,inter))
{
return online(l1,inter);
}
else
return false;
}

本库题目(必须完成)
479 计算几何相关基本知识
480 多边形面积和周长
481 点在多边形上的判定
207  凸包测试
659 求两个凸多边形的交集面积
660 求凸多边形的面积
661 求平面凸包
662 屠龙传说-屠龙枪卷
1391 nh 
1392 国王的一周
334 Meksikanac
402 折纸
1946  Pick
1947 Svjetlost

可以应付大部分zoj上的几何题 Sample TextSample Text

zoj上的计算几何题
Vol I
1010 by pandahyx
1032 by javaman
1037 by Vegetable Bird
1041 by javaman
1081 by Vegetable Bird
1090 by Vegetable Bird

Vol II
1104 by javaman
1123 by javaman
1139 by Vegetable Bird
1165 by javaman
1199 by Vegetable Bird

Vol V
1426 by Vegetable Bird
1439 by Vegetable Bird
1460 by Vegetable Bird
1472 by Vegetable Bird

Vol VI
1597 by Vegetable Bird

Vol VII
1608 by Vegetable Bird
1648 by Vegetable Bird

Vol XII
2102 by pandahyx
2107 by pandahyx
2157 by pandahyx

Vol XIII
2234 by pandahyx

Vol XIV
2318 by ahyangyi
2394 by qherlyt

Vol XV
2403 by Vegetable Bird

网络流专题

2022-07-13 15:11:20 By root

网络流算法总结
http://www.cnblogs.com/longdouhzt/archive/2012/05/20/2510753.html
https://www.cnblogs.com/yangsongyi/p/10167165.html
总体上来说,最大流算法分为两大类:增广路 (Augmenting Path) 和预流推进重标号 (Push Relabel) 。也有算法同时借鉴了两者的长处,如 Improved SAP 。本篇主要介绍增广路类算法,思想、复杂度及实际运行效率比较,并试图从中选择一种兼顾代码复杂度和运行效率的较好方案。以下我们将会看到,有时理论分析的时间复杂度并不能很好的反映一种算法的实际效率。

  1. Ford - Fulkerson 方法 所有增广路算法的基础都是 Ford - Fulkerson 方法。称之为方法而不是算法是因为 Ford - Fulkerson 只提供了一类思想,在此之上的具体操作可有不同的实现方案。

给定一个有向网络 G(V,E) 以及源点 s 终点 t ,FF 方法描述如下:

Ford-Fulkerson 方法 (G,s,t) 1 将各边上流量 f 初始化为 0 2 while 存在一条增广路径 p 3 do 沿路径 p 增广流量 f 4 return f 假设有向网络 G 中边 (i,j) 的容量为 c(i,j) ,当前流量为 f(i,j) ,则此边的剩余流量即为 r(i,j) = c(i,j) - f(i,j) ,其反向边的剩余流量为 r(j,i) = f(i,j) 。有向网中所有剩余流量 r(i,j) > 0 的边构成残量网络 Gf ,增广路径p即是残量网络中从源点 s 到终点 t 的路径。

沿路径 p 增广流量 f 的操作基本都是相同的,各算法的区别就在于寻找增广路径 p 的方法不同。例如可以寻找从 s 到 t 的最短路径,或者流量最大的路径。

  1. Edmonds - Karp 算法 Shortest Augmenting Path (SAP) 是每次寻找最短增广路的一类算法,Edmonds - Karp 算法以及后来著名的 Dinic 算法都属于此。SAP 类算法可统一描述如下:

复制代码 Shortest Augmenting Path 1 x <-- 0 2 while 在残量网络 Gx 中存在增广路 s ~> t 3 do 找一条最短的增广路径 P 4 delta <-- min{rij:(i,j) 属于 P} 5 沿 P 增广 delta 大小的流量 6 更新残量网络 Gx 7 return x 复制代码 在无权边的有向图中寻找最短路,最简单的方法就是广度优先搜索 (BFS),E-K 算法就直接来源于此。每次用一遍 BFS 寻找从源点 s 到终点 t 的最短路作为增广路径,然后增广流量 f 并修改残量网络,直到不存在新的增广路径。

E-K 算法的时间复杂度为 O(VE2),由于 BFS 要搜索全部小于最短距离的分支路径之后才能找到终点,因此可以想象频繁的 BFS 效率是比较低的。实践中此算法使用的机会较少。

  1. Dinic 算法 BFS 寻找终点太慢,而 DFS 又不能保证找到最短路径。1970年 Dinic 提出一种思想,结合了 BFS 与 DFS 的优势,采用构造分层网络的方法可以较快找到最短增广路,此算法又称为阻塞流算法 (Blocking Flow Algorithm)。

首先定义分层网络 AN(f)。在残量网络中从源点 s 起始进行 BFS,这样每个顶点在 BFS 树中会得到一个距源点 s 的距离 d,如 d(s) = 0,直接从 s 出发可到达的点距离为 1,下一层距离为2 ... 。称所有具有相同距离的顶点位于同一层,在分层网络中,只保留满足条件 d(i) + 1 = d(j) 的边,这样在分层网络中的任意路径就成为到达此顶点的最短路径。

Dinic 算法每次用一遍 BFS 构建分层网络 AN(f),然后在 AN(f) 中一遍 DFS 找到所有到终点 t 的路径增广;之后重新构造 AN(f),若终点 t 不在 AN(f) 中则算法结束。DFS 部分算法可描述如下:

复制代码 1 p <-- s 2 while s 的出度 > 0 do 3 u <-- p.top 4 if u != t then 5 if u 的出度 > 0 then 6 设 (u,v) 为 AN(f) 中一条边 7 p <-- p, v 8 else 9 从 p 和 AN(f) 中删除点 u 以及和 u 连接的所有边 10 else 11 沿 p 增广 12 令 p.top 为从 s 沿 p 可到达的最后顶点 13 end while 复制代码

实际代码中不必真的用一个图来存储分层网络,只需保存每个顶点的距离标号并在 DFS 时判断 dist[i] + 1 = dist[j] 即可。Dinic 的时间复杂度为 O(V2E)。由于较少的代码量和不错的运行效率,Dinic 在实践中比较常用。具体代码可参考 DD_engi 网络流算法评测包中的标程,这几天 dinic 算法的实现一共看过和比较过将近 10 个版本了,DD 写的那个在效率上是数一数二的,逻辑上也比较清晰。

  1. Improved SAP 算法 本次介绍的重头戏。通常的 SAP 类算法在寻找增广路时总要先进行 BFS,BFS 的最坏情况下复杂度为 O(E),这样使得普通 SAP 类算法最坏情况下时间复杂度达到了 O(VE2)。为了避免这种情况,Ahuja 和 Orlin 在1987年提出了Improved SAP 算法,它充分利用了距离标号的作用,每次发现顶点无出弧时不是像 Dinic 算法那样到最后进行 BFS,而是就地对顶点距离重标号,这样相当于在遍历的同时顺便构建了新的分层网络,每轮寻找之间不必再插入全图的 BFS 操作,极大提高了运行效率。国内一般把这个算法称为 SAP...显然这是不准确的,毕竟从字面意思上来看 E-K 和 Dinic 都属于 SAP,我还是习惯称为 ISAP 或改进的 SAP 算法。

    与 Dinic 算法不同,ISAP 中的距离标号是每个顶点到达终点 t 的距离。同样也不需显式构造分层网络,只要保存每个顶点的距离标号即可。程序开始时用一个反向 BFS 初始化所有顶点的距离标号,之后从源点开始,进行如下三种操作:(1)当前顶点 i 为终点时增广 (2) 当前顶点有满足 dist[i] = dist[j] + 1 的出弧时前进 (3) 当前顶点无满足条件的出弧时重标号并回退一步。整个循环当源点 s 的距离标号 dist[s] >= n 时结束。对 i 点的重标号操作可概括为 dist[i] = 1 + min{dist[j] : (i,j)属于残量网络Gf}。具体算法描述如下:

复制代码 algorithm Improved-Shortest-Augmenting-Path 1 f <-- 0 2 从终点 t 开始进行一遍反向 BFS 求得所有顶点的起始距离标号 d(i) 3 i <-- s 4 while d(s) < n do 5 if i = t then // 找到增广路 6 Augment 7 i <-- s // 从源点 s 开始下次寻找 8 if Gf 包含从 i 出发的一条允许弧 (i,j) 9 then Advance(i) 10 else Retreat(i) // 没有从 i 出发的允许弧则回退 11 return f

procedure Advance(i) 1 设 (i,j) 为从 i 出发的一条允许弧 2 pi(j) <-- i // 保存一条反向路径,为回退时准备 3 i <-- j // 前进一步,使 j 成为当前结点

procedure Retreat(i) 1 d(i) <-- 1 + min{d(j):(i,j)属于残量网络Gf} // 称为重标号的操作 2 if i != s 3 then i <-- pi(i) // 回退一步

procedure Augment 1 pi 中记录为当前找到的增广路 P 2 delta <-- min{rij:(i,j)属于P} 3 沿路径 P 增广 delta 的流量 4 更新残量网络 Gf 复制代码

算法中的允许弧是指在残量网络中满足 dist[i] = dist[j] + 1 的弧。Retreat 过程中若从 i 出发没有弧属于残量网络 Gf 则把顶点距离重标号为 n 。

虽然 ISAP 算法时间复杂度与 Dinic 相同都是 O(V2E),但在实际表现中要好得多。要提的一点是关于 ISAP 的一个所谓 GAP 优化。由于从 s 到 t 的一条最短路径的顶点距离标号单调递减,且相邻顶点标号差严格等于1,因此可以预见如果在当前网络中距离标号为 k (0 <= k < n) 的顶点数为 0,那么可以知道一定不存在一条从 s 到 t 的增广路径,此时可直接跳出主循环。在我的实测中,这个优化是绝对不能少的,一方面可以提高速度,另外可增强 ISAP 算法时间上的稳定性,不然某些情况下 ISAP 会出奇的费时,而且大大慢于 Dinic 算法。相对的,初始的一遍 BFS 却是可有可无,因为 ISAP 可在循环中自动建立起分层网络。实测加不加 BFS 运行时间差只有 5% 左右,代码量可节省 15~20 行。

  1. 最大容量路径算法 (Maximum Capacity Path Algorithm) 1972年还是那个 E-K 组合提出的另一种最大流算法。每次寻找增广路径时不找最短路径,而找容量最大的。可以预见,此方法与 SAP 类算法相比可更快逼近最大流,从而降低增广操作的次数。实际算法也很简单,只用把前面 E-K 算法的 BFS 部分替换为一个类 Dijkstra 算法即可。USACO 4.2 节的说明详细介绍了此算法,这里就不详述了。

    时间复杂度方面。BFS 是 O(E),简单 Dijkstra 是 O(V2),因此效果可想而知。但提到 Dijkstra 就不能不提那个 Heap 优化,虽然 USACO 的算法例子中没有用 Heap ,我自己还是实现了一个加 Heap 的版本,毕竟 STL 的优先队列太好用了不加白不加啊。效果也是非常明显的,但比起 Dinic 或 ISAP 仍然存在海量差距,这里就不再详细介绍了。

  2. Capacity Scaling Algorithm 不知道怎么翻比较好,索性就这么放着吧。叫什么的都有,容量缩放算法、容量变尺度算法等,反正就那个意思。类似于二分查找的思想,寻找增广路时不必非要局限于寻找最大容量,而是找到一个可接受的较大值即可,一方面有效降低寻找增广路时的复杂度,另一方面增广操作次数也不会增加太多。时间复杂度 O(E2logU) 实际效率嘛大约稍好于最前面 BFS 的 E-K 算法,稀疏图时表现较优,但仍然不敌 Dinic 与 ISAP。

  3. 算法效率实测! 重头戏之二,虽然引用比较多,哎~

    首先引用此篇强文 《Maximum Flow: Augmenting Path Algorithms Comparison》

    对以上算法在稀疏图、中等稠密图及稠密图上分别进行了对比测试。直接看结果吧:

稀疏图:

ISAP 轻松拿下第一的位置,图中最左边的 SAP 应该指的是 E-K 算法,这里没有比较 Dinic 算法是个小遗憾吧,他把 Dinic 与 SAP 归为一类了。最大流量路径的简单 Dijkstra 实现实在是太失败了 - -,好在 Heap 优化后还比较能接受……可以看到 Scaling 算法也有不错的表现。

中等稠密图:

ISAP 依然领先。最大流量算法依然不太好过……几个 Scaling 类算法仍然可接受。

稠密图:

ISAP……你无敌了!这次可以看出 BFS 的 Scaling 比 DFS 实现好得多,而且几乎与 Improved Scaling 不相上下,代码量却差不少。看来除 ISAP 外 BFS Scaling 也是个不错的选择。

最后补个我自己实测的图,比较算法有很多是 DD 网络流算法评测包里的标程,评测系统用的 Cena,评测数据为 DD ditch 数据生成程序生成的加强版数据:

我一共写了 7 个版本的 ISAP ,Dinic 算法也写了几个递归版的但效率都不高,只放上来一个算了。从上图来看似乎 ISAP 全面超越了号称最大流最快速算法的 HLPP,但在另外一台机器上测试结果与此却不大相同,有时 ISAP 与 HLPP 基本持平,有时又稍稍慢一些。在这种差距非常小的情况下似乎编译器的效果也比较明显。那个 HLPP 是用 PASCAL 写的,我不知在 Win32 平台下编译代码效率如何,至少我的几个 ISAP 用 VC2008 + SP1 编译比用 g++ 要强不少,也可能是参数设置的问题。

不过这些都是小事,关键是见证了 ISAP 的实际效率。从上面可以看出不加 GAP 优化的 ISAP 有几个测试点干脆无法通过,而不加 BFS 却无甚大影响,递归与非递归相差在 7% 左右的样子。综合以上表现,推荐采用 ISAP 不加 BFS,非递归 + GAP 优化的写法,Ditch 这道题一共也才 80 行左右代码。要提的一点是 GAP 优化用递归来表现的话不如 while 循环来得直接。另外,看起来 HLPP 表现确实很优秀,有机会也好好研究一下吧,预流推进重标号算法也是一大类呢……

最后附上一个 ISAP + GAP + BFS 的非递归版本代码,网络采用邻接表 + 反向弧指针:

#include<cstdio>
#include<memory>
using namespace std;
const int maxnode = 1024;
const int infinity = 2100000000;
struct edge{
    int ver;    // vertex
    int cap;    // capacity
    int flow;   // current flow in this arc
    edge *next; // next arc
    edge *rev;  // reverse arc
    edge(){}
    edge(int Vertex, int Capacity, edge *Next)
        :ver(Vertex), cap(Capacity), flow(0), next(Next), rev((edge*)NULL){}
    void* operator new(size_t, void *p){
        return p;
    }
}*Net[maxnode];
int dist[maxnode]= {0}, numbs[maxnode] = {0}, src, des, n;
void rev_BFS(){
    int Q[maxnode], head = 0, tail = 0;
    for(int i=1; i<=n; ++i){
        dist[i] = maxnode;
        numbs[i] = 0;
    }
    Q[tail++] = des;
    dist[des] = 0;
    numbs[0] = 1;
    while(head != tail){
        int v = Q[head++];
        for(edge *e = Net[v]; e; e = e->next){
            if(e->rev->cap == 0 || dist[e->ver] < maxnode )continue;
            dist[e->ver] = dist[v] + 1;
            ++numbs[dist[e->ver]];
            Q[tail++] = e->ver;
        }
    }
}
int maxflow(){
    int u, totalflow = 0;
    edge *CurEdge[maxnode], *revpath[maxnode];
    for(int i=1; i < = n; ++i)CurEdge[i] = Net[i];
    u = src;
    while(dist[src] < n){
        if(u == des){    // find an augmenting path
            int augflow = infinity;
            for(int i = src; i != des; i = CurEdge[i]->ver)
                augflow = min(augflow, CurEdge[i]->cap);
            for(int i = src; i != des; i = CurEdge[i]->ver){
                CurEdge[i]->cap -= augflow;
                CurEdge[i]->rev->cap += augflow;
                CurEdge[i]->flow += augflow;
                CurEdge[i]->rev->flow -= augflow;
            }
            totalflow += augflow;
            u = src;
        }
        edge *e;
        for(e = CurEdge[u]; e; e = e->next)
            if(e->cap > 0 && dist[u] == dist[e->ver] + 1)break;
        if(e){    // find an admissible arc, then Advance
            CurEdge[u] = e;
            revpath[e->ver] = e->rev;
            u = e->ver;
        } else {    // no admissible arc, then relabel this vertex
            if(0 == (--numbs[dist[u]]))break;    // GAP cut, Important!
            CurEdge[u] = Net[u];
            int mindist = n;
            for(edge *te = Net[u]; te; te = te->next)
                if(te->cap > 0)mindist = min(mindist, dist[te->ver]);
            dist[u] = mindist + 1;
            ++numbs[dist[u]];
            if(u != src)
                u = revpath[u]->ver;    // Backtrack
        }
    }
    return totalflow;
}
int main(){
    int m, u, v, w;
    freopen("ditch.in", "r", stdin);
    freopen("ditch.out", "w", stdout);
    while(scanf("%d%d", &m, &n) != EOF){    // POJ 1273 need this while loop
        edge *buffer = new edge[2*m];
        edge *data = buffer;
        src = 1; des = n;
        while(m--){
            scanf("%d%d%d", &u, &v, &w);
            Net[u] = new((void*) data++) edge(v, w, Net[u]);
            Net[v] = new((void*) data++) edge(u, 0, Net[v]);
            Net[u]->rev = Net[v];
            Net[v]->rev = Net[u];
        }
        rev_BFS();
        printf("%d\n", maxflow());
        delete [] buffer;
    }
    return 0;
}

网络流建模汇总
最小割模型在信息学竞赛中的应用

网络流题单
最大流:
20- K-联赛
458 - 最大流(不解释)
63 - 乒乓混双配对(匹配可做)
264 - 小行星(匹配可做)
267 - 机器生产调度(匹配可做)
587 - 最大流加强(最裸)
618 - SLICICE(推荐)
622 - 跳舞(貌似可以贪心?)
663 - 求容量有上下界的最小流(拓展)
716 - 狼和羊的故事(算经典问题了)
826 - 耙耳朵系列之家务活
397 - 蜥蜴
842 - 最大流输路径(有讨论意义,lwj表示自身没有经历过这个问题的讨论)
844 - 逃跑
855 - 猪
1322 - 飞行员配对方案问题
1323- 最小路径覆盖问题(最大流拓展)
1324 - 魔术球问题(据说可以贪心?)
1325 - 骑士共存问题(匹配问题,匈牙利复杂度不能接受)
1560 - 士兵占领

费用流题目
426 - 志愿者招募
455 - 修车
459 - 最小费用流(不解释) 已做
1214 - 科比的比赛(较难)
747 - 探险计划(综合,难)
460 - 任务安排(easy) 已做

最小割、最大权闭合图题目
418 - 最大获利 已做
702 - 猫狗大战 已做
848 - 植物大战僵尸
1266 - 海拔(对偶图转化,难)
1329 - 最小生成树

root Avatar