2021CCPC桂林站 F-Illuminations II(计算几何)

2021CCPC桂林站 F-Illuminations II

Problem Statement

你得到两个凸包, 外凸包有$n$个点, 内凸包有$m$个点. 保证内凸包严格在外凸包内. 现在我们在外凸包边界放一个光源, 求照亮内凸包边界的期望长度.

Solution

由于期望的线性性: 我们考虑求内凸包每一条边的贡献.

• 下面:我们不妨设用两个点来表示一条边: 例如$(u_i,v_i)$表示起点为$u_i$终点为$v_i$的边, 这里表述和点坐标不同. 计算几何中也常用两个点来表示一个直线. 下面我们就用这种表示法来表述.

设内凸包(逆时针顺序记录点)上一个边$(u_i,v_i)$, 外凸包的周长为$\text{Sum}$, $(u_i,v_i)$与外凸包的交点为$(a_i,b_i)$, 那么外凸包$(a_i,b_i)$之间且在$(u_i,v_i)$右侧的折线距离设为$Le{n}_{(a_i,b_i)}$. 那么边$(u_i,v_i)$对答案的贡献就是. $E_{(u_i,v_i)}=\frac{di{s}_{(u_i,v_i)}\times Le{n}_{(a_i,b_i)}}{Sum}$. 其中你可以理解为随机变量为$X_{(u_i,v_i)}=di{s}_{(u_i,v_i)}$其概率为$P_{(u_i,v_i)}=\frac{Le{n}_{(a_i,b_i)}}{Sum}$. 那么有期望公式即可得到$E_{(u_i,v_i)}=X_{(u_i,v_i)}{P}_{(u_i,v_i)}$.

那么答案就是: $Ans={\sum }_{i=1}^m{E}_{(u_i,v_i)}=\frac{1}{Sum}{\sum }_{i=1}^m[di{s}_{(u_i,v_i)}\times Le{n}_{(a_i,b_i)}]$.

• 注意请使用long double.
• 注意内凸包可能与外凸包有交点, 请使用叉积判断是否与外凸包相交.
• 特别地: 对于内凸包与外凸包的边相交, 当且仅当内凸包的这一条边与外凸包对应边的叉积值是异号的.

Code(Based On Kuangbin)

# define Fast_IO std::ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
# include "unordered_map"
# include "algorithm"
# include "iostream"
# include "cstdlib"
# include "cstring"
# include "cstdio"
# include "vector"
# include "bitset"
# include "queue"
# include "cmath"
# include "map"
# include "set"

using namespace std;

const double eps = 1e-8;
const double pi = acos(-1.0);
const int maxp =1e5+10;

int sgn(double x){return (fabs(x)<eps)?0:(x<0?-1:1);}
inline double sqr(double x){return x*x;}
struct Point{
	double x,y;
	Point(){}
	Point(double _x,double _y){x = _x,y = _y;}
	void input(){scanf("%lf%lf",&x,&y);}
	void output(){printf("%.2lf %.2lf\n",x,y);}
	bool operator == (Point b)const{return sgn(x-b.x)==0 && sgn(y-b.y)==0;}
	bool operator < (Point b)const{return sgn(x-b.x)== 0?sgn(y-b.y)<0:x<b.x;}
	Point operator -(const Point &b)const{return Point(x-b.x,y-b.y);} //叉积
	double operator ^(const Point &b)const{return x*b.y - y*b.x;} //点积
	double operator *(const Point &b)const{return x*b.x + y*b.y;} //返回长度
	double len(){return hypot(x,y);} //库函数->返回长度的平方
	double len2(){return x*x + y*y;} //返回两点的距离
	double distance(Point p){return hypot(x-p.x,y-p.y);}
	Point operator +(const Point &b)const{return Point(x+b.x,y+b.y);}
	Point operator *(const double &k)const{return Point(x*k,y*k);}
	Point operator /(const double &k)const{return Point(x/k,y/k);}
	//计算pa和pb的夹角->就是求这个点看a,b 所成的夹角
	//测试 LightOJ1203
	double rad(Point a,Point b){Point p=*this;return fabs(atan2(fabs((a-p)^(b-p)),(a-p)*(b-p)));}
	//化为长度为r的向量
	Point trunc(double r){return !sgn(len())?*this:Point(x*r,y*r)/len();}
	Point rotleft(){return Point(-y,x);} //逆时针旋转90度
	Point rotright(){return Point(y,-x);} //顺时针旋转90度
	Point rotate(Point p,double angle){ //绕着p点逆时针旋转angle
		Point v=(*this)-p;
		double c=cos(angle),s=sin(angle);
		return Point(p.x+v.x*c-v.y*s,p.y+v.x*s+v.y*c);
	}
};

struct Line{
	Point s,e;
	Line(){}
	Line(Point _s,Point _e){s=_s,e=_e;}
	bool operator==(Line v){return (s == v.s)&&(e == v.e);}
	//`根据一个点和倾斜角angle确定直线,0<=angle<pi`
	Line(Point p,double angle){s=p,e=sgn(angle-pi/2)==0?(s+Point(0,1)):(s+Point(1,tan(angle)));}
	Line(double a,double b,double c){ //ax+by+c=0
		if(sgn(a)==0) s=Point(0,-c/b),e=Point(1,-c/b);
		else if(sgn(b)==0) s=Point(-c/a,0),e=Point(-c/a,1);
		else s=Point(0,-c/b),e=Point(1,(-c-a)/b);
	}
	void input(){s.input();e.input();}
	void adjust(){if(e<s)swap(s,e);}
	double length(){return s.distance(e);} //求线段长度
	double angle(){ //返回直线倾斜角 0<=angle<pi
		double k=atan2(e.y-s.y,e.x-s.x);
		k+=sgn(k)<0?pi:0,k-=sgn(k-pi)==0?pi:0;
		return k;
	}//点和直线关系:1-在左侧;2-在右侧;3-在直线上
	int relation(Point p){
		int c=sgn((p-s)^(e-s));
		return (c<0)?1:(c>0?2:3);
	}//点在线段上的判断
	bool pointonseg(Point p){return sgn((p-s)^(e-s))==0 && sgn((p-s)*(p-e))<=0;}
	//两向量平行(对应直线平行或重合)
	bool parallel(Line v){return sgn((e-s)^(v.e-v.s))==0;}
	//两线段相交判断:2-规范相交;1-非规范相交;0-不相交
	int segcrossseg(Line v){
		static int d1,d2,d3,d4;
		d1=sgn((e-s)^(v.s-s)),d2=sgn((e-s)^(v.e-s));
		d3=sgn((v.e-v.s)^(s-v.s)),d4=sgn((v.e-v.s)^(e-v.s));
		if((d1^d2)==-2 && (d3^d4)==-2)return 2;
		return (d1==0 && sgn((v.s-s)*(v.s-e))<=0) ||
			(d2==0 && sgn((v.e-s)*(v.e-e))<=0) ||
			(d3==0 && sgn((s-v.s)*(s-v.e))<=0) ||
			(d4==0 && sgn((e-v.s)*(e-v.e))<=0);
	}//直线和线段相交判断:2-规范相交;1-非规范相交;0-不相交`
	int linecrossseg(Line v){
		static int d1,d2;
		d1=sgn((e-s)^(v.s-s)),d2=sgn((e-s)^(v.e-s));
		return (d1^d2)==-2?2:(d1==0||d2==0);
	}//两直线关系:0-平行;1-重合;2-相交
	int linecrossline(Line v){return (*this).parallel(v)?v.relation(s)==3:2;}
	//求两直线的交点:要保证两直线不平行或重合
	Point crosspoint(Line v){
		double a1=(v.e-v.s)^(s-v.s);
		double a2=(v.e-v.s)^(e-v.s);
		return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
	}//点到直线的距离
	double dispointtoline(Point p){return fabs((p-s)^(e-s))/length();}//点到线段的距离
	double dispointtoseg(Point p){
		if(sgn((p-s)*(e-s))<0 || sgn((p-e)*(s-e))<0) return min(p.distance(s),p.distance(e));
		return dispointtoline(p);
	}//返回线段到线段的距离:前提是两线段不相交，相交距离就是0了`
	double dissegtoseg(Line v){
		return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));
	}//返回点p在直线上的投影
	Point lineprog(Point p){return s+(((e-s)*((e-s)*(p-s)))/((e-s).len2()));}
	//返回点p关于直线的对称点
	Point symmetrypoint(Point p){Point q=lineprog(p);return Point(2*q.x-p.x,2*q.y-p.y);}
};

const int maxm=2e5+10;
int n,m;
Point P[maxm], p[maxm];
bool chkleft(Point a, Point b, Point c, Point d){return sgn((b-a)^(c-a))<0 && sgn((b-a)^(d-a))>=0;} //判断(c,d)是否在(a,b)左侧 
bool chkright(Point a, Point b, Point c, Point d){return sgn((b-a)^(d-a))>0 && sgn((b-a)^(c-a))<=0;} //判断(c,d)是否在(a,b)右侧 
double Dis[maxm],Sum;

int main(){
	static int i,j,k;
	scanf("%d%d",&n,&m);
	for(i=0;i<=n-1;i++) P[i].input();
	for(i=0;i<=m-1;i++) p[i].input();
	for(i=0;i<=n-1;i++) Dis[i] = P[i].distance(P[(i+1)%n]), Sum += Dis[i];
	double now=0,tot=0;
	i=0,j=0,k=0;
	while(!chkleft(p[(i+1)%m],p[i],P[j],P[(j+1)%n])) j=(j+1)%n;
	k=j;
	double l,r;
	while(i<m){
		l=P[j].distance(Line(p[i], p[(i+1)%m]).crosspoint(Line(P[j],P[(j+1)%n])));
		while(!chkright(p[i], p[(i+1)%m],P[k],P[(k+1)%n])) now+=Dis[k],k=(k+1)%n;
		r= P[k].distance(Line(p[i],p[(i+1)%m]).crosspoint(Line(P[k],P[(k+1)%n])));
		tot+=p[i].distance(p[(i+1)%m])*(now-l+r);i++;
		while(!chkleft(p[(i+1)%m],p[i],P[j],P[(j+1)%n])) now-=Dis[j],j=(j+1)%n;
	}printf("%.10lf",tot/Sum);
	return 0; 
}

Link

Link1: CF Gym 103409 F-Illuminations II

Link2-Kuangbin: Kuangbin

Link3-算法: marx97