http://acm.hdu.edu.cn/showproblem.php?pid=4327
题意:给你一些点( x, y ), 0 < x < 1 , 0 < y < 1.再给你一个概率分布函数,表示每个坐标点的概率分布,显然整个矩形的分布概率之和为1
另外有一个规则:每个点都只管辖离自己最近的区域,即和其他所有点的中垂线所切割成的区域。
求每块区域分布的概率是多少
解法:显然,每个点所管辖的区域可以用这个点与其他点得中垂线去切割平面,这个点所在的一侧为有效半平面,所以要先判断点与直线的方位再去切割平面,这里我写了两个切割的函数,
容易看出,题目所给的概率分布函数是空间中的一个平面,所以可以想象以半平面切割后的区域为底,最顶上是一个斜面的立体(这个空间立体要仔细想想)
那么这个立体的体积就是我们要求的答案。
具体在求的时候可以先求出底下棱柱的体积,最上面的那个角单独计算。
咋一看,最上面的部分什么形状都没有,没办法算啊,仔细一看就知道了,可以将它切割成若干个四棱锥,锥顶始终是最矮的那个点,这个需要大家发挥空间想象能力。
这个题搞了我一个下午+晚上的时间啊,最后发现错在半平面交 模板,我的模板是n^2的,以前写的时候都没有注意一个很小的细节,在这道题上暴露出来了,发现问题,及时解决总是好的,go on AC!
ps:半平面交学习的话可以去这里,写的很好
#include <cmath>
#include <cstdio>
#include<algorithm>
using namespace std;
const double INF = 1000000000.0;
const int maxn = 1010;
const double eps = 1e-12;
inline double sgn(double x) {return fabs(x)<eps?0:(x>0?1:-1);}
struct Point{
double x,y;
Point(double tx=0,double ty=0){x=tx;y=ty;}
bool operator == (const Point& t) const {
return sgn(x-t.x)==0 && sgn(y-t.y)==0;
}
Point operator - (const Point& t) const {
Point tmp;
tmp.x = x - t.x;
tmp.y = y - t.y;
return tmp;
}
}p[maxn],tmp[maxn],pp[maxn],GP;
struct Seg{Point s,e;};
struct Line {
double a, b, c;
};//中垂线
double cross(Point a,Point b,Point c){return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);}
inline Point intersect(Point x,Point y,double a,double b,double c){
double u = fabs(a * x.x + b * x.y + c);
double v = fabs(a * y.x + b * y.y + c);
return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}
int n,ban_tot;
void CUT1(double a,double b,double c){//points must be clockwise a*x+b*y+c>=0表示点在半平面内
int i,tot=0;
for(int i = 1; i <= ban_tot; ++i){
if(a*p[i].x + b*p[i].y + c >= eps) pp[++tot] = p[i];
else {
if(a*p[i-1].x + b*p[i-1].y + c > eps){
pp[++tot] = intersect(p[i],p[i-1],a,b,c);
}
if(a*p[i+1].x + b*p[i+1].y + c > eps){
pp[++tot] = intersect(p[i],p[i+1],a,b,c);
}
}
} ban_tot=tot;
pp[tot+1]=pp[1];pp[0]=pp[tot];
memcpy(p,pp,sizeof(pp));
}
void CUT2(double a,double b,double c){//a*x+b*y+c<=0表示点在半平面内
int i,tot=0;
for(int i = 1; i <= ban_tot; ++i){
if(!(a*p[i].x + b*p[i].y + c > eps) ) pp[++tot] = p[i]; //点在半平面内或者边界上
else {
if(a*p[i-1].x + b*p[i-1].y + c < -eps){ //上个点完全在半平面内
pp[++tot] = intersect(p[i],p[i-1],a,b,c);
}
if(a*p[i+1].x + b*p[i+1].y + c < -eps){ //下个点完全在半平面内
pp[++tot] = intersect(p[i],p[i+1],a,b,c);
}
}
} ban_tot=tot;
pp[tot+1]=pp[1];pp[0]=pp[tot];//两端都扩展一个点
memcpy(p,pp,sizeof(pp));
}
Line Turn(Point s, Point e) { // 线段转直线表达式
Line ln;
ln.a = s.y - e.y;
ln.b = e.x - s.x;
ln.c = s.x*e.y - e.x*s.y;
return ln;
}
Line make(Point a,Point b)//求a,b中垂线的直线方程
{
double x0=(a.x+b.x)/2;
double y0=(a.y+b.y)/2;
Line tmp=Turn(a,b);
Line ans;
ans.a=tmp.b;
ans.b=-tmp.a;
ans.c=tmp.a*y0-tmp.b*x0;
return ans;
}
Line ln[maxn];
inline double PPdis(Point a, Point b) { // 点点距离
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline double PLdis(Point p,Point l1,Point l2){ // 点线距离
return fabs(cross(p,l1,l2))/PPdis(l1,l2);
}
double calc(Point *p,int n)
{
if(n<3) return 0;
double area=0,V=0;
for(int i=0;i<n;i++) area+=p[(i+1)%n].y*(p[i].x-p[(i+2)%n].x);
area/=2;
area=fabs(area);
double H=10;
int pos=0;
for(int i=0;i<n;i++)
{
if(2-p[i].x-p[i].y<H)
{
H=2-p[i].x-p[i].y;
pos=i;
}
}
V+=area*H;
for(int i=pos+1;;i++)
{
int id1=i%n;
int id2=(i+1)%n;
if(id2==pos) break;
double h=PLdis(p[pos],p[id1],p[id2]);
double s=((2-p[id1].x-p[id1].y-H) + (2-p[id2].x-p[id2].y-H)) * PPdis(p[id1],p[id2])/2;
V+=s*h/3;
}
return V;
}
double ans[110];
int main(){
int t,ca=1,n;
scanf("%d",&t);
while(t--)
{
scanf("%d",&n);
for(int i=1;i<=n;i++) scanf("%lf%lf",&tmp[i].x,&tmp[i].y);
for(int i=1;i<=n;i++)
{
p[1]=Point(0,0);
p[2]=Point(0,1);
p[3]=Point(1,1);
p[4]=Point(1,0);
p[0]=p[4];
p[5]=p[1];//首尾都要加上一个点,保证边界的合法点不会漏掉
ban_tot=4;
double x0=tmp[i].x,y0=tmp[i].y;
for(int j=1;j<=n;j++)
{
if(i==j) continue;
Line mid_line = make(tmp[i],tmp[j]);
double a=mid_line.a,b=mid_line.b,c=mid_line.c;
if(a*x0+b*y0+c >= eps)
{
CUT1(a,b,c);
}
else
{
CUT2(a,b,c);
}
}
double tmpv=calc(p,ban_tot);
ans[i]=tmpv;
}
printf("Case #%d:\n",ca++);
for(int i=1;i<=n;i++)
{
printf("%.6lf\n",ans[i]);
}
}
return 0;
}