有错指出,长期更新
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<complex>
#include<vector>
#define db double
using namespace std ;
namespace Complex{
typedef complex <double > Point;
typedef Point Vector;
db dot (Vector a,Vector b){return real(conj(a)*b);}
db cross(Vector a,Vector b){return imag(conj(a)*b);}
Vector rotate(Vector a,db rad){return a*exp (Point(0 ,rad));}
}
namespace geometry{
db eps=1e-8 ;
int dcmp(db x){if (fabs (x)<eps)return 0 ;return x<0 ?-1 :1 ;}
const db pi=acos (-1 );
struct Point{
db x,y;
Point(db a=0 ,db b=0 ):x(a),y(b){}
void in(){scanf ("%lf%lf" ,&x,&y);}
};
typedef Point Vector;
Vector operator +(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
Vector operator -(Point a,Point b){return Vector(a.x-b.x,a.y-b.y);}
Vector operator *(Vector a,db b){return Vector(a.x*b,a.y*b); }
Vector operator /(Vector a,db b){return Vector(a.x/b,a.y/b); }
bool operator < (Point a,Point b){return a.x<b.x||(a.x==b.x&&a.y<b.y); }
bool operator ==(Point a,Point b){return dcmp(a.x-b.x)==0 &&dcmp(a.y-b.y)==0 ;}
db cross (Vector a,Vector b){return a.x*b.y-a.y*b.x;}
db dot (Vector a,Vector b){return a.x*b.x+a.y*b.y;}
db length (Vector a) {return sqrt (dot(a,a)); }
db angle (Vector a,Vector b){return acos (dot(a,b)/length(a)/length(b));}
db area (Point a,Point b,Point c){return cross(b-a,c-a); }
db area2 (Point &a,Point &b,Point &c){return cross(b-a,c-a); }
Vector normal(Vector a) {db l=length(a);return Vector(-a.y/l,a.x/l);}
Vector rotate(Point a,db rad){return Vector(a.x*cos (rad)-a.y*sin (rad),a.x*sin (rad)+a.y*cos (rad));}
Point getlinecut(Point a,Vector u,Point b,Vector v){
Vector w=a-b;
db t=cross(v,w)/cross(u,v);
return a+u*t;
}
db distoline(Point p,Point a,Point b){
Vector v1=b-a,v2=p-a;
return fabs (cross(v1,v2)/length(v1));
}
db distosegm(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);
if (dcmp(dot(v1,v3))>0 ) return length(v3);
return fabs (cross(v1,v2)/length(v1));
}
Point getlinepro(Point p,Point a,Point b){
Vector v=b-a;
return a+v*(dot(v,p-a)/dot(v,v));
}
bool segmpropcut(Point a1,Point a2,Point b1,Point b2){
db c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
db c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
return dcmp(c1)*dcmp(c2)<0 &&dcmp(c3)*dcmp(c4)<0 ;
}
bool onsegm(Point p,Point a,Point b){
return dcmp(cross(a-p,b-p))==0 &&dcmp(dot(a-p,b-p))<0 ;
}
db polyarea(Point *p,int n){
db ans=0 ;
for (int i=2 ;i<=n-1 ;i++)ans+=area(p[1 ],p[i],p[i+1 ]);
return ans/2 ;
}
db pangle(Vector v){return atan2 (v.y,v.x);}
struct line{
Vector v;
Point p;
db ang;
line(){}
line(Point a,Vector b):p(a),v(b){ang=atan2 (v.y,v.x);}
Point getp(db t){return p+v*t;}
bool operator <(line a)const {return ang<a.ang;}
};
Point getlinecut(line a,line b){return getlinecut(a.p,a.v,b.p,b.v);}
struct segment{
Point s,t;
Vector v;
db ang;
segment(){}
segment(Point a,Point b):s(a),t(b){v=b-a;ang=atan2 (v.y,v.x);}
};
bool segmpropcut(segment a,segment b){return segmpropcut(a.s,a.t,b.s,b.t);}
bool onsegm(Point p,segment seg){return onsegm(p,seg.s,seg.t);}
struct circle{
Point o;
db r;
circle(Point a,db b):o(a),r(b){}
Point getp(db rad){return Point(o.x+cos (rad)*r,o.y+sin (rad)*r);}
};
int getlinecirclecut(line l,circle cir,vector <Point> &ans){
db a=l.v.x,b=l.p.x-cir.o.x,c=l.v.y,d=l.p.y-cir.o.y;
db e=a*a+c*c,f=2 *(a*b+c*d),g=b*b+d*d-cir.r*cir.r;
db delta=f*f-4 *e*g;
if (dcmp(delta)<0 ) return 0 ;
if (dcmp(delta)==0 ) {
db t1=-f/(2 *e);
ans.push_back(l.getp(t1));
return 1 ;
}
db t1=(-f-sqrt (delta))/(2 *e);ans.push_back(l.getp(t1));
db t2=(-f+sqrt (delta))/(2 *e);ans.push_back(l.getp(t2));
return 2 ;
}
int getlinecirclecut(line l,circle cir,db &t1,db &t2,vector <Point> &ans){
db a=l.v.x,b=l.p.x-cir.o.x,c=l.v.y,d=l.p.y-cir.o.y;
db e=a*a+c*c,f=2 *(a*b+c*d),g=b*b+d*d-cir.r*cir.r;
db delta=f*f-4 *e*g;
if (dcmp(delta)<0 ) return 0 ;
if (dcmp(delta)==0 ) {
t1=t2=-f/(2 *e);
ans.push_back(l.getp(t1));
return 1 ;
}
t1=(-f-sqrt (delta))/(2 *e);ans.push_back(l.getp(t1));
t2=(-f+sqrt (delta))/(2 *e);ans.push_back(l.getp(t2));
return 2 ;
}
int getcircut(circle a,circle b,vector <Point> &ans){
db d=length(a.o-b.o);
if (dcmp(d)==0 ){if (dcmp(a.r-b.r)==0 ) return -1 ;return 0 ;}
if (dcmp(a.r+b.r-d)<0 ) return 0 ;
if (dcmp(fabs (a.r-b.r)-d)>0 ) return 0 ;
db rad=pangle(a.o-b.o);
db drad=acos ((a.r*a.r+d*d-b.r*b.r)/(2 *a.r*d));
Point p1=a.getp(rad-drad),p2=a.getp(rad+drad);
ans.push_back(p1);
if (p1==p2) return 1 ;
ans.push_back(p2);
return 2 ;
}
int gettan(Point a,circle c,Vector *ans){
Vector u=c.o-a;
db dist=length(u);
if (dist<c.r) return 0 ;
if (dcmp(dist-c.r)==0 ){
ans[0 ]=rotate(u,pi/2 );
return 1 ;
}
db ang=asin (c.r/dist);
ans[0 ]=rotate(u,-ang);
ans[1 ]=rotate(u,+ang);
return 2 ;
}
int gettan(Point a,circle c,vector <line> &ans){
Vector u=c.o-a;
db dist=length(u);
if (dist<c.r) return 0 ;
if (dcmp(dist-c.r)==0 ){
ans.push_back(line(a,rotate(u,pi/2 )));
return 1 ;
}
db ang=asin (c.r/dist);
ans.push_back(line(a,rotate(u,-ang)));
ans.push_back(line(a,rotate(u,+ang)));
return 2 ;
}
db dis2(Point a,Point b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);}
int getcirtan(circle c1,circle c2,Point *a,Point *b){
int cnt=0 ;
if (c1.r<c2.r) swap(c1,c2),swap(a,b);
Point o1=c1.o,o2=c2.o;
db d2=dis2(o1,o2);
db rdiff=c1.r-c2.r;
if (d2<rdiff*rdiff) return 0 ;
db base=pangle(o2-o1);
if (d2==0 &&c1.r==c2.r) return -1 ;
if (d2==rdiff*rdiff) {
a[++cnt]=c1.getp(base);b[cnt]=c2.getp(base);return 1 ;
}
db rsum=c1.r+c2.r;
db ang=acos ((c1.r-c2.r)/sqrt (d2));
a[++cnt]=c1.getp(base+ang);b[cnt]=c2.getp(base+ang);
a[++cnt]=c1.getp(base-ang);b[cnt]=c2.getp(base-ang);
if (d2==rsum*rsum){
ang=acos ((c1.r+c2.r)/sqrt (d2));
a[++cnt]=c1.getp(base+ang);b[cnt]=c2.getp(base+ang);
a[++cnt]=c1.getp(base-ang);b[cnt]=c2.getp(base-ang);
}
return cnt;
}
int getcirtan(circle a,circle b,vector <line> &ans){
Point *a1=new Point[5 ];
Point *b1=new Point[5 ];
int now=getcirtan(a,b,a1,b1);
for (int i=1 ;i<=now;i++){ans.push_back(line(a1[i],a1[i]-b1[i]));}
delete [] a1;
delete [] b1;
return now;
}
circle outcircle(Point p1,Point p2,Point p3){
db Bx = p2.x-p1.x, By = p2.y-p1.y;
db Cx = p3.x-p1.x, Cy = p3.y-p1.y;
db D = 2 *(Bx*Cy-By*Cx);
db cx = (Cy*(Bx*Bx+By*By) - By*(Cx*Cx+Cy*Cy))/D + p1.x;
db cy = (Bx*(Cx*Cx+Cy*Cy) - Cx*(Bx*Bx+By*By))/D + p1.y;
Point p = Point(cx, cy);
return circle(p, length(p1-p));
}
circle inscircle(Point p1,Point p2,Point p3){
db a = length(p2-p3);
db b = length(p3-p1);
db c = length(p1-p2);
Point p = (p1*a+p2*b+p3*c)/(a+b+c);
return circle(p, distoline(p, p1, p2));
}
struct Poly{
Point *p;
db are;
int sze,newsze;
bool isare;
void resize(){
if (sze<newsze-1 ) return ;
Point *now=p;
newsze=newsze<<1 |1 ;
p=new Point[newsze];
for (int i=1 ;i<=sze;i++)p[i]=now[i];
delete [] now;
}
Poly(){isare=0 ;}
Poly(int size){p=new Point[size];newsze=size;isare=0 ;}
void insert(Point a){
resize();
p[++sze]=a;isare=0 ;
}
int size(){return sze;}
db getarea(){
if (isare) return are;
are=0 ;
for (int i=2 ;i<=sze-1 ;i++){
are+=area(p[1 ],p[i],p[i+1 ]);
}
are/=2 ;
return are;
}
Point &operator [](int a){return p[a];}
void clear(){sze=0 ;delete [] p;p=new Point[5 ];newsze=5 ;}
~Poly(){delete [] p;}
};
struct VPoly{
vector <Point> p;
void insert(Point a){p.push_back(a);}
Point &operator [](int a){return p[a-1 ];}
int size(){return p.size();}
void clear(){p.clear();}
};
int ispointinpoly(Point p,Poly poly){
int wn=0 ;
int n=poly.size();
for (int i=1 ;i<=n;i++){
if (onsegm(p,poly[i],poly[(i+1 )%n])) return -1 ;
int k =dcmp(cross(poly[(i+1 )%n]-poly[i],p-poly[i]));
int d1=dcmp(poly[i].y-p.y);
int d2=dcmp(poly[(i+1 )%n].y-p.y);
if (k>0 &&d1<=0 &&d2>0 ) wn++;
if (k>0 &&d2<=0 &&d1>0 ) wn--;
}
if (wn!=0 ) return 1 ;
return 0 ;
}
int convexhull(Point *p,int n,Point *ans){
sort(p+1 ,p+n+1 );
int m=0 ;
for (int i=1 ;i<=n;i++){
while (m>2 &&dcmp(cross(ans[m-1 ]-ans[m-2 ],p[i]-ans[m-2 ]))<=0 ) m--;
ans[++m]=p[i];
}
int k=m;
for (int i=n-1 ;i>=1 ;i--){
while (m>=k&&dcmp(cross(ans[m-1 ]-ans[m-2 ],p[i]-ans[m-2 ]))<=0 ) m--;
ans[++m]=p[i];
}
return m;
}
int convexhull_noside(Point *p,int n,Point *ans){
sort(p+1 ,p+n+1 );
int m=0 ;
for (int i=1 ;i<=n;i++){
while (m>2 &&dcmp(cross(ans[m-1 ]-ans[m-2 ],p[i]-ans[m-2 ]))<0 ) m--;
ans[++m]=p[i];
}
int k=m;
for (int i=n-1 ;i>=1 ;i--){
while (m>=k&&dcmp(cross(ans[m-1 ]-ans[m-2 ],p[i]-ans[m-2 ]))<0 ) m--;
ans[++m]=p[i];
}
return m;
}
int convexhull_zero(Point *p,int n,Point *ans){
sort(p,p+n);
int m=0 ;
for (int i=0 ;i<n;i++){
while (m>1 &&dcmp(cross(ans[m-1 ]-ans[m-2 ],p[i]-ans[m-2 ]))<=0 ) m--;
ans[m++]=p[i];
}
int k=m;
for (int i=n-2 ;i>=0 ;i--){
while (m>=k&&dcmp(cross(ans[m-1 ]-ans[m-2 ],p[i]-ans[m-2 ]))<=0 ) m--;
ans[m++]=p[i];
}
if (n>1 )m--;
return m;
}
int convexhull_zero_noside(Point *p,int n,Point *ans){
sort(p,p+n);
int m=0 ;
for (int i=0 ;i<n;i++){
while (m>1 &&dcmp(cross(ans[m-1 ]-ans[m-2 ],p[i]-ans[m-2 ]))<=0 ) m--;
ans[m++]=p[i];
}
int k=m;
for (int i=n-2 ;i>=0 ;i--){
while (m>k&&dcmp(cross(ans[m-1 ]-ans[m-2 ],p[i]-ans[m-2 ]))<=0 ) m--;
ans[m++]=p[i];
}
if (n>1 )m--;
return m;
}
bool onleft(line l,Point p){
return dcmp(cross(l.v,p-l.p))>0 ;
}
bool onleft2(line l,Point p){
return dcmp(cross(l.v,p-l.p))>=0 ;
}
int halfcut(line *l,int n,Point *ans){
sort(l,l+n);
int fi,la;
Point *p= new Point[n];
line *q= new line[n];
q[fi=la=0 ]=l[0 ];
for (int i=1 ;i<n;i++){
while (fi<la&&!onleft(l[i],p[la-1 ])) la--;
while (fi<la&&!onleft(l[i],p[fi])) fi++;
q[++la]=l[i];
if (fabs (cross(q[la].v,q[la-1 ].v))<eps){
la--;
if (onleft(q[la],l[i].p)) q[la]=l[i];
}
if (fi<la) p[la-1 ]=getlinecut(q[la-1 ],q[la]);
}
while (fi<la&&!onleft(q[fi],p[la-1 ]))la--;
if (la-fi<=1 ) return 0 ;
p[la]=getlinecut(q[la],q[fi]);
int m=0 ;
for (int i=fi;i<=la;i++) ans[m++]=p[i];
delete [] p;
delete [] q;
return m;
}
int halfcut_one(line *l,int n,Point *ans){
sort(l+1 ,l+n+1 );
int fi,la;
Point *p= new Point[n];
line *q= new line[n];
q[fi=la=1 ]=l[1 ];
for (int i=2 ;i<=n;i++){
while (fi<la&&!onleft(l[i],p[la-1 ])) la--;
while (fi<la&&!onleft(l[i],p[fi])) fi++;
q[++la]=l[i];
if (fabs (cross(q[la].v,q[la].v))<eps){
la--;
if (onleft(q[la],l[i].p)) q[la]=l[i];
}
if (fi<la) p[la-1 ]=getlinecut(q[la-1 ],q[la]);
}
while (fi<la&&!onleft(q[fi],p[la-1 ]))la--;
if (la-fi<=1 ) return 0 ;
p[la]=getlinecut(q[la],q[fi]);
int m=0 ;
for (int i=fi;i<=la;i++) ans[++m]=p[i];
delete [] p;
delete [] q;
return m;
}
db dsts(Point a,Point b,Point c,Point d){
return min(min(distosegm(c,a,b),distosegm(d,a,b)),min(distosegm(a,c,d),distosegm(b,c,d)));
}
db get_two_poly_min_dis(Point *p1,Point *p2,int n,int m){
int pos1=0 ,pos2=0 ;
db ans=1e90 ;
for (int i=0 ;i<n;i++){
if (p1[i].y>=p1[pos1].y){
if (p1[i].y>p1[pos1].y) pos1=i;
else if (p1[i].x>p1[pos1].x) pos1=i;
}
}
for (int i=0 ;i<m;i++){
if (p2[i].y<=p2[pos2].y){
if (p2[i].y<p2[pos2].y) pos2=i;
else if (p2[i].x<p2[pos2].x) pos2=i;
}
}
db ls;
for (int i=0 ;i<n;i++){
while (ls=(cross(p1[(pos1+1 )%n]-p1[pos1],p2[(pos2+1 )%m]-p1[pos1])-cross(p1[(pos1+1 )%n]-p1[pos1],p2[pos2]-p1[pos1]))>eps)
pos2=(pos2+1 )%m;
if (ls<-eps) ans=min(ans,distosegm(p2[pos2],p1[(pos1+1 )%n],p1[pos1]));
else ans=min(ans,dsts(p1[(pos1+1 )%n],p1[pos1],p2[(pos2+1 )%m],p2[pos2]));
pos1=(pos1+1 )%n;
}
return ans;
}
db dis(Point a,Point b){return length(a-b);}
db RotatingCaliper(int m,Point *ch)
{
int q = 1 ;
ch[m] = ch[0 ];
db ans = 0 ;
for (int p = 0 ; p < m; p++)
{
while (fabs (cross(ch[q+1 ]-ch[p+1 ], ch[p]-ch[p+1 ])) > fabs (cross(ch[q]-ch[p+1 ], ch[p]-ch[p+1 ]))) q = (q+1 )%m;
ans = max(ans,max(dis(ch[p],ch[q]),dis(ch[p+1 ],ch[q+1 ])));
ans = max(ans,max(dis(ch[p],ch[q+1 ]),dis(ch[p+1 ],ch[q])));
}
return ans;
}
db torad(db deg){return deg/180 *pi;}
void get_coord(db r,db lat,db lng,db &x,db &y,db &z){
lat=torad(lat);
lng=torad(lng);
x=r*cos (lat)*cos (lng);
y=r*cos (lat)*sin (lng);
z=r*sin (lat);
}
struct Point3{
db x,y,z;
Point3(db a=0 ,db b=0 ,db c=0 ):x(a),y(b),z(c){}
};
typedef Point3 Vector3;
Vector3 operator +(Vector3 a,Vector3 b){return Vector3(a.x+b.x,a.y+b.y,a.z+b.z);}
Vector3 operator -(Point3 a,Point3 b){return Vector3(a.x-b.x,a.y-b.y,a.z-b.z);}
Vector3 operator *(Vector3 a,db b){return Vector3(a.x*b, a.y*b, a.z*b) ;}
Vector3 operator /(Vector3 a,db b){return Vector3(a.x/b, a.y/b, a.z/b) ;}
db dot3(Vector3 a,Vector3 b){return a.x*b.x+a.y*b.y+a.z*b.z;}
Vector3 cross(Vector3 a,Vector3 b){return Vector3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);}
}
using namespace geometry;
int main(){
return 0 ;
}