题意:给出小于10个点形成的凸多边形 和一个半径为r 可以移动的圆 求圆心在何处的面积交最大,面积为多少
题解:三分套三分求出圆心位置,再用圆与多边形面积求交
1 #include<bits/stdc++.h> 2 #define inf 1000000000000 3 #define M 100009 4 #define eps 1e-12 5 #define PI acos(-1.0) 6 using namespace std; 7 struct node 8 { 9 double x,y; 10 node(){} 11 node(double xx,double yy) 12 { 13 x=xx; 14 y=yy; 15 } 16 node operator -(node s) 17 { 18 return node(x-s.x,y-s.y); 19 } 20 node operator +(node s) 21 { 22 return node(x+s.x,y+s.y); 23 } 24 double operator *(node s) 25 { 26 return x*s.x+y*s.y; 27 } 28 double operator ^(node s) 29 { 30 return x*s.y-y*s.x; 31 } 32 }p[M],a[M]; 33 double max(double a,double b) 34 { 35 return a>b?a:b; 36 } 37 double min(double a,double b) 38 { 39 return a<b?a:b; 40 } 41 double len(node a) 42 { 43 return sqrt(a*a); 44 } 45 double dis(node a,node b)//两点之间的距离 46 { 47 return len(b-a); 48 } 49 double cross(node a,node b,node c)//叉乘 50 { 51 return (b-a)^(c-a); 52 } 53 double dot(node a,node b,node c)//点乘 54 { 55 return (b-a)*(c-a); 56 } 57 int judge(node a,node b,node c)//判断c是否在ab线段上（前提是c在直线ab上） 58 { 59 if(c.x>=min(a.x,b.x) 60 &&c.x<=max(a.x,b.x) 61 &&c.y>=min(a.y,b.y) 62 &&c.y<=max(a.y,b.y)) 63 return 1; 64 return 0; 65 } 66 double area(node b,node c,double r) 67 { 68 node a(0.0,0.0); 69 if (dis(b,c)<eps) return 0.0; 70 double h=fabs(cross(a,b,c))/dis(b,c); 71 if(dis(a,b)>r-eps&&dis(a,c)>r-eps)//两个端点都在圆的外面则分为两种情况 72 { 73 double angle=acos(dot(a,b,c)/dis(a,b)/dis(a,c)); 74 if(h>r-eps) 75 { 76 return 0.5*r*r*angle; 77 } 78 else if(dot(b,a,c)>0&&dot(c,a,b)>0) 79 { 80 double angle1=2*acos(h/r); 81 return 0.5*r*r*fabs(angle-angle1)+0.5*r*r*sin(angle1); 82 } 83 else 84 { 85 return 0.5*r*r*angle; 86 } 87 } 88 else if(dis(a,b)<r+eps&&dis(a,c)<r+eps)//两个端点都在圆内的情况 89 { 90 return 0.5*fabs(cross(a,b,c)); 91 } 92 else//一个端点在圆上一个端点在圆内的情况 93 { 94 if(dis(a,b)>dis(a,c))//默认b在圆内 95 { 96 swap(b,c); 97 } 98 if(fabs(dis(a,b))<eps)//ab距离为0直接返回0 99 { 100 return 0.0; 101 } 102 if(dot(b,a,c)<eps) 103 { 104 double angle1=acos(h/dis(a,b)); 105 double angle2=acos(h/r)-angle1; 106 double angle3=acos(h/dis(a,c))-acos(h/r); 107 return 0.5*dis(a,b)*r*sin(angle2)+0.5*r*r*angle3; 108 109 } 110 else 111 { 112 double angle1=acos(h/dis(a,b)); 113 double angle2=acos(h/r); 114 double angle3=acos(h/dis(a,c))-angle2; 115 return 0.5*r*dis(a,b)*sin(angle1+angle2)+0.5*r*r*angle3; 116 } 117 } 118 } 119 int n; 120 double x,y,h,x1,yy,R; 121 double gets(node O)//求圆与多边形面积并 122 { 123 for (int i=1;i<=n+1;i++) p[i]=a[i]; 124 for(int i=1;i<=n+1;i++) p[i]=p[i]-O; 125 O=node(0,0); 126 double sum=0; 127 for(int i=1;i<=n;i++) 128 { 129 int j=i+1; 130 double s=area(p[i],p[j],R); 131 if(cross(O,p[i],p[j])>0) sum+=s; else sum-=s; 132 } 133 return fabs(sum); 134 } 135 int dcmp(double x) 136 { 137 if(x > eps) return 1; 138 return x < -eps ? -1 : 0; 139 } 140 double calc(double x) 141 { 142 double l=1000000,r=-100000,m1,m2,f1=0,f2=0; 143 for (int i=1;i<=n;i++) 144 { 145 if (dcmp((a[i].x-x) * (a[i+1].x-x))<=0) 146 { 147 node tmp; 148 tmp.x=a[i].x+(a[i+1].x-a[i].x)/(a[i+1].x-a[i].x)*(x-a[i].x); 149 tmp.y=a[i].y+(a[i+1].y-a[i].y)/(a[i+1].x-a[i].x)*(x-a[i].x); 150 l=min(l,tmp.y); 151 r=max(r,tmp.y); 152 } 153 } 154 while (l+eps<r) 155 { 156 m1=(l+l+r)/3.0; 157 node bom=node(x,m1); 158 f1=gets(bom); 159 m2=(l+r+r)/3.0; 160 bom=node(x,m2); 161 f2=gets(bom); 162 if (f1>f2) r=m2;else l=m1; 163 } 164 return f1; 165 } 166 int main() 167 { 168 scanf("%d%lf",&n,&R); 169 double l=1000,r=-1000,m1,m2,f1,f2; 170 for(int i=1;i<=n;i++) 171 { 172 scanf("%lf%lf",&a[i].x,&a[i].y); 173 l=min(l,a[i].x); 174 r=max(r,a[i].x); 175 } 176 a[n+1]=a[1]; 177 while (l+eps<r) 178 { 179 m1=(l+l+r)/3.0;f1=calc(m1); 180 m2=(l+r+r)/3.0;f2=calc(m2); 181 if (f1>f2) r=m2;else l=m1; 182 } 183 printf("%.6lf ",f1); 184 return 0; 185 }