Fukushima Nuclear Blast - UVa 11978 求圆和多边形面积的并

One of the most disastrous nuclear accidents in history has taken place at Japan’s Fukushima Daiichi nuclear plant. Reactor buildings have been rocked by explosions, caused after damage was sustained from a massive earthquake and tsunami on 11^th march and thus releasing dangerous radiation of unspecified proportions into the air.

The radiation is spreading and it's not a threat only to Japan but also to other countries. If the level of radiation reaches a high level, it will surely cause harms to human health as well as the environment.

You, one of the great programmers in the city, have planned to find the zones that are vulnerable to radiation. If you can measure the level of radiation in a certain city, further actions can be taken to prevent the people from radiation related health consequences. So, at first you want to find the time in which a certain percentage of an area is under radiation threat.

So, at first you modeled the map of a city in 2D space as a simple polygon [1]. You denoted the origin of the explosion as a single point. From this origin, the radiation spreads circularly. You plotted the map such that in each unit of time the radius of the radiation grows one unit. Now, you want to find the time when P% of the area of a city is under threat of radiation.

Input

Input starts with an integer T (≤ 70), denoting the number of test cases.

Each test case starts with a blank line. The next line contains an integer n (3 ≤ n ≤ 5000) denoting the number of vertices of the polygon. Each of the next n lines will contain two integers x_i, y_i denoting the co-ordinate of a vertex of the polygon. The vertices will be given in anticlockwise order. The next line contains 3 integers r_x, r_y and P (0 < P ≤ 100), where (r_x, r_y) denotes the co-ordinate of the origin of explosion and P denotes the percentage. All values for the co-ordinates are between -200 to 200 (inclusive).

Output

For each case, print the case number and the time when exactly P% of the total area is under threat of radiation. Round the time to nearest integer.

Sample Input

Output for Sample Input

2

3

-5 0

5 0

0 5

0 0 100

4

0 0

5 0

3 1

5 6

0 0 17

Case 1: 5

Case 2: 2

题意：给你一个多边形，给你圆的圆心，问当圆的半径为多少时，覆盖多边形面积的P%。

思路：首先求出多边形面积的大小，然后二分圆的半径，每次求一个圆和多边形面积的并。

AC代码如下：

#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>#include<vector>#define double long doubleusing namespace std;const double PI=acos(-1.0);const double eps=1e-10;struct Point{    double x,y;    Point(double x=0,double y=0):x(x),y(y){}};typedef Point Vector;Vector operator + (Vector A,Vector B){return Vector(A.x+B.x,A.y+B.y);}Vector operator - (Vector A,Vector B){return Vector(A.x-B.x,A.y-B.y);}Vector operator * (Vector A,double p){return Vector(A.x*p,A.y*p);}Vector operator / (Vector A,double p){return Vector(A.x/p,A.y/p);}int dcmp(double x){return (x>eps)-(x<-eps);}bool operator <(const Point &A,const Point &B){return dcmp(A.x-B.x)<0 || (dcmp(A.x-B.x)==0 && dcmp(A.y-B.y)<0);}bool operator == (const Point &A,const Point B){return dcmp(A.x-B.x)==0 && dcmp(A.y-B.y)==0;}double Dot(Vector A,Vector B){return A.x*B.x+A.y*B.y;}double Length(Vector A){return sqrt(Dot(A,A));}double Angle(Vector A,Vector B){return acos(Dot(A,B)/Length(A)/Length(B));}double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}double Area2(Point A,Point B,Point C){return Cross(B-A,C-A);}Point GetLineIntersection(Point P,Vector v,Point Q,Vector w){    Vector u=P-Q;    double t=Cross(w,u)/Cross(v,w);    return P+v*t;}double DistanceToLine(Point P,Point A,Point B){    Vector v1=B-A,v2=P-A;    return fabs(Cross(v1,v2)/Length(v1));}double DistanceToSegment(Point P,Point A,Point B){    if(A==B) return Length(P-A);    Vector v1=B-A,v2=P-A,v3=P-B;    if(dcmp(Dot(v1,v2))<0)      return Length(v2);    else if(dcmp(Dot(v1,v3))>0)      return Length(v3);    else      return fabs(Cross(v1,v2)/Length(v1));}Point GetLineProjection(Point P,Point A,Point B){    Vector v=B-A;    return A+v*(Dot(v,P-A)/Dot(v,v));}bool OnSegment(Point p,Point a1,Point a2){    return dcmp(Cross(a1-p,a2-p))==0 && dcmp(Dot(a1-p,a2-p))<=0;}double PolygonArea(Point *p,int n){    double area=0;    for(int i=1;i<n-1;i++)    {        area+=Area2(p[0],p[i],p[i+1]);        //printf("i=%d %.2f\n",i,area);    }    return area/=2;}struct Circle{    Point c;    double r;    Circle(Point c=Point(0,0),double r=0):c(c),r(r){}    Point point(double a){return Point(c.x+r*cos(a),c.y+r*sin(a));}};struct Line{    Point p;    Vector v;    Line(Point p=Point(0,0),Vector v=Vector(0,1)):p(p),v(v){}    Point point(double t){return Point(p+v*t);}};int GetLineCircleIntersection(Line L,Circle C,double &t1,double &t2,vector<Point> &sol){    double a=L.v.x;    double b=L.p.x-C.c.x;    double c=L.v.y;    double d=L.p.y-C.c.y;    double e=a*a+c*c;    double f=2*(a*b+c*d);    double g=b*b+d*d-C.r*C.r;    double delta=f*f-4*e*g;    if(dcmp(delta)==0)    {        t1=t2=-f/(2*e);        sol.push_back(L.point(t1));        return 1;    }    else    {        t1=(-f-sqrt(delta))/(2*e);        t2=(-f+sqrt(delta))/(2*e);        sol.push_back(L.point(t1));        sol.push_back(L.point(t2));        return 2;    }}Point O=Point(0,0);double common_area(Circle C,Point A,Point B){    if(A==B) return 0;    if(A==C.c || B==C.c) return 0;    double OA=Length(A-C.c),OB=Length(B-C.c);    double d=DistanceToLine(O,A,B);    int sg=dcmp(Cross(A,B));    if(sg==0) return 0;    double angle=Angle(A,B);    if(dcmp(OA-C.r)<=0 && dcmp(OB-C.r)<=0)      return Cross(A,B)/2;    else if(dcmp(OA-C.r)>=0 && dcmp(OB-C.r)>=0 && dcmp(d-C.r)>=0)      return sg*C.r*C.r*angle/2;    else if(dcmp(OA-C.r)>=0 && dcmp(OB-C.r)>=0 && dcmp(d-C.r)<0)    {        Point prj=GetLineProjection(O,A,B);        if(OnSegment(prj,A,B))        {            vector<Point> p;            Line L=Line(A,B-A);            double t1,t2;            GetLineCircleIntersection(L,C,t1,t2,p);            double s1=C.r*C.r*angle/2;            double s2=C.r*C.r*Angle(p[0],p[1])/2;            s2-=fabs(Cross(p[0],p[1])/2);            s1-=s2;            return sg*s1;        }        else          return sg*C.r*C.r*angle/2;    }    else    {        if(dcmp(OB-C.r)<0)        {            Point temp=A;            A=B;            B=temp;        }        Point inter_point;        double t1,t2;        Line L=Line(A,B-A);        vector<Point> inter;        GetLineCircleIntersection(L,C,t1,t2,inter);        if(OnSegment(inter[0],A,B))          inter_point=inter[0];        else          inter_point=inter[1];        double s=fabs(Cross(inter_point,A)/2);        s+=C.r*C.r*Angle(inter_point,B)/2;        return s*sg;    }}int T,t,n;Point p[5010];double Area,P,S;int main(){    int i,j,k;    double l,r,mi;    scanf("%d",&T);    for(t=1;t<=T;t++)    {        scanf("%d",&n);        for(i=0;i<n;i++)           scanf("%Lf%Lf",&p[i].x,&p[i].y);        p[n]=p[0];        scanf("%Lf%Lf%Lf",&O.x,&O.y,&P);        for(i=0;i<=n;i++)           p[i]=p[i]-O;        O=Point(0,0);        Area=fabs(PolygonArea(p,n))*P/100;        l=0;r=10000;        while(r-l>0.0001)        {            mi=(l+r)/2;            Circle C(O,mi);            S=0;            for(i=0;i<n;i++)               S+=common_area(C,p[i],p[i+1]);            S=fabs(S);            if(dcmp(S-Area)<0)              l=mi;            else              r=mi;        }        printf("Case %d: %.0Lf\n",t,l);    }}