zoj 3720 计算几何，期望

链接：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=5043

题意：给一个矩形，矩形内一个多边形，朝矩形内投掷落在某点上，由于磁力的作用，该点会移到离它最近的整点处，如果这个整点（x，y）在多边形内，则有分数S=A*x+B*y，求分数的平均期望。

思路：枚举矩形内的整点，判断它是否在多边形内，如果在，那么有分数，计算出来。由于每个整点的吸收范围是[x-0.5,x+0.5],[y-0.5,y+0.5]，所以相应的概率P为这个小矩形的面积除以大矩形的面积。期望即是Si 乘以P，再求和。期望公式 ∑ i*Pi .

#include<cstdio>
#include<algorithm>
#include<vector>
#include<iostream>
#include<cmath>
#include<cstdlib>
using namespace std;
const int maxn=20+5;
struct point
{
    int x,y;
    point(int x=0,int y=0):x(x),y(y){}
    point operator+(const point& a)const {return point(x+a.x, y+a.y);}
    point operator-(const point& a)const {return point(x-a.x, y-a.y);}
    point operator*(const int& k)const {return point(x*k, y*k);}
    int operator*(const point& a)const {return x*a.x + y*a.y;}
    int operator^(const point& a)const {return x*a.y - a.x*y;}
    int cross(point a,point b){return a-*this ^ b-*this;}
    int cross(point b,point c,point d){return b-*this ^ d-c;}
}v[maxn],v1[maxn];
double P0,Q0,P1,Q1;
int n,A,B;
int cmp1(point a,point b)
{
    if( a.x != b.x )
       return a.x<b.x;
    else
       return a.y<b.y;
}
int cmp2(point a,point b)
{
    if( a.y != b.y )
       return a.y<b.y;
    else
       return a.x<b.x;
}
int InPolygon(point a)
{
    int nCross = 0;
    for(int i = 0 ; i < n; i++) {
        int j = (i+1) % n;
        if((a.y-v1[i].y)*(a.y-v1[j].y) <= 0) {
           int dir = v1[j].cross(v1[i] , a);
           if(dir==0 &&
              (a.x>=v1[i].x && a.x<=v1[j].x || a.x >= v1[j].x && a.x <= v1[i].x)) return 1;
           if(dir > 0 && v1[i].y > v1[j].y && a.y != v1[i].y ||
              dir < 0 && v1[i].y < v1[j].y && a.y != v1[j].y)  nCross++;
        }
    }
    return nCross&1;
}
void sol()
{
    double ans=0.0,area=(P1-P0)*(Q1-Q0);
    sort(v,v+n,cmp1);int lx=v[0].x,rx=v[n-1].x;
    sort(v,v+n,cmp2);int ly=v[0].y,ry=v[n-1].y;
    for(int i=lx;i<=rx;i++)
       for(int j=ly;j<=ry;j++)
       {
           point p=point(i,j);
           if(InPolygon(p))
           {
              double l=min(P1,i+0.5)-max(P0,i-0.5);
              double w=min(Q1,j+0.5)-max(Q0,j-0.5);
              ans+=(A*i+B*j)*l*w;
           }
       }
     printf("%.3lf\n",ans/area);
}

int main()
{
    while(~scanf("%lf%lf%lf%lf",&P0,&Q0,&P1,&Q1))
    {
        scanf("%d%d%d",&n,&A,&B);
        for(int i=0;i<n;i++)
        {
            scanf("%d%d",&v[i].x,&v[i].y);
            v1[i].x=v[i].x,v1[i].y=v[i].y;
        }
        sol();
    }
    return 0;
}