题解 P4502 【[ZJOI2018]保镖】

xzyxzy

2019-02-08 11:18:21

Solution

## 追求更佳阅读体验:[BLOG](https://www.cnblogs.com/xzyxzy/p/10349399.html)。记得推荐&点赞哦! # [ZJOI2018]保镖 Tags:题解 --- ## 题意 [链接](https://www.luogu.org/problemnew/show/P4502) 初始在平面上有一些点,九条可怜随机出现在一个矩形内的任意一点。若九条可怜出现在$O$点,则平面上所有的点都从$P_i$移动到$P'_i$,使得$P'_i$在射线$OP_i$上,且满足$|OP_i|*|OP'_i|=1$。现在给定矩形范围,求这些点移动后所构成的凸包的期望点数。 $n\le 2000,x,y\le 10^5$,精度要求绝对误差或相对误差不超过$10^{-7}$。 ## 题解 ### **前言** 神仙不可做题终于被杠下来了!撒花! 不得不说九老师这个多合一是出的真的牛逼!(~~比lalaxu不知道高明到哪里去了~~) 首先感谢[Ez3real](https://loj.ac/submission/101720)的代码框架(不过LOJ两人AC代码一样什么鬼)和[yuhaoxiang](https://yhx-12243.github.io/OI-transit/records/lg4502%3Bloj2530%3Bsoj111.html)的题解(这个网站很慢)。 ### **Part 1 前置知识:圆与矩形的面积交** 在[计算几何基础](https://www.cnblogs.com/xzyxzy/p/10033130.html)里有。 ### **Part 2 前置知识:三维凸包** 在[三维凸包](https://www.cnblogs.com/xzyxzy/p/10225804.html)里有。 ### **Part 3 前置知识:欧拉公式** 在[Pick定理、欧拉公式和圆的反演](https://www.cnblogs.com/xzyxzy/p/10241872.html)里有。 ### **Part 4 前置知识:反演** 这道题显然是要求平面上的点关于$O$的、以$1$为反演幂(反演半径)的反形的凸包期望点数。 至于反演是什么可以看这个:[Pick定理、欧拉公式和圆的反演](https://www.cnblogs.com/xzyxzy/p/10241872.html) ### **Part 5 前置知识:Voronoi图** 又称**泰森多边形**。 大概就是一个平面划分,平面上的每个点划分到离它最近的关键点上。 Wiki有张十分形象的图 ![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_%E6%B3%B0%E6%A3%AE%E5%A4%9A%E8%BE%B9%E5%BD%A2.gif) ### **Part 6 前置知识:Delaunay三角剖分** #### **三角剖分** 感性理解一下就是 ![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_%E4%B8%89%E8%A7%92%E5%89%96%E5%88%86.jpg) Delaunay三角剖分是一种有着优秀性质的三角剖分。 **定理:对于任何一种三角剖分,三角形个数和外围凸包点数之和为2n-2**。 这里凸包是严格凸的,也就是没有三点共线情况。 考虑用欧拉公式证明:设凸包上的点数为$k$,三角形个数为$F-1$,则有 $$V-E+F=n-((F-1)*3+k)/2+F=2$$凸包上的边算了一次,三角形上的边算了两次、 即$$k+F=2n-1,k+F-1=2n-2$$ <br> #### **立体Delaunay三角剖分** 当然我们目前只考虑平面的Delaunay三角剖分,至于立体的可以看看这张图,本文不会涉及。 ![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_Delaunay%E4%B8%89%E8%A7%92%E5%89%96%E5%88%86.png) (图片来源于网络) <br> #### **Delaunay三角剖分和泰森多边形** Delaunay三角剖分和泰森多边形是对偶图。 对偶图是什么呢,看下面的构造方法: 泰森多边形的交点一般属于三个区域,将这三个区域的标志点连起来,就得到了一个原图的三角剖分。 ![](https://img-blog.csdn.net/20141126144248328?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvbWFrZW5vdGhpbmc=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center) (图片来源于网络) 上图中,实线是泰森多边形,虚线链接,得到标志点的一个Delaunay三角剖分。 <br> #### **Delaunay三角剖分的性质** 由于其美妙的构造,可以得到一些美妙的性质: - 平面上的点集有且仅有唯一的Delaunay三角剖分(除出现四点共圆的情况,这时泰森多边形有顶点属于四个区域)。 - 任意一个Delaunay三角形的外接圆不包含点集中的其他点。(称为Delaunay三角形的空圆性质)。 - Delaunay三角剖分相比其他的三角剖分,所有三角形的最小角最大 <br> #### **Delaunay三角剖分的构造** 以下内容摘自百度百科 **Bowyer-Watson算法** 1. 构造超级三角形(类似半平面交中的超级平面)。 2. 插入点$P$,找到点$P$的影响三角形(外接圆包括点$P$的三角形),删除影响三角形的公共边,并把$P$向这些影响到的点连边。 3. 对原图进行优化 4. 重复$2$直到所有点插入完毕 步骤2图示: ![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_Delaunay%E4%B8%89%E8%A7%92%E5%89%96%E5%88%86.jpg) 步骤3图示:转变连接对角线的方式使其满足空圆性质 ![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_Delaunay%E4%B8%89%E8%A7%92%E5%89%96%E5%88%86%E5%B1%80%E9%83%A8%E4%BC%98%E5%8C%96.jpg) 和求三维凸包类似的复杂度分析,复杂度大概是$\cal O(n^2)$ ### **Part 7 初步转化** 凸包点数,只能整体地去求,由于三角剖分的定理,我们可以转而求三角剖分的三角形个数。 三角形个数是可以用期望算的。 每一个Delaunay三角形对应一个外接圆,我们称为Delaunay圆。所以题目又转化为算Delaunay圆的期望个数。若Delaunay圆的期望个数为$num$,答案就是$Ans=2n-2-num$。 定义支配圆为包含点集中所有点的圆。Delaunay圆内不包含除Delaunay三角形三个顶点外的其他任何点,所以支配圆与Delaunay圆恰好相反。 则以下结论成立 - **对于Delaunay圆,若反演中心在圆内,其反形是支配圆;否则反形还是Delaunay圆。** - **对于支配圆,若反演中心在圆内,其反形是Delaunay圆;否则反形还是支配圆。** 这题不用考虑反演中心在圆周上的情况(概率为0) [这里有yuhaoxiang大佬做的一个Geogebra演示文件](https://www.geogebra.org/m/begkbgrn),我把两种情况截图下来是这样的: Delaunay圆,反演中心在圆内,反形是支配圆 ![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_Delaunay%E5%9C%86%E7%9A%84%E5%8F%8D%E6%BC%942.png) Delaunay圆,反演中心在远外,反形还是Delaunay圆 ![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_Delaunay%E5%9C%86%E7%9A%84%E5%8F%8D%E6%BC%941.png) 题目转化为:求**原图中Delaunay圆的个数×反演中心在圆内的概率+支配圆的个数×反演中心在圆外的概率** 若其期望为$E$,则$Ans=2n-2-$反形是Delaunay圆的概率$=2+E$ (这句话看完$Part 8$再回来看)对于这$n$个点,每个点一定至少会是一个Delaunay圆对应的其中一个顶点,所以这$n$个点每个点都会出现在凸包上,则凸包一共有$2n-4$个面,所以Delaunay圆+支配圆$=2n-4$。 ### **Part 8 进一步转化** 考虑圆的方程$$x^2+y^2+Dx+Ey+F=0$$若令$z=x^2+y^2$,可以得到$Dx+Ey+z+F=0$,**这是空间坐标系中的一个平面方程** $z=x^2+y^2$长这样:![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_%E5%9C%86%E5%8F%98%E5%B9%B3%E9%9D%A21.png) 例如圆$(x-1)^2+(y-1)^2=1$可以表示为$-2x-2y+z+1=0$,长这样:![](http://images.cnblogs.com/cnblogs_com/xzyxzy/1374475/o_%E5%9C%86%E5%8F%98%E5%B9%B3%E9%9D%A22.png) 自己做了一个动画:[网址](https://www.geogebra.org/classic/xjmqkytp)。 那么如果一个点$(x,y)$在圆上,则$(x,y,x^2+y^2)$在该圆对应的平面上 同理,在圆外或圆内,对应着在平面的一侧。具体来说,在其平面上面表示在圆外,在平面下面表示在圆内。 于是,把所有点映射到三维坐标系中,求凸包,**下凸面对应Delaunay圆,上凸面对应支配圆。** ### **Part 9 实现过程** 首先读入所有的点并进行随机微小扰动,使得不存在多点共圆以及最后求出三维凸包中不存在与$z $轴平行的凸面。 然后求解三维凸包,这里采用的是增量法。 对于每个凸面,得到对应的三个点、求出其外接圆。 - 如果其为上凸面,则其为支配圆,只有在反演中心在圆外,贡献答案 - 如果其为下凸面,则其为Delaunay圆,只有反演中心在圆内,贡献答案 那么就是算给定矩形和圆形的面积交,这个在前置知识$Part 1$里啦 ### **Part 10 总结&代码** 写了一年终于写完了! 完结撒花!! 这题考场上一定要果断丢,没有部分分。 这题出得很好,考察知识点全面。很巧妙的地方是:巧妙地把圆转化成三维空间的平面,从而把平面问题转化为三维凸包问题。巧妙地运用三角剖分,把求凸包顶点期望个数变为求圆的期望个数。 不得不说,orz jiry!!! **Code** ```cpp #include<iostream> #include<cmath> #define db __float128 #define orzjiry_2 19491001 using namespace std; const db eps=1e-10; db Rand() {return 1.0*rand()/RAND_MAX;} int sign(db x) {return x<-eps?-1:(x>eps);} struct v2 { db x,y; v2 operator + (v2 a) {return (v2){x+a.x,y+a.y};} v2 operator - (v2 a) {return (v2){x-a.x,y-a.y};} v2 operator / (db t) {return (v2){x/t,y/t};} v2 operator ^ (db t) {return (v2){x*t,y*t};} db operator * (v2 a) {return x*a.y-y*a.x;} db operator & (v2 a) {return x*a.x+y*a.y;} db dis() {return sqrt((double)(x*x+y*y));} db dis2() {return x*x+y*y;} void rot() {db t=x;x=-y;y=t;} }jir[4]; struct v3 { db x,y,z; v3 operator + (v3 a) {return (v3){x+a.x,y+a.y,z+a.z};} v3 operator - (v3 a) {return (v3){x-a.x,y-a.y,z-a.z};} v3 operator * (v3 a) {return (v3){y*a.z-z*a.y,z*a.x-x*a.z,x*a.y-y*a.x};} db operator & (v3 a) {return x*a.x+y*a.y+z*a.z;} void shake() {x+=Rand()*1e-10,y+=Rand()*1e-10,z+=Rand()*1e-10;} }P[2100]; struct Face { int v[3]; v3 Normal() {return (P[v[1]]-P[v[0]])*(P[v[2]]-P[v[0]]);} }F[8100],C[8100]; int n,cnt; namespace TAT2D { v2 Cross(v2 a1,v2 a2,v2 b1,v2 b2) { v2 a=a2-a1,b=b2-b1,c=b1-a1; db t=(b*c)/(b*a); return a1+(a^t); } int cmp(db a,db b) {return sign(a-b);} db rad(v2 p1,v2 p2) {return atan2(double(p1*p2),double(p1&p2));} db Calc(db r,v2 p1,v2 p2) { v2 e=(p1-p2)/(p1-p2).dis(),e1=e;e.rot(); v2 mid=Cross(p1,p2,(v2){0,0},e),d1=mid; if(d1.dis()>r) return r*r*rad(p1,p2)/2; db d=sqrt(double(r*r-d1.dis2())); v2 w1=mid+(e1^d),w2=mid-(e1^d); int b1=cmp(p1.dis2(),r*r)==1,b2=cmp(p2.dis2(),r*r)==1; if(b1&&b2) { if(sign((p1-w1)&(p2-w1))<=0) return r*r*(rad(p1,w1)+rad(w2,p2))/2+(w1*w2)/2; else return r*r*rad(p1,p2)/2; } if(b1) return (r*r*rad(p1,w1)+w1*p2)/2; if(b2) return (p1*w2+r*r*rad(w2,p2))/2; return p1*p2/2; } db intersect(v2 O,db r) { db res=0; for(int i=0;i<4;i++) res+=Calc(r,jir[i]-O,jir[(i+1)%4]-O); return res; } } namespace TAT3D { bool vis[2100][2100]; int see(Face a,v3 b) {return ((b-P[a.v[0]])&a.Normal())>0;} void Convex() { for(int i=0;i<n;i++) P[i].shake(); int cc=-1;cnt=-1; F[++cnt]=(Face){0,1,2}; F[++cnt]=(Face){2,1,0}; for(int i=3;i<n;i++) { for(int j=0,v;j<=cnt;j++) { if(!(v=see(F[j],P[i]))) C[++cc]=F[j]; for(int k=0;k<3;k++) vis[F[j].v[k]][F[j].v[(k+1)%3]]=v; } for(int j=0;j<=cnt;j++) for(int k=0;k<3;k++) { int x=F[j].v[k],y=F[j].v[(k+1)%3]; if(vis[x][y]&&!vis[y][x]) C[++cc]=(Face){x,y,i}; } for(int j=0;j<=cc;j++) F[j]=C[j]; cnt=cc;cc=-1; } } } int main() { //Part 1 输入以及初步转化 srand(orzjiry_2); double xx,yy; cin>>n>>xx>>yy;jir[0]=(v2){xx,yy}; cin>> xx>>yy;jir[2]=(v2){xx,yy}; jir[1]=(v2){jir[2].x,jir[0].y}; jir[3]=(v2){jir[0].x,jir[2].y}; db S=(jir[2].x-jir[0].x)*(jir[2].y-jir[0].y),Ans=0; for(int i=0;i<n;i++) { double x,y;cin>>x>>y; P[i]=(v3){x,y,x*x+y*y}; } //Part 2 计算三维凸包并求出反演后支配圆的期望数量 TAT3D::Convex(); for(int i=0;i<=cnt;i++) { v3 o=F[i].Normal(); v2 a1=(v2){P[F[i].v[0]].x,P[F[i].v[0]].y}; v2 a2=(v2){P[F[i].v[1]].x,P[F[i].v[1]].y}; v2 c=(a1+a2)/2.0,d=a2-a1;d.rot(); a1=c;a2=c+d; v2 b1=(v2){P[F[i].v[1]].x,P[F[i].v[1]].y}; v2 b2=(v2){P[F[i].v[2]].x,P[F[i].v[2]].y}; c=(b1+b2)/2.0,d=b2-b1;d.rot(); b1=c;b2=c+d; v2 O=TAT2D::Cross(a1,a2,b1,b2); d=(v2){P[F[i].v[0]].x,P[F[i].v[0]].y}; db r=(O-d).dis(); if(o.z>0) Ans+=S-TAT2D::intersect(O,r); else Ans+=TAT2D::intersect(O,r); } printf("%.11f\n",(double)(Ans/S+2)); } ```