1 #include<vector>   2 #include<list>   3 #include<map>   4 #include<set>   5 #include<deque>   6 #include<queue>   7 #include<stack>   8 #include<bitset>   9 #include<algorithm>  10 #include<functional>  11 #include<numeric>  12 #include<utility>  13 #include<iostream>  14 #include<sstream>  15 #include<iomanip>  16 #include<cstdio>  17 #include<cmath>  18 #include<cstdlib>  19 #include<cctype>  20 #include<string>  21 #include<cstring>  22 #include<cstdio>  23 #include<cmath>  24 #include<cstdlib>  25 #include<ctime>  26 #include<climits>  27 #include<complex>  28 #define mp make_pair  29 #define pb push_back  30 using namespace std;  31 const double eps=1e-8;  32 const double pi=acos(-1.0);  33 const double inf=1e20;  34 const int maxp=1111;  35 int dblcmp(double d)  36 {  37     if (fabs(d)<eps)return 0;  38     return d>eps?1:-1;  39 }  40 inline double sqr(double x){return x*x;}  41 struct point  42 {  43     double x,y;  44     point(){}  45     point(double _x,double _y):  46     x(_x),y(_y){};  47     void input()  48     {  49         scanf("%lf%lf",&x,&y);  50     }  51     void output()  52     {  53         printf("%.2f %.2f\n",x,y);  54     }  55     bool operator==(point a)const  56     {  57         return dblcmp(a.x-x)==0&&dblcmp(a.y-y)==0;  58     }  59     bool operator<(point a)const  60     {  61         return dblcmp(a.x-x)==0?dblcmp(y-a.y)<0:x<a.x;  62     }  63     double len()  64     {  65         return hypot(x,y);  66     }  67     double len2()  68     {  69         return x*x+y*y;  70     }  71     double distance(point p)  72     {  73         return hypot(x-p.x,y-p.y);  74     }  75     point add(point p)  76     {  77         return point(x+p.x,y+p.y);  78     }  79     point sub(point p)  80     {  81         return point(x-p.x,y-p.y);  82     }  83     point mul(double b)  84     {  85         return point(x*b,y*b);  86     }  87     point div(double b)  88     {  89         return point(x/b,y/b);  90     }  91     double dot(point p)  92     {  93         return x*p.x+y*p.y;  94     }  95     double det(point p)  96     {  97         return x*p.y-y*p.x;  98     }  99     double rad(point a,point b) 100     { 101         point p=*this; 102         return fabs(atan2(fabs(a.sub(p).det(b.sub(p))),a.sub(p).dot(b.sub(p)))); 103     } 104     point trunc(double r) 105     { 106         double l=len(); 107         if (!dblcmp(l))return *this; 108         r/=l; 109         return point(x*r,y*r); 110     } 111     point rotleft() 112     { 113         return point(-y,x); 114     } 115     point rotright() 116     { 117         return point(y,-x); 118     } 119     point rotate(point p,double angle)//绕点p逆时针旋转angle角度 120     { 121         point v=this->sub(p); 122         double c=cos(angle),s=sin(angle); 123         return point(p.x+v.x*c-v.y*s,p.y+v.x*s+v.y*c); 124     } 125 }; 126 struct line 127 { 128     point a,b; 129     line(){} 130     line(point _a,point _b) 131     { 132         a=_a; 133         b=_b; 134     } 135     bool operator==(line v) 136     { 137         return (a==v.a)&&(b==v.b); 138     } 139     //倾斜角angle 140     line(point p,double angle) 141     { 142         a=p; 143         if (dblcmp(angle-pi/2)==0) 144         { 145             b=a.add(point(0,1)); 146         } 147         else 148         { 149             b=a.add(point(1,tan(angle))); 150         } 151     } 152     //ax+by+c=0 153     line(double _a,double _b,double _c) 154     { 155         if (dblcmp(_a)==0) 156         { 157             a=point(0,-_c/_b); 158             b=point(1,-_c/_b); 159         } 160         else if (dblcmp(_b)==0) 161         { 162             a=point(-_c/_a,0); 163             b=point(-_c/_a,1); 164         } 165         else 166         { 167             a=point(0,-_c/_b); 168             b=point(1,(-_c-_a)/_b); 169         } 170     } 171     void input() 172     { 173         a.input(); 174         b.input(); 175     } 176     void adjust() 177     { 178         if (b<a)swap(a,b); 179     } 180     double length() 181     { 182         return a.distance(b); 183     } 184     double angle()//直线倾斜角 0<=angle<180 185     { 186         double k=atan2(b.y-a.y,b.x-a.x); 187         if (dblcmp(k)<0)k+=pi; 188         if (dblcmp(k-pi)==0)k-=pi; 189         return k; 190     } 191     //点和线段关系 192     //1 在逆时针 193     //2 在顺时针 194     //3 平行 195     int relation(point p) 196     { 197         int c=dblcmp(p.sub(a).det(b.sub(a))); 198         if (c<0)return 1; 199         if (c>0)return 2; 200         return 3; 201     } 202     bool pointonseg(point p) 203     { 204         return dblcmp(p.sub(a).det(b.sub(a)))==0&&dblcmp(p.sub(a).dot(p.sub(b)))<=0; 205     } 206     bool parallel(line v) 207     { 208         return dblcmp(b.sub(a).det(v.b.sub(v.a)))==0; 209     } 210     //2 规范相交 211     //1 非规范相交 212     //0 不相交 213     int segcrossseg(line v) 214     { 215         int d1=dblcmp(b.sub(a).det(v.a.sub(a))); 216         int d2=dblcmp(b.sub(a).det(v.b.sub(a))); 217         int d3=dblcmp(v.b.sub(v.a).det(a.sub(v.a))); 218         int d4=dblcmp(v.b.sub(v.a).det(b.sub(v.a))); 219         if ((d1^d2)==-2&&(d3^d4)==-2)return 2; 220         return (d1==0&&dblcmp(v.a.sub(a).dot(v.a.sub(b)))<=0|| 221                 d2==0&&dblcmp(v.b.sub(a).dot(v.b.sub(b)))<=0|| 222                 d3==0&&dblcmp(a.sub(v.a).dot(a.sub(v.b)))<=0|| 223                 d4==0&&dblcmp(b.sub(v.a).dot(b.sub(v.b)))<=0); 224     } 225     int linecrossseg(line v)//*this seg v line 226     { 227         int d1=dblcmp(b.sub(a).det(v.a.sub(a))); 228         int d2=dblcmp(b.sub(a).det(v.b.sub(a))); 229         if ((d1^d2)==-2)return 2; 230         return (d1==0||d2==0); 231     } 232     //0 平行 233     //1 重合 234     //2 相交 235     int linecrossline(line v) 236     { 237         if ((*this).parallel(v)) 238         { 239             return v.relation(a)==3; 240         } 241         return 2; 242     } 243     point crosspoint(line v) 244     { 245         double a1=v.b.sub(v.a).det(a.sub(v.a)); 246         double a2=v.b.sub(v.a).det(b.sub(v.a)); 247         return point((a.x*a2-b.x*a1)/(a2-a1),(a.y*a2-b.y*a1)/(a2-a1)); 248     } 249     double dispointtoline(point p) 250     { 251         return fabs(p.sub(a).det(b.sub(a)))/length(); 252     } 253     double dispointtoseg(point p) 254     { 255         if (dblcmp(p.sub(b).dot(a.sub(b)))<0||dblcmp(p.sub(a).dot(b.sub(a)))<0) 256         { 257             return min(p.distance(a),p.distance(b)); 258         } 259         return dispointtoline(p); 260     } 261     point lineprog(point p) 262     { 263         return a.add(b.sub(a).mul(b.sub(a).dot(p.sub(a))/b.sub(a).len2())); 264     } 265     point symmetrypoint(point p) 266     { 267         point q=lineprog(p); 268         return point(2*q.x-p.x,2*q.y-p.y); 269     } 270 }; 271 struct circle 272 { 273     point p; 274     double r; 275     circle(){} 276     circle(point _p,double _r): 277     p(_p),r(_r){}; 278     circle(double x,double y,double _r): 279     p(point(x,y)),r(_r){}; 280     circle(point a,point b,point c)//三角形的外接圆 281     { 282         p=line(a.add(b).div(2),a.add(b).div(2).add(b.sub(a).rotleft())).crosspoint(line(c.add(b).div(2),c.add(b).div(2).add(b.sub(c).rotleft()))); 283         r=p.distance(a); 284     } 285     circle(point a,point b,point c,bool t)//三角形的内切圆 286     { 287         line u,v; 288         double m=atan2(b.y-a.y,b.x-a.x),n=atan2(c.y-a.y,c.x-a.x); 289         u.a=a; 290         u.b=u.a.add(point(cos((n+m)/2),sin((n+m)/2))); 291         v.a=b; 292         m=atan2(a.y-b.y,a.x-b.x),n=atan2(c.y-b.y,c.x-b.x); 293         v.b=v.a.add(point(cos((n+m)/2),sin((n+m)/2))); 294         p=u.crosspoint(v); 295         r=line(a,b).dispointtoseg(p); 296     } 297     void input() 298     { 299         p.input(); 300         scanf("%lf",&r); 301     } 302     void output() 303     { 304         printf("%.2lf %.2lf %.2lf\n",p.x,p.y,r); 305     } 306     bool operator==(circle v) 307     { 308         return ((p==v.p)&&dblcmp(r-v.r)==0); 309     } 310     bool operator<(circle v)const 311     { 312         return ((p<v.p)||(p==v.p)&&dblcmp(r-v.r)<0); 313     } 314     double area() 315     { 316         return pi*sqr(r); 317     } 318     double circumference() 319     { 320         return 2*pi*r; 321     } 322     //0 圆外 323     //1 圆上 324     //2 圆内 325     int relation(point b) 326     { 327         double dst=b.distance(p); 328         if (dblcmp(dst-r)<0)return 2; 329         if (dblcmp(dst-r)==0)return 1; 330         return 0; 331     } 332     int relationseg(line v) 333     { 334         double dst=v.dispointtoseg(p); 335         if (dblcmp(dst-r)<0)return 2; 336         if (dblcmp(dst-r)==0)return 1; 337         return 0; 338     } 339     int relationline(line v) 340     { 341         double dst=v.dispointtoline(p); 342         if (dblcmp(dst-r)<0)return 2; 343         if (dblcmp(dst-r)==0)return 1; 344         return 0; 345     } 346     //过a b两点 半径r的两个圆 347     int getcircle(point a,point b,double r,circle&c1,circle&c2) 348     { 349         circle x(a,r),y(b,r); 350         int t=x.pointcrosscircle(y,c1.p,c2.p); 351         if (!t)return 0; 352         c1.r=c2.r=r; 353         return t; 354     } 355     //与直线u相切 过点q 半径r1的圆 356     int getcircle(line u,point q,double r1,circle &c1,circle &c2) 357     { 358         double dis=u.dispointtoline(q); 359         if (dblcmp(dis-r1*2)>0)return 0; 360         if (dblcmp(dis)==0) 361         { 362             c1.p=q.add(u.b.sub(u.a).rotleft().trunc(r1)); 363             c2.p=q.add(u.b.sub(u.a).rotright().trunc(r1)); 364             c1.r=c2.r=r1; 365             return 2; 366         } 367         line u1=line(u.a.add(u.b.sub(u.a).rotleft().trunc(r1)),u.b.add(u.b.sub(u.a).rotleft().trunc(r1))); 368         line u2=line(u.a.add(u.b.sub(u.a).rotright().trunc(r1)),u.b.add(u.b.sub(u.a).rotright().trunc(r1))); 369         circle cc=circle(q,r1); 370         point p1,p2; 371         if (!cc.pointcrossline(u1,p1,p2))cc.pointcrossline(u2,p1,p2); 372         c1=circle(p1,r1); 373         if (p1==p2) 374         { 375             c2=c1;return 1; 376         } 377         c2=circle(p2,r1); 378         return 2; 379     } 380     //同时与直线u,v相切 半径r1的圆 381     int getcircle(line u,line v,double r1,circle &c1,circle &c2,circle &c3,circle &c4) 382     { 383         if (u.parallel(v))return 0; 384         line u1=line(u.a.add(u.b.sub(u.a).rotleft().trunc(r1)),u.b.add(u.b.sub(u.a).rotleft().trunc(r1))); 385         line u2=line(u.a.add(u.b.sub(u.a).rotright().trunc(r1)),u.b.add(u.b.sub(u.a).rotright().trunc(r1))); 386         line v1=line(v.a.add(v.b.sub(v.a).rotleft().trunc(r1)),v.b.add(v.b.sub(v.a).rotleft().trunc(r1))); 387         line v2=line(v.a.add(v.b.sub(v.a).rotright().trunc(r1)),v.b.add(v.b.sub(v.a).rotright().trunc(r1))); 388         c1.r=c2.r=c3.r=c4.r=r1; 389         c1.p=u1.crosspoint(v1); 390         c2.p=u1.crosspoint(v2); 391         c3.p=u2.crosspoint(v1); 392         c4.p=u2.crosspoint(v2); 393         return 4; 394     } 395     //同时与不相交圆cx,cy相切 半径为r1的圆 396     int getcircle(circle cx,circle cy,double r1,circle&c1,circle&c2) 397     { 398         circle x(cx.p,r1+cx.r),y(cy.p,r1+cy.r); 399         int t=x.pointcrosscircle(y,c1.p,c2.p); 400         if (!t)return 0; 401         c1.r=c2.r=r1; 402         return t; 403     } 404     int pointcrossline(line v,point &p1,point &p2)//求与线段交要先判断relationseg 405     { 406         if (!(*this).relationline(v))return 0; 407         point a=v.lineprog(p); 408         double d=v.dispointtoline(p); 409         d=sqrt(r*r-d*d); 410         if (dblcmp(d)==0) 411         { 412             p1=a; 413             p2=a; 414             return 1; 415         } 416         p1=a.sub(v.b.sub(v.a).trunc(d)); 417         p2=a.add(v.b.sub(v.a).trunc(d)); 418         return 2; 419     } 420     //5 相离 421     //4 外切 422     //3 相交 423     //2 内切 424     //1 内含 425     int relationcircle(circle v) 426     { 427         double d=p.distance(v.p); 428         if (dblcmp(d-r-v.r)>0)return 5; 429         if (dblcmp(d-r-v.r)==0)return 4; 430         double l=fabs(r-v.r); 431         if (dblcmp(d-r-v.r)<0&&dblcmp(d-l)>0)return 3; 432         if (dblcmp(d-l)==0)return 2; 433         if (dblcmp(d-l)<0)return 1; 434     } 435     int pointcrosscircle(circle v,point &p1,point &p2) 436     { 437         int rel=relationcircle(v); 438         if (rel==1||rel==5)return 0; 439         double d=p.distance(v.p); 440         double l=(d+(sqr(r)-sqr(v.r))/d)/2; 441         double h=sqrt(sqr(r)-sqr(l)); 442         p1=p.add(v.p.sub(p).trunc(l).add(v.p.sub(p).rotleft().trunc(h))); 443         p2=p.add(v.p.sub(p).trunc(l).add(v.p.sub(p).rotright().trunc(h))); 444         if (rel==2||rel==4) 445         { 446             return 1; 447         } 448         return 2; 449     } 450     //过一点做圆的切线 (先判断点和圆关系) 451     int tangentline(point q,line &u,line &v) 452     { 453         int x=relation(q); 454         if (x==2)return 0; 455         if (x==1) 456         { 457             u=line(q,q.add(q.sub(p).rotleft())); 458             v=u; 459             return 1; 460         } 461         double d=p.distance(q); 462         double l=sqr(r)/d; 463         double h=sqrt(sqr(r)-sqr(l)); 464         u=line(q,p.add(q.sub(p).trunc(l).add(q.sub(p).rotleft().trunc(h)))); 465         v=line(q,p.add(q.sub(p).trunc(l).add(q.sub(p).rotright().trunc(h)))); 466         return 2; 467     } 468     double areacircle(circle v) 469     { 470         int rel=relationcircle(v); 471         if (rel>=4)return 0.0; 472         if (rel<=2)return min(area(),v.area()); 473         double d=p.distance(v.p); 474         double hf=(r+v.r+d)/2.0; 475         double ss=2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d)); 476         double a1=acos((r*r+d*d-v.r*v.r)/(2.0*r*d)); 477         a1=a1*r*r; 478         double a2=acos((v.r*v.r+d*d-r*r)/(2.0*v.r*d)); 479         a2=a2*v.r*v.r; 480         return a1+a2-ss; 481     } 482     double areatriangle(point a,point b) 483     { 484         if (dblcmp(p.sub(a).det(p.sub(b))==0))return 0.0; 485         point q[5]; 486         int len=0; 487         q[len++]=a; 488         line l(a,b); 489         point p1,p2; 490         if (pointcrossline(l,q[1],q[2])==2) 491         { 492             if (dblcmp(a.sub(q[1]).dot(b.sub(q[1])))<0)q[len++]=q[1]; 493             if (dblcmp(a.sub(q[2]).dot(b.sub(q[2])))<0)q[len++]=q[2]; 494         } 495         q[len++]=b; 496         if (len==4&&(dblcmp(q[0].sub(q[1]).dot(q[2].sub(q[1])))>0))swap(q[1],q[2]); 497         double res=0; 498         int i; 499         for (i=0;i<len-1;i++) 500         { 501             if (relation(q[i])==0||relation(q[i+1])==0) 502             { 503                 double arg=p.rad(q[i],q[i+1]); 504                 res+=r*r*arg/2.0; 505             } 506             else 507             { 508                 res+=fabs(q[i].sub(p).det(q[i+1].sub(p))/2.0); 509             } 510         } 511         return res; 512     } 513 }; 514 struct polygon 515 { 516     int n; 517     point p[maxp]; 518     line l[maxp]; 519     void input() 520     { 521         n=4; 522         for (int i=0;i<n;i++) 523         { 524             p[i].input(); 525         } 526     } 527     void add(point q) 528     { 529         p[n++]=q; 530     } 531     void getline() 532     { 533         for (int i=0;i<n;i++) 534         { 535             l[i]=line(p[i],p[(i+1)%n]); 536         } 537     } 538     struct cmp 539     { 540         point p; 541         cmp(const point &p0){p=p0;} 542         bool operator()(const point &aa,const point &bb) 543         { 544             point a=aa,b=bb; 545             int d=dblcmp(a.sub(p).det(b.sub(p))); 546             if (d==0) 547             { 548                 return dblcmp(a.distance(p)-b.distance(p))<0; 549             } 550             return d>0; 551         } 552     }; 553     void norm() 554     { 555         point mi=p[0]; 556         for (int i=1;i<n;i++)mi=min(mi,p[i]); 557         sort(p,p+n,cmp(mi)); 558     } 559     void getconvex(polygon &convex) 560     { 561         int i,j,k; 562         sort(p,p+n); 563         convex.n=n; 564         for (i=0;i<min(n,2);i++) 565         { 566             convex.p[i]=p[i]; 567         } 568         if (n<=2)return; 569         int &top=convex.n; 570         top=1; 571         for (i=2;i<n;i++) 572         { 573             while (top&&convex.p[top].sub(p[i]).det(convex.p[top-1].sub(p[i]))<=0) 574                 top--; 575             convex.p[++top]=p[i]; 576         } 577         int temp=top; 578         convex.p[++top]=p[n-2]; 579         for (i=n-3;i>=0;i--) 580         { 581             while (top!=temp&&convex.p[top].sub(p[i]).det(convex.p[top-1].sub(p[i]))<=0) 582                 top--; 583             convex.p[++top]=p[i]; 584         } 585     } 586     bool isconvex() 587     { 588         bool s[3]; 589         memset(s,0,sizeof(s)); 590         int i,j,k; 591         for (i=0;i<n;i++) 592         { 593             j=(i+1)%n; 594             k=(j+1)%n; 595             s[dblcmp(p[j].sub(p[i]).det(p[k].sub(p[i])))+1]=1; 596             if (s[0]&&s[2])return 0; 597         } 598         return 1; 599     } 600     //3 点上 601     //2 边上 602     //1 内部 603     //0 外部 604     int relationpoint(point q) 605     { 606         int i,j; 607         for (i=0;i<n;i++) 608         { 609             if (p[i]==q)return 3; 610         } 611         getline(); 612         for (i=0;i<n;i++) 613         { 614             if (l[i].pointonseg(q))return 2; 615         } 616         int cnt=0; 617         for (i=0;i<n;i++) 618         { 619             j=(i+1)%n; 620             int k=dblcmp(q.sub(p[j]).det(p[i].sub(p[j]))); 621             int u=dblcmp(p[i].y-q.y); 622             int v=dblcmp(p[j].y-q.y); 623             if (k>0&&u<0&&v>=0)cnt++; 624             if (k<0&&v<0&&u>=0)cnt--; 625         } 626         return cnt!=0; 627     } 628     //1 在多边形内长度为正 629     //2 相交或与边平行 630     //0 无任何交点 631     int relationline(line u) 632     { 633         int i,j,k=0; 634         getline(); 635         for (i=0;i<n;i++) 636         { 637             if (l[i].segcrossseg(u)==2)return 1; 638             if (l[i].segcrossseg(u)==1)k=1; 639         } 640         if (!k)return 0; 641         vector<point>vp; 642         for (i=0;i<n;i++) 643         { 644             if (l[i].segcrossseg(u)) 645             { 646                 if (l[i].parallel(u)) 647                 { 648                     vp.pb(u.a); 649                     vp.pb(u.b); 650                     vp.pb(l[i].a); 651                     vp.pb(l[i].b); 652                     continue; 653                 } 654                 vp.pb(l[i].crosspoint(u)); 655             } 656         } 657         sort(vp.begin(),vp.end()); 658         int sz=vp.size(); 659         for (i=0;i<sz-1;i++) 660         { 661             point mid=vp[i].add(vp[i+1]).div(2); 662             if (relationpoint(mid)==1)return 1; 663         } 664         return 2; 665     } 666     //直线u切割凸多边形左侧 667     //注意直线方向 668     void convexcut(line u,polygon &po) 669     { 670         int i,j,k; 671         int &top=po.n; 672         top=0; 673         for (i=0;i<n;i++) 674         { 675             int d1=dblcmp(p[i].sub(u.a).det(u.b.sub(u.a))); 676             int d2=dblcmp(p[(i+1)%n].sub(u.a).det(u.b.sub(u.a))); 677             if (d1>=0)po.p[top++]=p[i]; 678             if (d1*d2<0)po.p[top++]=u.crosspoint(line(p[i],p[(i+1)%n])); 679         } 680     } 681     double getcircumference() 682     { 683         double sum=0; 684         int i; 685         for (i=0;i<n;i++) 686         { 687             sum+=p[i].distance(p[(i+1)%n]); 688         } 689         return sum; 690     } 691     double getarea() 692     { 693         double sum=0; 694         int i; 695         for (i=0;i<n;i++) 696         { 697             sum+=p[i].det(p[(i+1)%n]); 698         } 699         return fabs(sum)/2; 700     } 701     bool getdir()//1代表逆时针 0代表顺时针 702     { 703         double sum=0; 704         int i; 705         for (i=0;i<n;i++) 706         { 707             sum+=p[i].det(p[(i+1)%n]); 708         } 709         if (dblcmp(sum)>0)return 1; 710         return 0; 711     } 712     point getbarycentre() 713     { 714         point ret(0,0); 715         double area=0; 716         int i; 717         for (i=1;i<n-1;i++) 718         { 719             double tmp=p[i].sub(p[0]).det(p[i+1].sub(p[0])); 720             if (dblcmp(tmp)==0)continue; 721             area+=tmp; 722             ret.x+=(p[0].x+p[i].x+p[i+1].x)/3*tmp; 723             ret.y+=(p[0].y+p[i].y+p[i+1].y)/3*tmp; 724         } 725         if (dblcmp(area))ret=ret.div(area); 726         return ret; 727     } 728     double areaintersection(polygon po) 729     { 730     } 731     double areaunion(polygon po) 732     { 733         return getarea()+po.getarea()-areaintersection(po); 734     } 735     double areacircle(circle c) 736     { 737         int i,j,k,l,m; 738         double ans=0; 739         for (i=0;i<n;i++) 740         { 741             int j=(i+1)%n; 742             if (dblcmp(p[j].sub(c.p).det(p[i].sub(c.p)))>=0) 743             { 744                 ans+=c.areatriangle(p[i],p[j]); 745             } 746             else 747             { 748                 ans-=c.areatriangle(p[i],p[j]); 749             } 750         } 751         return fabs(ans); 752     } 753     //多边形和圆关系 754     //0 一部分在圆外 755     //1 与圆某条边相切 756     //2 完全在圆内 757     int relationcircle(circle c) 758     { 759         getline(); 760         int i,x=2; 761         if (relationpoint(c.p)!=1)return 0; 762         for (i=0;i<n;i++) 763         { 764             if (c.relationseg(l[i])==2)return 0; 765             if (c.relationseg(l[i])==1)x=1; 766         } 767         return x; 768     } 769     void find(int st,point tri[],circle &c) 770     { 771         if (!st) 772         { 773             c=circle(point(0,0),-2); 774         } 775         if (st==1) 776         { 777             c=circle(tri[0],0); 778         } 779         if (st==2) 780         { 781             c=circle(tri[0].add(tri[1]).div(2),tri[0].distance(tri[1])/2.0); 782         } 783         if (st==3) 784         { 785             c=circle(tri[0],tri[1],tri[2]); 786         } 787     } 788     void solve(int cur,int st,point tri[],circle &c) 789     { 790         find(st,tri,c); 791         if (st==3)return; 792         int i; 793         for (i=0;i<cur;i++) 794         { 795             if (dblcmp(p[i].distance(c.p)-c.r)>0) 796             { 797                 tri[st]=p[i]; 798                 solve(i,st+1,tri,c); 799             } 800         } 801     } 802     circle mincircle()//点集最小圆覆盖 803     { 804         random_shuffle(p,p+n); 805         point tri[4]; 806         circle c; 807         solve(n,0,tri,c); 808         return c; 809     } 810     int circlecover(double r)//单位圆覆盖 811     { 812         int ans=0,i,j; 813         vector<pair<double,int> >v; 814         for (i=0;i<n;i++) 815         { 816             v.clear(); 817             for (j=0;j<n;j++)if (i!=j) 818             { 819                 point q=p[i].sub(p[j]); 820                 double d=q.len(); 821                 if (dblcmp(d-2*r)<=0) 822                 { 823                     double arg=atan2(q.y,q.x); 824                     if (dblcmp(arg)<0)arg+=2*pi; 825                     double t=acos(d/(2*r)); 826                     v.push_back(make_pair(arg-t+2*pi,-1)); 827                     v.push_back(make_pair(arg+t+2*pi,1)); 828                 } 829             } 830             sort(v.begin(),v.end()); 831             int cur=0; 832             for (j=0;j<v.size();j++) 833             { 834                 if (v[j].second==-1)++cur; 835                 else --cur; 836                 ans=max(ans,cur); 837             } 838         } 839         return ans+1; 840     } 841     int pointinpolygon(point q)//点在凸多边形内部的判定 842     { 843         if (getdir())reverse(p,p+n); 844         if (dblcmp(q.sub(p[0]).det(p[n-1].sub(p[0])))==0) 845         { 846             if (line(p[n-1],p[0]).pointonseg(q))return n-1; 847             return -1; 848         } 849         int low=1,high=n-2,mid; 850         while (low<=high) 851         { 852             mid=(low+high)>>1; 853             if (dblcmp(q.sub(p[0]).det(p[mid].sub(p[0])))>=0&&dblcmp(q.sub(p[0]).det(p[mid+1].sub(p[0])))<0) 854             { 855                 polygon c; 856                 c.p[0]=p[mid]; 857                 c.p[1]=p[mid+1]; 858                 c.p[2]=p[0]; 859                 c.n=3; 860                 if (c.relationpoint(q))return mid; 861                 return -1; 862             } 863             if (dblcmp(q.sub(p[0]).det(p[mid].sub(p[0])))>0) 864             { 865                 low=mid+1; 866             } 867             else 868             { 869                 high=mid-1; 870             } 871         } 872         return -1; 873     } 874 }; 875 struct polygons 876 { 877     vector<polygon>p; 878     polygons() 879     { 880         p.clear(); 881     } 882     void clear() 883     { 884         p.clear(); 885     } 886     void push(polygon q) 887     { 888         if (dblcmp(q.getarea()))p.pb(q); 889     } 890     vector<pair<double,int> >e; 891     void ins(point s,point t,point X,int i) 892     { 893         double r=fabs(t.x-s.x)>eps?(X.x-s.x)/(t.x-s.x):(X.y-s.y)/(t.y-s.y); 894         r=min(r,1.0);r=max(r,0.0); 895         e.pb(mp(r,i)); 896     } 897     double polyareaunion() 898     { 899         double ans=0.0; 900         int c0,c1,c2,i,j,k,w; 901         for (i=0;i<p.size();i++) 902         { 903             if (p[i].getdir()==0)reverse(p[i].p,p[i].p+p[i].n); 904         } 905         for (i=0;i<p.size();i++) 906         { 907             for (k=0;k<p[i].n;k++) 908             { 909                 point &s=p[i].p[k],&t=p[i].p[(k+1)%p[i].n]; 910                 if (!dblcmp(s.det(t)))continue; 911                 e.clear(); 912                 e.pb(mp(0.0,1)); 913                 e.pb(mp(1.0,-1)); 914                 for (j=0;j<p.size();j++)if (i!=j) 915                 { 916                     for (w=0;w<p[j].n;w++) 917                     { 918                         point a=p[j].p[w],b=p[j].p[(w+1)%p[j].n],c=p[j].p[(w-1+p[j].n)%p[j].n]; 919                         c0=dblcmp(t.sub(s).det(c.sub(s))); 920                         c1=dblcmp(t.sub(s).det(a.sub(s))); 921                         c2=dblcmp(t.sub(s).det(b.sub(s))); 922                         if (c1*c2<0)ins(s,t,line(s,t).crosspoint(line(a,b)),-c2); 923                         else if (!c1&&c0*c2<0)ins(s,t,a,-c2); 924                         else if (!c1&&!c2) 925                         { 926                             int c3=dblcmp(t.sub(s).det(p[j].p[(w+2)%p[j].n].sub(s))); 927                             int dp=dblcmp(t.sub(s).dot(b.sub(a))); 928                             if (dp&&c0)ins(s,t,a,dp>0?c0*((j>i)^(c0<0)):-(c0<0)); 929                             if (dp&&c3)ins(s,t,b,dp>0?-c3*((j>i)^(c3<0)):c3<0); 930                         } 931                     } 932                 } 933                 sort(e.begin(),e.end()); 934                 int ct=0; 935                 double tot=0.0,last; 936                 for (j=0;j<e.size();j++) 937                 { 938                     if (ct==2)tot+=e[j].first-last; 939                     ct+=e[j].second; 940                     last=e[j].first; 941                 } 942                 ans+=s.det(t)*tot; 943             } 944         } 945         return fabs(ans)*0.5; 946     } 947 }; 948 const int maxn=500; 949 struct circles 950 { 951     circle c[maxn]; 952     double ans[maxn];//ans[i]表示被覆盖了i次的面积 953     double pre[maxn]; 954     int n; 955     circles(){} 956     void add(circle cc) 957     { 958         c[n++]=cc; 959     } 960     bool inner(circle x,circle y) 961     { 962         if (x.relationcircle(y)!=1)return 0; 963         return dblcmp(x.r-y.r)<=0?1:0; 964     } 965     void init_or()//圆的面积并去掉内含的圆 966     { 967         int i,j,k=0; 968         bool mark[maxn]={0}; 969         for (i=0;i<n;i++) 970         { 971             for (j=0;j<n;j++)if (i!=j&&!mark[j]) 972             { 973                 if ((c[i]==c[j])||inner(c[i],c[j]))break; 974             } 975             if (j<n)mark[i]=1; 976         } 977         for (i=0;i<n;i++)if (!mark[i])c[k++]=c[i]; 978         n=k; 979     } 980     void init_and()//圆的面积交去掉内含的圆 981     { 982         int i,j,k=0; 983         bool mark[maxn]={0}; 984         for (i=0;i<n;i++) 985         { 986             for (j=0;j<n;j++)if (i!=j&&!mark[j]) 987             { 988                 if ((c[i]==c[j])||inner(c[j],c[i]))break; 989             } 990             if (j<n)mark[i]=1; 991         } 992         for (i=0;i<n;i++)if (!mark[i])c[k++]=c[i]; 993         n=k; 994     } 995     double areaarc(double th,double r) 996     { 997         return 0.5*sqr(r)*(th-sin(th)); 998     } 999     void getarea()1000     {1001         int i,j,k;1002         memset(ans,0,sizeof(ans));1003         vector<pair<double,int> >v;1004         for (i=0;i<n;i++)1005         {1006             v.clear();1007             v.push_back(make_pair(-pi,1));1008             v.push_back(make_pair(pi,-1));1009             for (j=0;j<n;j++)if (i!=j)1010             {1011                 point q=c[j].p.sub(c[i].p);1012                 double ab=q.len(),ac=c[i].r,bc=c[j].r;1013                 if (dblcmp(ab+ac-bc)<=0)1014                 {1015                     v.push_back(make_pair(-pi,1));1016                     v.push_back(make_pair(pi,-1));1017                     continue;1018                 }1019                 if (dblcmp(ab+bc-ac)<=0)continue;1020                 if (dblcmp(ab-ac-bc)>0) continue;1021                 double th=atan2(q.y,q.x),fai=acos((ac*ac+ab*ab-bc*bc)/(2.0*ac*ab));1022                 double a0=th-fai;1023                 if (dblcmp(a0+pi)<0)a0+=2*pi;1024                 double a1=th+fai;1025                 if (dblcmp(a1-pi)>0)a1-=2*pi;1026                 if (dblcmp(a0-a1)>0)1027                 {1028                     v.push_back(make_pair(a0,1));1029                     v.push_back(make_pair(pi,-1));1030                     v.push_back(make_pair(-pi,1));1031                     v.push_back(make_pair(a1,-1));1032                 }1033                 else1034                 {1035                     v.push_back(make_pair(a0,1));1036                     v.push_back(make_pair(a1,-1));1037                 }1038             }1039             sort(v.begin(),v.end());1040             int cur=0;1041             for (j=0;j<v.size();j++)1042             {1043                 if (cur&&dblcmp(v[j].first-pre[cur]))1044                 {1045                     ans[cur]+=areaarc(v[j].first-pre[cur],c[i].r);1046                     ans[cur]+=0.5*point(c[i].p.x+c[i].r*cos(pre[cur]),c[i].p.y+c[i].r*sin(pre[cur])).det(point(c[i].p.x+c[i].r*cos(v[j].first),c[i].p.y+c[i].r*sin(v[j].first)));1047                 }1048                 cur+=v[j].second;1049                 pre[cur]=v[j].first;1050             }1051         }1052         for (i=1;i<=n;i++)1053         {1054             ans[i]-=ans[i+1];1055         }1056     }1057 };1058 struct halfplane:public line1059 {1060     double angle;1061     halfplane(){}1062     //表示向量 a->b逆时针(左侧)的半平面1063     halfplane(point _a,point _b)1064     {1065         a=_a;1066         b=_b;1067     }1068     halfplane(line v)1069     {1070         a=v.a;1071         b=v.b;1072     }1073     void calcangle()1074     {1075         angle=atan2(b.y-a.y,b.x-a.x);1076     }1077     bool operator<(const halfplane &b)const1078     {1079         return angle<b.angle;1080     }1081 };1082 struct halfplanes1083 {1084     int n;1085     halfplane hp[maxp];1086     point p[maxp];1087     int que[maxp];1088     int st,ed;1089     void push(halfplane tmp)1090     {1091         hp[n++]=tmp;1092     }1093     void unique()1094     {1095         int m=1,i;1096         for (i=1;i<n;i++)1097         {1098             if (dblcmp(hp[i].angle-hp[i-1].angle))hp[m++]=hp[i];1099             else if (dblcmp(hp[m-1].b.sub(hp[m-1].a).det(hp[i].a.sub(hp[m-1].a))>0))hp[m-1]=hp[i];1100         }1101         n=m;1102     }1103     bool halfplaneinsert()1104     {1105         int i;1106         for (i=0;i<n;i++)hp[i].calcangle();1107         sort(hp,hp+n);1108         unique();1109         que[st=0]=0;1110         que[ed=1]=1;1111         p[1]=hp[0].crosspoint(hp[1]);1112         for (i=2;i<n;i++)1113         {1114             while (st<ed&&dblcmp((hp[i].b.sub(hp[i].a).det(p[ed].sub(hp[i].a))))<0)ed--;1115             while (st<ed&&dblcmp((hp[i].b.sub(hp[i].a).det(p[st+1].sub(hp[i].a))))<0)st++;1116             que[++ed]=i;1117             if (hp[i].parallel(hp[que[ed-1]]))return false;1118             p[ed]=hp[i].crosspoint(hp[que[ed-1]]);1119         }1120         while (st<ed&&dblcmp(hp[que[st]].b.sub(hp[que[st]].a).det(p[ed].sub(hp[que[st]].a)))<0)ed--;1121         while (st<ed&&dblcmp(hp[que[ed]].b.sub(hp[que[ed]].a).det(p[st+1].sub(hp[que[ed]].a)))<0)st++;1122         if (st+1>=ed)return false;1123         return true;1124     }1125     void getconvex(polygon &con)1126     {1127         p[st]=hp[que[st]].crosspoint(hp[que[ed]]);1128         con.n=ed-st+1;1129         int j=st,i=0;1130         for (;j<=ed;i++,j++)1131         {1132             con.p[i]=p[j];1133         }1134     }1135 };1136 struct point31137 {1138     double x,y,z;1139     point3(){}1140     point3(double _x,double _y,double _z):1141     x(_x),y(_y),z(_z){};1142     void input()1143     {1144         scanf("%lf%lf%lf",&x,&y,&z);1145     }1146     void output()1147     {1148         printf("%.2lf %.2lf %.2lf\n",x,y,z);1149     }1150     bool operator==(point3 a)1151     {1152         return dblcmp(a.x-x)==0&&dblcmp(a.y-y)==0&&dblcmp(a.z-z)==0;1153     }1154     bool operator<(point3 a)const1155     {1156         return dblcmp(a.x-x)==0?dblcmp(y-a.y)==0?dblcmp(z-a.z)<0:y<a.y:x<a.x;1157     }1158     double len()1159     {1160         return sqrt(len2());1161     }1162     double len2()1163     {1164         return x*x+y*y+z*z;1165     }1166     double distance(point3 p)1167     {1168         return sqrt((p.x-x)*(p.x-x)+(p.y-y)*(p.y-y)+(p.z-z)*(p.z-z));1169     }1170     point3 add(point3 p)1171     {1172         return point3(x+p.x,y+p.y,z+p.z);1173     }1174     point3 sub(point3 p)1175     {1176         return point3(x-p.x,y-p.y,z-p.z);1177     }1178     point3 mul(double d)1179     {1180         return point3(x*d,y*d,z*d);1181     }1182     point3 div(double d)1183     {1184         return point3(x/d,y/d,z/d);1185     }1186     double dot(point3 p)1187     {1188         return x*p.x+y*p.y+z*p.z;1189     }1190     point3 det(point3 p)1191     {1192         return point3(y*p.z-p.y*z,p.x*z-x*p.z,x*p.y-p.x*y);1193     }1194     double rad(point3 a,point3 b)1195     {1196         point3 p=(*this);1197         return acos(a.sub(p).dot(b.sub(p))/(a.distance(p)*b.distance(p)));1198     }1199     point3 trunc(double r)1200     {1201         r/=len();1202         return point3(x*r,y*r,z*r);1203     }1204     point3 rotate(point3 o,double r)1205     {1206     }1207 };1208 struct line31209 {1210     point3 a,b;1211     line3(){}1212     line3(point3 _a,point3 _b)1213     {1214         a=_a;1215         b=_b;1216     }1217     bool operator==(line3 v)1218     {1219         return (a==v.a)&&(b==v.b);1220     }1221     void input()1222     {1223         a.input();1224         b.input();1225     }1226     double length()1227     {1228         return a.distance(b);1229     }1230     bool pointonseg(point3 p)1231     {1232         return dblcmp(p.sub(a).det(p.sub(b)).len())==0&&dblcmp(a.sub(p).dot(b.sub(p)))<=0;1233     }1234     double dispointtoline(point3 p)1235     {1236         return b.sub(a).det(p.sub(a)).len()/a.distance(b);1237     }1238     double dispointtoseg(point3 p)1239     {1240         if (dblcmp(p.sub(b).dot(a.sub(b)))<0||dblcmp(p.sub(a).dot(b.sub(a)))<0)1241         {1242             return min(p.distance(a),p.distance(b));1243         }1244         return dispointtoline(p);1245     }1246     point3 lineprog(point3 p)1247     {1248         return a.add(b.sub(a).trunc(b.sub(a).dot(p.sub(a))/b.distance(a)));1249     }1250     point3 rotate(point3 p,double ang)//p绕此向量逆时针arg角度1251     {1252         if (dblcmp((p.sub(a).det(p.sub(b)).len()))==0)return p;1253         point3 f1=b.sub(a).det(p.sub(a));1254         point3 f2=b.sub(a).det(f1);1255         double len=fabs(a.sub(p).det(b.sub(p)).len()/a.distance(b));1256         f1=f1.trunc(len);f2=f2.trunc(len);1257         point3 h=p.add(f2);1258         point3 pp=h.add(f1);1259         return h.add((p.sub(h)).mul(cos(ang*1.0))).add((pp.sub(h)).mul(sin(ang*1.0)));1260     }1261 };1262 struct plane1263 {1264     point3 a,b,c,o;1265     plane(){}1266     plane(point3 _a,point3 _b,point3 _c)1267     {1268         a=_a;1269         b=_b;1270         c=_c;1271         o=pvec();1272     }1273     plane(double _a,double _b,double _c,double _d)1274     {1275         //ax+by+cz+d=01276         o=point3(_a,_b,_c);1277         if (dblcmp(_a)!=0)1278         {1279             a=point3((-_d-_c-_b)/_a,1,1);1280         }1281         else if (dblcmp(_b)!=0)1282         {1283             a=point3(1,(-_d-_c-_a)/_b,1);1284         }1285         else if (dblcmp(_c)!=0)1286         {1287             a=point3(1,1,(-_d-_a-_b)/_c);1288         }1289     }1290     void input()1291     {1292         a.input();1293         b.input();1294         c.input();1295         o=pvec();1296     }1297     point3 pvec()1298     {1299         return b.sub(a).det(c.sub(a));1300     }1301     bool pointonplane(point3 p)//点是否在平面上1302     {1303         return dblcmp(p.sub(a).dot(o))==0;1304     }1305     //0 不在1306     //1 在边界上1307     //2 在内部1308     int pointontriangle(point3 p)//点是否在空间三角形abc上1309     {1310         if (!pointonplane(p))return 0;1311         double s=a.sub(b).det(c.sub(b)).len();1312         double s1=p.sub(a).det(p.sub(b)).len();1313         double s2=p.sub(a).det(p.sub(c)).len();1314         double s3=p.sub(b).det(p.sub(c)).len();1315         if (dblcmp(s-s1-s2-s3))return 0;1316         if (dblcmp(s1)&&dblcmp(s2)&&dblcmp(s3))return 2;1317         return 1;1318     }1319     //判断两平面关系1320     //0 相交1321     //1 平行但不重合1322     //2 重合1323     bool relationplane(plane f)1324     {1325         if (dblcmp(o.det(f.o).len()))return 0;1326         if (pointonplane(f.a))return 2;1327         return 1;1328     }1329     double angleplane(plane f)//两平面夹角1330     {1331         return acos(o.dot(f.o)/(o.len()*f.o.len()));1332     }1333     double dispoint(point3 p)//点到平面距离1334     {1335         return fabs(p.sub(a).dot(o)/o.len());1336     }1337     point3 pttoplane(point3 p)//点到平面最近点1338     {1339         line3 u=line3(p,p.add(o));1340         crossline(u,p);1341         return p;1342     }1343     int crossline(line3 u,point3 &p)//平面和直线的交点1344     {1345         double x=o.dot(u.b.sub(a));1346         double y=o.dot(u.a.sub(a));1347         double d=x-y;1348         if (dblcmp(fabs(d))==0)return 0;1349         p=u.a.mul(x).sub(u.b.mul(y)).div(d);1350         return 1;1351     }1352     int crossplane(plane f,line3 &u)//平面和平面的交线1353     {1354         point3 oo=o.det(f.o);1355         point3 v=o.det(oo);1356         double d=fabs(f.o.dot(v));1357         if (dblcmp(d)==0)return 0;1358         point3 q=a.add(v.mul(f.o.dot(f.a.sub(a))/d));1359         u=line3(q,q.add(oo));1360         return 1;1361     }1362 };

  1 #include<vector>  2 #include<list>  3 #include<map>  4 #include<set>  5 #include<deque>  6 #include<queue>  7 #include<stack>  8 #include<bitset>  9 #include<algorithm> 10 #include<functional> 11 #include<numeric> 12 #include<utility> 13 #include<iostream> 14 #include<sstream> 15 #include<iomanip> 16 #include<cstdio> 17 #include<cmath> 18 #include<cstdlib> 19 #include<cctype> 20 #include<string> 21 #include<cstring> 22 #include<cstdio> 23 #include<cmath> 24 #include<cstdlib> 25 #include<ctime> 26 #include<climits> 27 #include<complex> 28 #define mp make_pair 29 #define pb push_back 30 using namespace std; 31 const double eps=1e-8;//精度 32 const double pi=acos(-1.0);//π 33 const double inf=1e20;//无穷大 34 const int maxp=1111;//最大点数 35 /* 36     判断d是否在精度内等于0 37 */ 38 int dblcmp(double d) 39 { 40     if (fabs(d)<eps)return 0; 41     return d>eps?1:-1; 42 } 43 /* 44     求x的平方 45 */ 46 inline double sqr(double x){return x*x;} 47 /* 48     点/向量 49 */ 50 struct point 51 { 52     double x,y; 53     point(){} 54     point(double _x,double _y):x(_x),y(_y){}; 55     //读入一个点 56     void input() 57     { 58         scanf("%lf%lf",&x,&y); 59     } 60     //输出一个点 61     void output() 62     { 63         printf("%.2f %.2f\n",x,y); 64     } 65     //判断两点是否相等 66     bool operator==(point a)const 67     { 68         return dblcmp(a.x-x)==0&&dblcmp(a.y-y)==0; 69     } 70     //判断两点大小 71     bool operator<(point a)const 72     { 73         return dblcmp(a.x-x)==0?dblcmp(y-a.y)<0:x<a.x; 74     } 75     //点到源点的距离/向量的长度 76     double len() 77     { 78         return hypot(x,y); 79     } 80     //点到源点距离的平方 81     double len2() 82     { 83         return x*x+y*y; 84     } 85     //两点间的距离 86     double distance(point p) 87     { 88         return hypot(x-p.x,y-p.y); 89     } 90     //向量加 91     point add(point p) 92     { 93         return point(x+p.x,y+p.y); 94     } 95     //向量减 96     point sub(point p) 97     { 98         return point(x-p.x,y-p.y); 99     }100     //向量乘101     point mul(double b)102     {103         return point(x*b,y*b);104     }105     //向量除106     point div(double b)107     {108         return point(x/b,y/b);109     }110     //点乘111     double dot(point p)112     {113         return x*p.x+y*p.y;114     }115     //叉乘116     double det(point p)117     {118         return x*p.y-y*p.x;119     }120     //XXXXXXX121     double rad(point a,point b)122     {123         point p=*this;124         return fabs(atan2(fabs(a.sub(p).det(b.sub(p))),a.sub(p).dot(b.sub(p))));125     }126     //截取长度r127     point trunc(double r)128     {129         double l=len();130         if (!dblcmp(l))return *this;131         r/=l;132         return point(x*r,y*r);133     }134     //左转90度135     point rotleft()136     {137         return point(-y,x);138     }139     //右转90度140     point rotright()141     {142         return point(y,-x);143     }144     //绕点p逆时针旋转angle角度145     point rotate(point p,double angle)146     {147         point v=this->sub(p);148         double c=cos(angle),s=sin(angle);149         return point(p.x+v.x*c-v.y*s,p.y+v.x*s+v.y*c);150     }151 };152 /*153     线段/直线154 */155 struct line156 {157     point a,b;158     line(){}159     line(point _a,point _b)160     {161         a=_a;162         b=_b;163     }164     //判断线段相等165     bool operator==(line v)166     {167         return (a==v.a)&&(b==v.b);168     }169     //点p做倾斜角为angle的射线170     line(point p,double angle)171     {172         a=p;173         if (dblcmp(angle-pi/2)==0)174         {175             b=a.add(point(0,1));176         }177         else178         {179             b=a.add(point(1,tan(angle)));180         }181     }182     //直线一般式ax+by+c=0183     line(double _a,double _b,double _c)184     {185         if (dblcmp(_a)==0)186         {187             a=point(0,-_c/_b);188             b=point(1,-_c/_b);189         }190         else if (dblcmp(_b)==0)191         {192             a=point(-_c/_a,0);193             b=point(-_c/_a,1);194         }195         else196         {197             a=point(0,-_c/_b);198             b=point(1,(-_c-_a)/_b);199         }200     }201     //读入一个线段202     void input()203     {204         a.input();205         b.input();206     }207     //校准线段两点208     void adjust()209     {210         if (b<a)swap(a,b);211     }212     //线段长度213     double length()214     {215         return a.distance(b);216     }217     //直线倾斜角 0<=angle<180218     double angle()219     {220         double k=atan2(b.y-a.y,b.x-a.x);221         if (dblcmp(k)<0)k+=pi;222         if (dblcmp(k-pi)==0)k-=pi;223         return k;224     }225     //点和线段关系226     //1 在逆时针227     //2 在顺时针228     //3 平行229     int relation(point p)230     {231         int c=dblcmp(p.sub(a).det(b.sub(a)));232         if (c<0)return 1;233         if (c>0)return 2;234         return 3;235     }236     //点是否在线段上237     bool pointonseg(point p)238     {239         return dblcmp(p.sub(a).det(b.sub(a)))==0&&dblcmp(p.sub(a).dot(p.sub(b)))<=0;240     }241     //两线是否平行242     bool parallel(line v)243     {244         return dblcmp(b.sub(a).det(v.b.sub(v.a)))==0;245     }246     //线段和线段关系247     //0 不相交248     //1 非规范相交249     //2 规范相交250     int segcrossseg(line v)251     {252         int d1=dblcmp(b.sub(a).det(v.a.sub(a)));253         int d2=dblcmp(b.sub(a).det(v.b.sub(a)));254         int d3=dblcmp(v.b.sub(v.a).det(a.sub(v.a)));255         int d4=dblcmp(v.b.sub(v.a).det(b.sub(v.a)));256         if ((d1^d2)==-2&&(d3^d4)==-2)return 2;257         return (d1==0&&dblcmp(v.a.sub(a).dot(v.a.sub(b)))<=0||258                 d2==0&&dblcmp(v.b.sub(a).dot(v.b.sub(b)))<=0||259                 d3==0&&dblcmp(a.sub(v.a).dot(a.sub(v.b)))<=0||260                 d4==0&&dblcmp(b.sub(v.a).dot(b.sub(v.b)))<=0);261     }262     //线段和直线v关系263     int linecrossseg(line v)//*this seg v line264     {265         int d1=dblcmp(b.sub(a).det(v.a.sub(a)));266         int d2=dblcmp(b.sub(a).det(v.b.sub(a)));267         if ((d1^d2)==-2)return 2;268         return (d1==0||d2==0);269     }270     //直线和直线关系271     //0 平行272     //1 重合273     //2 相交274     int linecrossline(line v)275     {276         if ((*this).parallel(v))277         {278             return v.relation(a)==3;279         }280         return 2;281     }282     //求两线交点283     point crosspoint(line v)284     {285         double a1=v.b.sub(v.a).det(a.sub(v.a));286         double a2=v.b.sub(v.a).det(b.sub(v.a));287         return point((a.x*a2-b.x*a1)/(a2-a1),(a.y*a2-b.y*a1)/(a2-a1));288     }289     //点p到直线的距离290     double dispointtoline(point p)291     {292         return fabs(p.sub(a).det(b.sub(a)))/length();293     }294     //点p到线段的距离295     double dispointtoseg(point p)296     {297         if (dblcmp(p.sub(b).dot(a.sub(b)))<0||dblcmp(p.sub(a).dot(b.sub(a)))<0)298         {299             return min(p.distance(a),p.distance(b));300         }301         return dispointtoline(p);302     }303     //XXXXXXXX304     point lineprog(point p)305     {306         return a.add(b.sub(a).mul(b.sub(a).dot(p.sub(a))/b.sub(a).len2()));307     }308     //点p关于直线的对称点309     point symmetrypoint(point p)310     {311         point q=lineprog(p);312         return point(2*q.x-p.x,2*q.y-p.y);313     }314 };315 316 /*317     多边形318 */319 struct polygon320 {321     int n;//点个数322     point p[maxp];//顶点323     line l[maxp];//边324     //读入一个多边形325     void input(int n=4)326     {327         for (int i=0;i<n;i++)328         {329             p[i].input();330         }331     }332     //添加一个点333     void add(point q)334     {335         p[n++]=q;336     }337     //取得边338     void getline()339     {340         for (int i=0;i<n;i++)341         {342             l[i]=line(p[i],p[(i+1)%n]);343         }344     }345     struct cmp346     {347         point p;348         cmp(const point &p0){p=p0;}349         bool operator()(const point &aa,const point &bb)350         {351             point a=aa,b=bb;352             int d=dblcmp(a.sub(p).det(b.sub(p)));353             if (d==0)354             {355                 return dblcmp(a.distance(p)-b.distance(p))<0;356             }357             return d>0;358         }359     };360     void norm()361     {362         point mi=p[0];363         for (int i=1;i<n;i++)mi=min(mi,p[i]);364         sort(p,p+n,cmp(mi));365     }366     //求凸包存入多边形convex367     void getconvex(polygon &convex)368     {369         int i,j,k;370         sort(p,p+n);371         convex.n=n;372         for (i=0;i<min(n,2);i++)373         {374             convex.p[i]=p[i];375         }376         if (n<=2)return;377         int &top=convex.n;378         top=1;379         for (i=2;i<n;i++)380         {381             while (top&&convex.p[top].sub(p[i]).det(convex.p[top-1].sub(p[i]))<=0)382                 top--;383             convex.p[++top]=p[i];384         }385         int temp=top;386         convex.p[++top]=p[n-2];387         for (i=n-3;i>=0;i--)388         {389             while (top!=temp&&convex.p[top].sub(p[i]).det(convex.p[top-1].sub(p[i]))<=0)390                 top--;391             convex.p[++top]=p[i];392         }393     }394     //判断是否凸多边形395     bool isconvex()396     {397         bool s[3];398         memset(s,0,sizeof(s));399         int i,j,k;400         for (i=0;i<n;i++)401         {402             j=(i+1)%n;403             k=(j+1)%n;404             s[dblcmp(p[j].sub(p[i]).det(p[k].sub(p[i])))+1]=1;405             if (s[0]&&s[2])return 0;406         }407         return 1;408     }409     //点与多边形关系410     //0 外部411     //1 内部412     //2 边上413     //3 点上414     int relationpoint(point q)415     {416         int i,j;417         for (i=0;i<n;i++)418         {419             if (p[i]==q)return 3;420         }421         getline();422         for (i=0;i<n;i++)423         {424             if (l[i].pointonseg(q))return 2;425         }426         int cnt=0;427         for (i=0;i<n;i++)428         {429             j=(i+1)%n;430             int k=dblcmp(q.sub(p[j]).det(p[i].sub(p[j])));431             int u=dblcmp(p[i].y-q.y);432             int v=dblcmp(p[j].y-q.y);433             if (k>0&&u<0&&v>=0)cnt++;434             if (k<0&&v<0&&u>=0)cnt--;435         }436         return cnt!=0;437     }438     //线段与多边形关系439     //0 无任何交点440     //1 在多边形内长度为正441     //2 相交或与边平行442     int relationline(line u)443     {444         int i,j,k=0;445         getline();446         for (i=0;i<n;i++)447         {448             if (l[i].segcrossseg(u)==2)return 1;449             if (l[i].segcrossseg(u)==1)k=1;450         }451         if (!k)return 0;452         vector<point>vp;453         for (i=0;i<n;i++)454         {455             if (l[i].segcrossseg(u))456             {457                 if (l[i].parallel(u))458                 {459                     vp.pb(u.a);460                     vp.pb(u.b);461                     vp.pb(l[i].a);462                     vp.pb(l[i].b);463                     continue;464                 }465                 vp.pb(l[i].crosspoint(u));466             }467         }468         sort(vp.begin(),vp.end());469         int sz=vp.size();470         for (i=0;i<sz-1;i++)471         {472             point mid=vp[i].add(vp[i+1]).div(2);473             if (relationpoint(mid)==1)return 1;474         }475         return 2;476     }477     //直线u切割凸多边形左侧478     //注意直线方向479     void convexcut(line u,polygon &po)480     {481         int i,j,k;482         int &top=po.n;483         top=0;484         for (i=0;i<n;i++)485         {486             int d1=dblcmp(p[i].sub(u.a).det(u.b.sub(u.a)));487             int d2=dblcmp(p[(i+1)%n].sub(u.a).det(u.b.sub(u.a)));488             if (d1>=0)po.p[top++]=p[i];489             if (d1*d2<0)po.p[top++]=u.crosspoint(line(p[i],p[(i+1)%n]));490         }491     }492     //取得周长493     double getcircumference()494     {495         double sum=0;496         int i;497         for (i=0;i<n;i++)498         {499             sum+=p[i].distance(p[(i+1)%n]);500         }501         return sum;502     }503     //取得面积504     double getarea()505     {506         double sum=0;507         int i;508         for (i=0;i<n;i++)509         {510             sum+=p[i].det(p[(i+1)%n]);511         }512         return fabs(sum)/2;513     }514     bool getdir()//1代表逆时针 0代表顺时针515     {516         double sum=0;517         int i;518         for (i=0;i<n;i++)519         {520             sum+=p[i].det(p[(i+1)%n]);521         }522         if (dblcmp(sum)>0)return 1;523         return 0;524     }525     //取得重心526     point getbarycentre()527     {528         point ret(0,0);529         double area=0;530         int i;531         for (i=1;i<n-1;i++)532         {533             double tmp=p[i].sub(p[0]).det(p[i+1].sub(p[0]));534             if (dblcmp(tmp)==0)continue;535             area+=tmp;536             ret.x+=(p[0].x+p[i].x+p[i+1].x)/3*tmp;537             ret.y+=(p[0].y+p[i].y+p[i+1].y)/3*tmp;538         }539         if (dblcmp(area))ret=ret.div(area);540         return ret;541     }542     //点在凸多边形内部的判定543     int pointinpolygon(point q)544     {545         if (getdir())reverse(p,p+n);546         if (dblcmp(q.sub(p[0]).det(p[n-1].sub(p[0])))==0)547         {548             if (line(p[n-1],p[0]).pointonseg(q))return n-1;549             return -1;550         }551         int low=1,high=n-2,mid;552         while (low<=high)553         {554             mid=(low+high)>>1;555             if (dblcmp(q.sub(p[0]).det(p[mid].sub(p[0])))>=0&&dblcmp(q.sub(p[0]).det(p[mid+1].sub(p[0])))<0)556             {557                 polygon c;558                 c.p[0]=p[mid];559                 c.p[1]=p[mid+1];560                 c.p[2]=p[0];561                 c.n=3;562                 if (c.relationpoint(q))return mid;563                 return -1;564             }565             if (dblcmp(q.sub(p[0]).det(p[mid].sub(p[0])))>0)566             {567                 low=mid+1;568             }569             else570             {571                 high=mid-1;572             }573         }574         return -1;575     }576 };