20200722 T2 ACT4！无限回转！【辛普森积分】

题目描述

$n$ 个点的凸多边形，逆时针方向给出，在内部随机一个起点，随机一个方向，一直往这个方向走直到碰到某一条边。求每条边被碰到的概率。

$n\le 1500,|x_i|,|y_i|\le 10^4$

题目分析

对于第 $i$ 条边，求它被碰到的概率。

随机过程有两个，一个是坐标，一个是角度，考虑积分，最后除以 (面积$\ast 2\pi$)

直接按坐标积分 $dx,dy$，用 $atan$ 求角度范围，两层积分复杂度太高

按照中垂线确定角度大小，算出圆弧的长度，积分角度，算圆弧复杂度太高。

考虑对于一条边，确定一个角度，那么从这个方向可以碰到这条边的面积是这样的：

假设这个面积是 $S(\theta )$，就是要求 $\int S(\theta )d\theta$

差分之后，就是求一条线右边与多边形交的面积。可以对其三角剖分，变成一个前缀和加上一个小三角形：

于是求一次的复杂度是 $O(\log n)$，然后就可以大力辛普森积分了。差分的减法可以直接利用下一个点的答案。

直接整可能会 $T$，需要一点常数优化，发现 $n,|x_i|,|y_i|$ 越大，小块的面积影响就越小，所以当 $n$ 大时 $eps$ 开小一点，$n=1500$ 时 $eps=1e-2$ 精度都很高。

复杂度 $O(n\ast Simpson(eps)\ast log)$

Code：

#include<bits/stdc++.h>
#define maxn 1505
using namespace std;
const double Pi = acos(-1);
int n;
double s[maxn],ans[maxn],lim[maxn],eps;
struct Point{
	double x,y; Point(double x=0,double y=0):x(x),y(y){}
	Point operator + (const Point &p){return Point(x+p.x,y+p.y);}
	Point operator - (const Point &p){return Point(x-p.x,y-p.y);}
	double operator * (const Point &p)const{return x*p.y-y*p.x;}
	double operator ^ (const Point &p)const{return x*p.x+y*p.y;}
	double len(){return sqrt(x*x+y*y);}
	bool operator < (const Point &p)const{return (*this)*p>0;}
}a[maxn*2],p[maxn],base;
double calc(double ang){
	double Sin=sin(ang),Cos=cos(ang);
	Point t = Point(base.x*Cos-base.y*Sin,base.y*Cos+base.x*Sin);
	int k = lower_bound(p+1,p+n,t)-p-1;
	return k==0 ? 0 : k==n-1 ? s[n-1] : s[k]+(p[k]*t)/(p[k]*t+t*p[k+1])*(p[k]*p[k+1]);
}
double Simpson(double l,double m,double r,double L,double M,double R,int dep){
	double m1=(l+m)/2,m2=(m+r)/2,M1=calc(m1),M2=calc(m2);
	double S=(L+4*M+R)/6*(r-l), S1=(L+4*M1+M)/6*(m-l), S2=(M+4*M2+R)/6*(r-m);
	if(fabs(S-S1-S2)<eps&&dep>=3) return S;
	return Simpson(l,m1,m,L,M1,M,dep+1)+Simpson(m,m2,r,M,M2,R,dep+1);
}
int main()
{
	freopen("tusk.in","r",stdin);
	freopen("tusk.out","w",stdout);
	scanf("%d",&n),eps=n<=5?1e-9:n<=100?1e-6:1e-2;
	for(int i=1;i<=n;i++) scanf("%lf%lf",&a[i].x,&a[i].y),a[n+i]=a[i];
	for(int i=1;i<=n;i++){
		for(int j=1;j<n;j++) p[j]=a[i+j]-a[i];
		for(int j=2;j<n;j++) s[j]=s[j-1]+p[j-1]*p[j];
		lim[i]=acos((p[1]^p[n-1])/p[1].len()/p[n-1].len());
		base=p[1];
		double L=calc(0),M=calc(lim[i]/2),R=calc(lim[i]);
		ans[i]=Simpson(0,lim[i]/2,lim[i],L,M,R,0);
	}
	double tmp=ans[1];
	for(int i=1;i<n;i++) ans[i]-=ans[i+1];
	ans[n]-=tmp;
	double S=0;
	for(int i=1;i<=n;i++) S+=a[i]*a[i+1];
	for(int i=1;i<=n;i++) ans[i]+=(Pi-lim[i])*S;
	for(int i=1;i<=n;i++) printf("%.10f\n",ans[i]/(2*Pi*S));
}

标算是再把面积分到每条边上，只需要积分那个小三角形的面积

还有一种做法是均分角度，每次将所有点旋转一个角度，然后每个点引一条射线，可以得到每个面积贡献到哪条边上，速度飞快，如图所示：

Code：

//code by dsfzsunjiawei
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<ctime>
#include<cstring>
#include<random>
#include<cmath>
#define errlog(...) fprintf(stderr,__VA_ARGS__)
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
std::mt19937 Rand((unsigned ll)(new int*));
inline ull Ranll()
{return(Rand()*(1llu<<32llu))+(Rand());}
inline ll random(ll l,ll r)
{return Ranll()%(r-l+1)+l;}
#define retrun return
#define cpnst const
int n;
int x[11111],y[11111];
struct pnt
{db x,y;pnt(db _x=0,db _y=0){x=_x,y=_y;}}r[11111];
db ans[11111];
db Fs(db L,db R,db H)
{return H*(R+L)/2;}
cpnst db eps=1e-8;
bool isequal(db x,db y){return x-y<eps&&y-x<eps;}
db Ht(pnt p1,pnt p2,db x)
{
	if(isequal(p1.x,x))return p1.y;
	if(isequal(p2.x,x))return p2.y;
	retrun p1.y+(p2.y-p1.y)*(x-p1.x)/(p2.x-p1.x);
}
void solve()
{
	register int i;
	db mnx=1e9;
	int mp=1,L,R;
	r[0]=r[n],r[n+1]=r[1];
	for(i=1;i<=n;i++)if(mnx>r[i].x)mnx=r[i].x,mp=i;
	db X=r[mp].x,YU=r[mp].y,YD=YU,Ls=0;
	L=mp,R=mp;
	for(i=2;i<=n;i++)
	{
		if(r[L-1].x<eps+r[R+1].x)
		{
			db det=r[L-1].y-Ht(r[R],r[R+1],r[L-1].x);
			ans[R]+=Fs(Ls,det,r[L-1].x-X);
			Ls=det,X=r[L-1].x,L=(L+n-2)%n+1;
		}else
		{
			db det=Ht(r[L],r[L-1],r[R+1].x)-r[R+1].y;
			ans[R]+=Fs(Ls,det,r[R+1].x-X);
			Ls=det,X=r[R+1].x,R=R%n+1;
		}
	}
}
void rota(db theta)
{
	register int i;
	db S=sin(theta),C=cos(theta);
	for(i=1;i<=n;i++)r[i]=pnt(r[i].x*C-r[i].y*S,r[i].y*C+r[i].x*S);
}
cpnst db pi=acos(-1);
int main()
{
	#ifndef StdIO
	freopen("tusk.in","r",stdin);
	freopen("tusk.out","w",stdout);
	#endif
	scanf("%d",&n);
	register int i;
	for(i=1;i<=n;i++)scanf("%d%d",x+i,y+i),r[i]=pnt(x[i],y[i]);
	for(i=1;i<=4954;i++)rota(pi/2477),solve();
	db SS=0;
	for(i=1;i<=n;i++)SS+=ans[i];
	for(i=1;i<=n;i++)printf("%.15lf\n",ans[i]/SS);
}