HDU 3982 半平面交+圆和凸多边形面积并

从这里看到的这个题。。很容易想到正解。

http://blog.csdn.net/zxy_snow/article/details/6739561

 

思路大概一样，就是不知道为什么我被卡精度了，，，

acos精度本来就不好，然后题目还要求输出百分比+五位小数，直接把精度卡了。

反正我写出来的当半径很大的时候误差会非常大，会达到3%左右。哎，查了一个下午，用几何画板模拟。真是恶心死了！

我早知道卡精度就不做了。

不知道“小媛”姐姐是怎么做的，实在不想看这个题了。

 

感觉我的思维和代码都是挺清晰的。。

尽管wa了，还是贴出来吧。。

 

View Code 
  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cstdlib>
  4 #include <algorithm>
  5 #include <cmath>
  6 #include <iostream>
  7 
  8 #define N 222222
  9 #define SIDE 11111
 10 #define EPS 1e-6
 11 #define BUG system("pause")
 12 
 13 const double PI=acos(-1.0);
 14 
 15 using namespace std;
 16 
 17 struct PO
 18 {
 19     double x,y,ag;
 20     bool fg;
 21     void prt() {printf("%lf     %lf    %d\n",x,y,fg);}
 22 }p[N],o,my,dp[N];
 23 
 24 struct LI
 25 {
 26     PO a,b; double ag;
 27     void in(double x1,double y1,double x2,double y2)
 28     {a.x=x1; a.y=y1; b.x=x2; b.y=y2;}
 29     void prt() {printf("%lf     %lf     %lf     %lf     %lf\n",a.x,a.y,b.x,b.y,ag);}
 30 }li[N],deq[N],qs;
 31 
 32 double r,ah;
 33 int n,cnt,tn,m;
 34 
 35 inline int dc(double x)
 36 {
 37     if(x>EPS) return 1;
 38     else if(x<-EPS) return -1;
 39     return 0;
 40 }
 41 
 42 inline PO operator -(PO a,PO b)
 43 {
 44     PO c; c.x=a.x-b.x; c.y=a.y-b.y;
 45     return c;
 46 }
 47 
 48 inline PO operator +(PO a,PO b)
 49 {
 50     PO c; c.x=a.x+b.x; c.y=a.y+b.y;
 51     return c;
 52 }
 53 
 54 inline PO operator *(PO a,double k)
 55 {
 56     a.x*=k; a.y*=k;
 57     return a;
 58 }
 59 
 60 inline PO operator /(PO a,double k)
 61 {
 62     a.x/=k; a.y/=k;
 63     return a;
 64 }
 65 
 66 inline double cross(PO a,PO b,PO c)
 67 {
 68     return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
 69 }
 70 
 71 inline double dot(PO a,PO b,PO c)
 72 {
 73     return (b.x-a.x)*(c.x-a.x)+(b.y-a.y)*(c.y-a.y);
 74 }
 75 
 76 inline double getlen(const LI &a)
 77 {
 78     return sqrt((a.b.x-a.a.x)*(a.b.x-a.a.x)+(a.b.y-a.a.y)*(a.b.y-a.a.y));
 79 }
 80 
 81 inline double getlen(const PO &a)
 82 {
 83     return sqrt(a.x*a.x+a.y*a.y);
 84 }
 85 
 86 inline double getdis(const PO &a,const PO &b)
 87 {
 88     return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
 89 }
 90 
 91 inline double getdisl(const PO &a,const LI &b)
 92 {
 93     PO t1=a-b.a,t2=a-b.b; 
 94     return fabs(cross(o,t1,t2))/getlen(b);
 95 }
 96 
 97 inline void read()
 98 {
 99     scanf("%lf%d",&r,&tn);
100     double a,b,c,d;
101     n=0;
102     for(int i=1;i<=tn;i++)
103     {
104         ++n;
105         scanf("%lf%lf%lf%lf",&li[i].a.x,&li[i].a.y,&li[i].b.x,&li[i].b.y);
106         if(dc(getdisl(o,li[i])-r)>0) n--; 
107     }
108     scanf("%lf%lf",&my.x,&my.y);
109 }
110 
111 inline void getdir()//规范每条线的方向 
112 {
113     for(int i=1;i<=n;i++)
114         if(dc(cross(li[i].a,li[i].b,my))<0) swap(li[i].a,li[i].b);
115 }
116 
117 inline void getangle()//计算所有直线的极角 
118 {
119     for(int i=1;i<=n;i++)
120         li[i].ag=atan2(li[i].b.y-li[i].a.y,li[i].b.x-li[i].a.x);
121     li[n+1].in(-SIDE,-SIDE,SIDE,-SIDE);
122     li[n+2].in(SIDE,-SIDE,SIDE,SIDE);
123     li[n+3].in(SIDE,SIDE,-SIDE,SIDE);
124     li[n+4].in(-SIDE,SIDE,-SIDE,-SIDE);
125     for(int i=n+1;i<=n+4;i++) li[i].ag=atan2(li[i].b.y-li[i].a.y,li[i].b.x-li[i].a.x);
126     n+=4;
127 }
128 
129 inline bool cmp(const LI &a,const LI &b)//极角排序，内侧在前 
130 {
131     if(dc(a.ag-b.ag)==0) return dc(cross(a.a,a.b,b.a))<0;
132     return a.ag>b.ag;
133 }
134 
135 inline PO getpoint(const LI &a,const LI &b)//直线交点 
136 {
137     double k1=cross(a.a,b.b,b.a),k2=cross(a.b,b.a,b.b);
138     PO c;
139     c=a.a+(a.b-a.a)*k1/(k1+k2);
140     return c;
141 }
142 
143 inline void getcut()//半平面交 
144 {
145     sort(li+1,li+1+n,cmp);
146     m=1;
147     for(int i=2;i<=n;i++)
148         if(dc(li[i].ag-li[m].ag)!=0)
149             li[++m]=li[i];
150     deq[1]=li[1]; deq[2]=li[2];
151     int bot=1,top=2;
152     for(int i=3;i<=m;i++)
153     {
154         while(bot<top&&dc(cross(li[i].a,li[i].b,getpoint(deq[top],deq[top-1])))<0) top--;
155         while(bot<top&&dc(cross(li[i].a,li[i].b,getpoint(deq[bot],deq[bot+1])))<0) bot++;
156         deq[++top]=li[i];
157     }
158     while(bot<top&&dc(cross(deq[bot].a,deq[bot].b,getpoint(deq[top],deq[top-1])))<0) top--;
159     while(bot<top&&dc(cross(deq[top].a,deq[top].b,getpoint(deq[bot],deq[bot+1])))<0) bot++;
160     if(bot==top) return;
161     cnt=0;
162     for(int i=bot;i<top;i++) p[++cnt]=getpoint(deq[i],deq[i+1]);
163     if(top-1>bot) p[++cnt]=getpoint(deq[bot],deq[top]);
164 }
165 
166 inline double getg(const PO &s,const PO &t)//弓形面积 
167 {
168     double res;
169     double cc=dot(o,s,t)/getlen(s)/getlen(t);//夹角的余弦值，不能用正弦 asin的范围是[-π/2, π/2] acos的范围是[0, π] 
170     cc=acos(cc);
171     if(dc(cross(o,s,t))<0) {cc=2*PI-cc;res=fabs(cross(o,s,t))/2;} //三角形 
172     else res=-fabs(cross(o,s,t)/2);
173     double la=r*cc; //弧长 
174     double ss=r*la/2;//扇形面积 
175     return res+ss;
176 }
177 
178 inline PO getf(const PO &a,const PO &b)//ab的右旋法向量 
179 {
180     PO zt=b-a,rt;
181     rt.x=zt.y; rt.y=-zt.x;
182     return rt;
183 }
184 
185 inline PO rotate(const PO &a,double sss,double ccc)
186 {
187     PO rt;
188     rt.x=a.x*ccc-a.y*sss;
189     rt.y=a.x*sss+a.y*ccc;
190     return rt;
191 }
192 
193 inline void getcpoint(const LI &a,PO &s,PO &t)//直线与圆的交点 
194 {
195     double h=getdisl(o,a);
196     double tt=sqrt(r*r-h*h);
197     double sss=tt/r,ccc=h/r;
198     PO f;
199     if(dc(cross(o,a.a,a.b))<0) f=getf(a.b,a.a);//判断顺逆时针旋转 
200     else f=getf(a.a,a.b);
201     f=f/getlen(a)*r;
202     s=rotate(f,-sss,ccc);
203     t=rotate(f,sss,ccc);
204 }
205 
206 inline bool onseg(const PO &a,const LI &b)//a在直线b上，求a是否在线段b上 
207 {
208     PO s1=b.a-a,s2=b.b-a;
209     if(dc(dot(o,s1,s2))<=0) return true;
210     return false;
211 }
212 
213 inline bool cmp1(const PO &a,const PO &b)
214 {
215     return a.ag<b.ag;
216 }
217 
218 inline bool cmp2(const PO &a,const PO &b)//unique的去重函数 
219 {
220     if(dc(a.ag-b.ag)==0) return true;
221     return false;
222 }
223 
224 inline void cir()
225 {
226     for(int i=1;i<=(cnt>>1);i++) swap(p[i],p[cnt-i+1]);//半平面交是顺时针的，转换为逆时针 
227     n=0;
228     for(int i=1;i<=cnt;i++)
229     {
230         int dis=dc(getdis(o,p[i])-r);
231         if(dis==-1)//凸多边形在圆内的顶点 
232         {
233             dp[++n]=p[i];
234             dp[n].fg=0;
235         }
236         else if(dis==0)//在圆上的顶点 
237         {
238             dp[++n]=p[i];
239             dp[n].fg=1;
240         }
241     }
242     p[cnt+1]=p[1];
243     PO s,t;LI pl;
244     for(int i=1;i<=cnt;i++)//凸多边形和圆的交点 
245     {
246         pl.a=p[i]; pl.b=p[i+1];
247         if(dc(getdisl(o,pl)-r)<0)//在圆内 
248         {
249             getcpoint(pl,s,t);
250             if(onseg(s,pl))
251             {
252                 dp[++n]=s;
253                 dp[n].fg=1;
254             }
255             if(onseg(t,pl))
256             {
257                 dp[++n]=t;
258                 dp[n].fg=1;
259             }
260         }
261     }
262     for(int i=1;i<=n;i++) dp[i].ag=atan2(dp[i].y-my.y,dp[i].x-my.x);
263     sort(dp+1,dp+1+n,cmp1); 
264     n=unique(dp+1,dp+1+n,cmp2)-dp-1;
265     dp[n+1]=dp[1];
266     ah=0;
267     for(int i=1;i<=n;i++)
268         if(dp[i].fg&&dp[i+1].fg) ah+=getg(dp[i],dp[i+1]);
269 }
270 
271 inline double getarea()
272 {
273     double res=0;
274     dp[n+1]=dp[1];
275     for(int i=1;i<=n;i++) res+=cross(o,dp[i],dp[i+1]);
276     return res/2;
277 }
278 
279 inline void go()
280 {
281     if(n==0)
282     {
283         printf("100.00000%%\n");
284         return;
285     }
286     getdir();
287     getangle();
288     getcut();
289     cir();
290     double area=fabs(getarea());
291     printf("%.5lf%%\n",(area+ah)/(PI*r*r)*100);
292     
293 }
294 
295 int main()
296 {
297     int cas; scanf("%d",&cas);
298     for(int i=1;i<=cas;i++)
299     {
300         printf("Case %d: ",i);
301         read(),go();
302     }
303     return 0;
304 }

（至今写的最长的计算几何题了。。。）

 

 

 

posted @ 2013-03-02 19:43  proverbs  阅读(589)  评论(0)  编辑  收藏  举报