某 SCOI 模拟赛 T3 无限回转（tusk）【辛普森积分】

题意

给出一个凸 $N$ 边形，在其内部随机选一个点，再随机选择一个方向，求从该点出发沿这个方向一直走最后碰到各条边的概率。要求绝对误差不超过 $10^{-6}$。$N\le 1500$，坐标范围 $10^4$。时限 1s。

（题目改编自 PE 613）

题解

枚举每条边。假如确定了方向，考虑有多大的面积朝这个方向走会走到这条边：

不断转两条射线，也就对应不同的角度时的概率。

两条直线不好处理，我们算其中一条直线左侧面积再算另一条最后相减会方便许多。

标算采用的是数学推导，对角度积分面积。

实际上通过均匀取角度求平均或者辛普森积分同样能做。先把所有点按照相对于射线的旋转中心的极角排序，并记录三角形面积的前缀和，然后确定角度后二分到与射线相交的线段，把前缀面积与多出来的一小块相加即为该角度对应的面积。

代码（辛普森积分）：

#include<bits/stdc++.h>
using namespace std;
int getint(){
	int ans=0,f=1;
	char c=getchar();
	while(c<'0'||c>'9'){
		if(c=='-')f=-1;
		c=getchar();
	}
	while(c>='0'&&c<='9'){
		ans=ans*10+c-'0';
		c=getchar();
	}
	return ans*f;
}
const int N=1510;
const double pi=3.14159265358979323846264,eps=1e-10;
bool eq(double a,double b){
	return abs(a-b)<eps;
}

struct point{
	double x,y;
	point(double x=0,double y=0):x(x),y(y){}
};
typedef point vec;
vec operator- (point a,point b){ return vec(a.x-b.x,a.y-b.y); }
point operator+ (point a,vec b){ return point(a.x+b.x,a.y+b.y); }
double cross(vec a,vec b){ return a.x*b.y-b.x*a.y; }
double dot(vec a,vec b){ return a.x*b.x+a.y*b.y; }
double atan2(point a){ return atan2(a.y,a.x); }
ostream& operator<< (ostream &os,const point p){
    os<<"("<<p.x<<","<<p.y<<")";
    return os;
}

struct line{
	point a,b;
	double k;
	line(){}
	line(point a,point b):a(a),b(b),k(atan2(b-a)){}
};
point inter(line a,line b){
	double s1=cross(b.a-a.a,b.b-a.a),s2=cross(b.b-a.b,b.a-a.b);
	if(eq(s1+s2,0)){
		cerr<<a.a<<" "<<a.b<<" | "<<b.a<<" "<<b.b<<endl;
		assert(0);
	}
	return point(a.a.x+(a.b.x-a.a.x)/(s1+s2)*s1,a.a.y+(a.b.y-a.a.y)/(s1+s2)*s1);
}
bool isleft(line l,point p){
	return cross(l.b-l.a,p-l.a)>-eps;
}

int n;
point p[N];
double ang[N],sum[N];
bool boo;

double calc(double x){
	int q=upper_bound(ang,ang+n,x)-ang-1;
	//cerr<<"calc "<<x<<" "<<q<<"  "<<ang[q]<<"? "<<ang[q+1]<<endl;
	double ans=sum[q];
	//if(q>=n-1||eq(x,ang[n-1]))cerr<<"[ "<<x<<" "<<x/pi*180<<" "<<ans<<endl;
	if(q>=n-1||eq(x,ang[n-1]))return ans;
	//if(q==0||eq(x,ang[1]))cerr<<"] "<<x<<" "<<x/pi*180<<" "<<0<<endl;
	if(q==0||eq(x,ang[1]))return 0;
	line l=line(p[0],p[0]+vec(cos(x),sin(x)));
	point c=inter(line(p[q],p[(q+1)%n]),l);
	//cerr<<"| "<<x<<" "<<q<<" "<<c<<" "<<ans<<"+"<<cross(c-p[0],p[q]-p[0])<<endl;
	ans+=cross(c-p[0],p[q]-p[0]);
	//cerr<<"> "<<x<<" "<<x/pi*180<<" "<<ans<<endl;
	return ans;
}
double calc(double dx,double y1,double y2,double y3){
	return (y1+4*y2+y3)*dx/6;
}
double simpeps=1e-7;
double simp(double l,double r,double mid,double yl,double yr,double ymid,double c){
	double lm=(l+mid)/2,ylm=calc(lm);
	double rm=(r+mid)/2,yrm=calc(rm);
	double c1=calc(mid-l,yl,ylm,ymid),c2=calc(r-mid,ymid,yrm,yr);
	//cerr<<"> "<<l<<" "<<mid<<" "<<r<<" "<<c1<<" "<<c2<<" "<<c<<endl;
	if(abs(c1+c2-c)<=simpeps)return c1+c2;
	return simp(l,mid,lm,yl,ymid,ylm,c1)+simp(mid,r,rm,ymid,yr,yrm,c2);
}
#define simp(l,r) simp(l,r,(l+r)/2,calc(l),calc(r),calc((l+r)/2),\
					   calc(r-l,calc(l),calc(r),calc((l+r)/2)))

int main(){
	n=getint();
	for(int i=0;i<n;i++)p[i].x=getint(),p[i].y=getint();
	double s=0;
	for(int i=0;i<n;i++)s+=cross(p[i],p[(i+1)%n]);
	simpeps=eps*s;
	for(int i=0;i<n;i++){
		//if(i%10==0)cerr<<i<<endl;
		//cerr<<"---- "<<i<<" ----"<<endl;
		for(int i=1;i<n;i++){
			ang[i]=atan2(p[i]-p[0]);
			if(i>1){
				while(ang[i]<ang[1])ang[i]+=pi*2;
				while(ang[i]>ang[1]+pi)ang[i]-=pi*2;
			}
			sum[i]=sum[i-1]+cross(p[i]-p[0],p[i-1]-p[0]);
		}
		sum[n]=sum[n-1];
		ang[0]=-pi*2;
		double simp1=simp(ang[1],ang[1]+pi);
		//cerr<<">> "<<simp1/pi/2<<endl;
		
		//cerr<<"|--- "<<i<<" ---|"<<endl;
		boo=1;
		p[n]=p[0];
		for(int i=0;i<n;i++)p[i]=p[i+1];
		ang[n-1]=atan2(p[n-1]-p[0]);
		//cerr<<"ang "<<n-1<<" "<<ang[n-1]<<endl;
		for(int i=1;i<n-1;i++){
			ang[i]=atan2(p[i]-p[0]);
			while(ang[i]<ang[n-1])ang[i]+=pi*2;
			while(ang[i]>ang[n-1]+pi)ang[i]-=pi*2;
			sum[i]=sum[i-1]+cross(p[i]-p[0],p[i-1]-p[0]);
			//cerr<<"ang "<<i<<" "<<ang[i]<<endl;
		}
		sum[n]=sum[n-1];
		ang[0]=-2*pi;
		double simp2=-simp(ang[n-1],ang[n-1]-pi);
		//cerr<<">> "<<simp2/pi/2<<endl;
		
		printf("%.10lf\n",(simp2-simp1)/s/pi/2);
	}
	return 0;
}