计算几何模板

Posted 2020-06-19 HugeGun

  1 #include<algorithm>
  2 #include<iostream>
  3 #include<cstdlib>
  4 #include<cstring>
  5 #include<cstdio>
  6 #include<string>
  7 #include<cmath>
  8 #include<ctime>
  9 #include<queue>
 10 #include<stack>
 11 #include<map>
 12 #include<set>
 13 #define rre(i,r,l) for(int i=(r);i>=(l);i--)
 14 #define re(i,l,r) for(int i=(l);i<=(r);i++)
 15 #define Clear(a,b) memset(a,b,sizeof(a))
 16 #define inout(x) printf("%d",(x))
 17 #define douin(x) scanf("%lf",&x)
 18 #define strin(x) scanf("%s",(x))
 19 #define LLin(x) scanf("%lld",&x)
 20 #define op operator
 21 #define CSC main
 22 typedef unsigned long long ULL;
 23 typedef const int cint;
 24 typedef long long LL;
 25 using namespace std;
 26 void inin(int &ret)
 27 {
 28     ret=0;int f=0;char ch=getchar();
 29     while(ch<‘0‘||ch>‘9‘){if(ch==‘-‘)f=1;ch=getchar();}
 30     while(ch>=‘0‘&&ch<=‘9‘)ret*=10,ret+=ch-‘0‘,ch=getchar();
 31     ret=f?-ret:ret;
 32 }
 33 const double eps=1e-8;
 34 const double pi=acos(-1);
 35 const double inf=2147483647;
 36 double f(const long double &a){return a*a;}
 37 inline int dcmp(const long double &a){return abs(a)<eps?0:a<0?-1:1;}
 38 struct xl
 39 {
 40     double x,y;
 41     xl(double x=0.,double y=0.):x(x),y(y){}
 42     void in(){douin(x),douin(y);}
 43     bool op == (const xl &rhs)const {return !dcmp(x-rhs.x)&&!dcmp(y-rhs.y);}
 44     bool op  < (const xl &rhs)const {return !dcmp(x-rhs.x)?y<rhs.y:x<rhs.x;}
 45     xl op + (const xl &rhs)const {return xl(x+rhs.x,y+rhs.y);}
 46     xl op - (const xl &rhs)const {return xl(x-rhs.x,y-rhs.y);}
 47     xl op * (const double &rhs)const {return xl(x*rhs,y*rhs);}
 48     xl op / (const double &rhs)const {return xl(x/rhs,y/rhs);}
 49     xl rotate(const double &rad)const {return xl(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad));}
 50     double len()const {return sqrt(f(x)+f(y));}
 51     double angle()const {return atan2(y,x);}
 52 };
 53 double D_(const xl &a,const xl &b){return a.x*b.x+a.y*b.y;}
 54 double X_(const xl &a,const xl &b){return a.x*b.y-a.y*b.x;}
 55 double dis(const xl &a,const xl &b){return sqrt(f(a.x-b.x)+f(a.y-b.y));}
 56 double dis2(const xl &a,const xl &b){return f(a.x-b.x)+f(a.y-b.y);}
 57 double angle(const xl &a,const xl &b){return acos(D_(a,b)/a.len()/b.len());}
 58 struct LINE
 59 {
 60     xl u,v;double rad;
 61     LINE(xl u=xl(0.,0.),xl v=xl(0.,0.)):u(u),v(v){}
 62     void js(){rad=(v-u).angle();}
 63     bool op < (const LINE &rhs)const {return rad<rhs.rad;}
 64     double k()const {return !dcmp(u.x-v.x)?inf:(v.y-u.y)/(v.x-u.x);}
 65     double len()const {return (v-u).len();}
 66 };
 67 xl getjd(const LINE &a,const LINE &b)
 68 {
 69     xl u=a.u-b.u,v1=a.v-a.u,v2=b.v-b.u;
 70     long double t=X_(v1,u)/X_(v1,v2);
 71     return a.u+v1*t;
 72 }
 73 bool pdxj(const LINE &a,const LINE &b)
 74 {
 75     if(!dcmp(X_(a.u,b.v-b.u)-X_(a.v,b.v-b.u)))return 0;
 76     return 1;
 77 }
 78 double distoline(const xl &a,const LINE &l)
 79 {
 80     return abs(X_(a,l.v-l.u))/l.len();
 81 }
 82 double distoseg(const xl &a,const LINE &l)
 83 {
 84     xl v1=a-l.u,v2=a-l.v,u=l.v-l.u;
 85     int t1=dcmp(X_(v1,u)),t2=dcmp(X_(v2,u));
 86     if(t1>0&&t2>0)return v2.len();
 87     if(t1<0&&t2<0)return v1.len();
 88     return distoline(a,l);
 89 }
 90 xl getpoty(const xl &a,const LINE &l)
 91 {
 92     xl v=l.v-l.u;
 93     return l.u+v*D_(a-l.u,v)/D_(v,v);
 94 }
 95 bool pdxjseg(const LINE &a,const LINE &b)
 96 {
 97     xl u1=a.u-b.u,u2=a.v-b.u,u=a.v-a.u;
 98     xl v1=b.u-a.u,v2=b.v-a.u,v=b.v-b.u;
 99     int t1=dcmp(X_(u1,v)),t2=dcmp(X_(u2,v));
100     int t3=dcmp(X_(v1,u)),t4=dcmp(X_(v2,v));
101     return t1*t2<=0&&t3*t4<=0;
102 }
103 bool onseg(const xl &a,const LINE &l)
104 {
105     return !dcmp(D_(a-l.u,a-l.v))&&dcmp(X_(a-l.u,a-l.v))<=0;
106 }
107 struct cir
108 {
109     xl c;double r;
110     cir(xl c=xl(0.,0.),double r=0):c(c),r(r){}
111     bool op < (const cir &rhs)const {return c==rhs.c?r<rhs.r:c<rhs.c;}
112     xl point (const double &rad)const {return c+xl(r*cos(rad),r*sin(rad));}
113 };
114 int getjd(const LINE &l,const cir &c,xl &ans1,xl &ans2)
115 {
116     double d=distoline(c.c,l);
117     if(dcmp(d-c.r)>0)return 0;
118     if(!dcmp(d-c.r))return ans1=getpoty(c.c,l),1;
119     xl temp=getpoty(c.c,l);
120     xl v=c.c-temp;swap(v.x,v.y);
121     double len=sqrt(f(c.r)-f(v.len()));
122     v=v*len/v.len();
123     ans1=temp+v,ans2=temp-v;
124     return 2;
125 }
126 int getjd(const cir &a,const cir &b,xl &ans1,xl &ans2)
127 {
128     double d=(b.c-a.c).len();
129     if(!dcmp(d))if(!dcmp(b.r-a.r))return inf;else return 0;else ;
130     if(dcmp(a.r+b.r-d)<0||dcmp(fabs(a.r-b.r)-d)>0)return 0;
131     double rad=(b.c-a.c).angle();
132     double da=acos((f(a.r)+f(d)-f(b.r))/(2*a.r*d));
133     ans1=a.point(rad-da);
134     xl ans3=a.point(rad+da);
135     if(ans1==ans3)return 1;
136     ans2=ans3;return 2;
137 }
138 int gettangents(cir &A,cir &B,LINE *ans)
139 {
140     if(A.c==B.c)if(!dcmp(A.r-B.r))return inf;
141     else return 0;else ;
142     if(A.r<B.r)swap(A,B);
143     double dd=f(A.c.x-B.c.x)+f(A.c.y-B.c.y);
144     double rdiff=A.r-B.r,rsum=A.r+B.r;
145     if(dcmp(dd-f(rdiff))<0)return 0;
146     xl u=B.c-A.c;
147     double base=atan2(B.c.y-A.c.y,B.c.x-A.c.x);
148     if(!dcmp(dd-f(rdiff)))
149     {
150         xl p=A.point(base);
151         ans[1]=LINE(p,p+u.rotate(pi/2));return 1;
152     }
153     double rad=acos((A.r-B.r)/sqrt(dd));
154     ans[1]=LINE(A.point(base+rad),B.point(base+rad));
155     ans[2]=LINE(A.point(base-rad),B.point(base-rad));
156     if(dcmp(dd-f(rsum))<0)return 2;
157     if(!dcmp(dd-f(rsum)))
158     {
159         xl p=A.point(base);
160         ans[3]=LINE(p,p+u.rotate(pi/2));
161         return 3;
162     }
163     rad=acos((A.r+B.r)/sqrt(dd));
164     ans[3]=LINE(A.point(base+rad),B.point(base+rad+pi));
165     ans[4]=LINE(A.point(base-rad),B.point(base-rad+pi));
166     return 4;
167 }
168 cir outoftriangle(const xl &A,const xl &B,const xl &C)
169 {
170     double bx=B.x-A.x,by=B.y-A.y;
171     double cx=C.x-A.x,cy=C.y-A.y;
172     double area=2*(bx*cy-by*cx);
173     double x=(cy*(f(bx)+f(by))-by*(f(cx)+f(cy)))/area+A.x;
174     double y=(bx*(f(cx)+f(cy))-cx*(f(bx)+f(by)))/area+A.y;
175     return cir(xl(x,y),(A-xl(x,y)).len());
176 }
177 cir inoftriangle(const xl &A,const xl &B,const xl &C)
178 {
179     double a=(B-C).len(),b=(C-A).len(),c=(A-B).len();
180     xl cc=(A*a+B*b+C*c)/(a+b+c);
181     return cir(cc,distoline(cc,LINE(A,B)));
182 }
183 struct pol
184 {
185     xl a[100010];int n;
186     void in(){inin(n);re(i,0,n-1)a[i].in();}
187     void Sort(){sort(a,a+n);}
188     xl& op [] (int x){return a[x];}
189 };
190 int pdinpol(const xl &p,pol &a)
191 {
192     int sum=0,n=a.n;
193     re(i,0,n)
194     {
195         if(onseg(p,LINE(a[i],a[(i+1)%n])))return -1;
196         int k=dcmp(X_(a[(i+1)%n]-a[i],p-a[i]));
197         int d1=dcmp(a[i].y-p.y);
198         int d2=dcmp(a[(i+1)%n].y-p.y);
199         if(k>0&&d1<=0&&d2>0)sum++;
200         if(k<0&&d2<=0&&d1>0)sum--;
201     }
202     if(sum!=0)return 1;
203     return 0;
204 }
205 void hull(pol &p,pol &ch)
206 {
207     int n=p.n,m=0;
208     p.Sort();
209     re(i,0,n-1)
210     {
211         while(m>1&&dcmp(X_(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<0)m--;
212         ch[m++]=p[i];
213     }
214     int k=m;
215     rre(i,n-1,0)
216     {
217         while(m>k&&dcmp(X_(ch[m-1]-ch[m-2],p[i]-ch[m-2]))<0)m--;
218         ch[m++]=p[i];
219     }
220     if(n>1)m--;
221     ch.n=m;
222 }
223 double lenofpol(pol &p)
224 {
225     double ret=0;
226     re(i,0,p.n-1)ret+=(p[(i+1)%p.n]-p[i]).len();
227     return ret;
228 }
229 double areaofpol(pol &p)
230 {
231     double ret=0;int n=p.n;
232     re(i,1,n-2)ret+=X_(p[i]-p[0],p[i+1]-p[0]);
233     return ret/2;
234 }
235 double diameterofhull(pol &p)
236 {
237     int n=p.n;double ret=0.;
238     if(n<=1)return 0.;
239     if(n==2)return dis2(p[0],p[1]);
240     int r=1;
241     re(l,0,n-1)while(1)
242     {
243         double d=X_(p[l+1]-p[l],p[r+1]-p[r]);
244         if(dcmp(d)<=0)
245         {
246             ret=max(ret,dis2(p[l],p[r]));
247             if(!dcmp(d))ret=max(ret,dis2(p[l],p[r+1]));
248             break;
249         }
250         r++,r%=n;
251     }
252     return ret;
253 }
254 bool onleft(const xl &a,const LINE &l){return X_(l.v-l.u,a-l.u)>0;}
255 bool halfplanej(LINE *l,int n,pol &ch)
256 {
257     re(i,0,n-1)l[i].js();
258     sort(l,l+n);
259     int ll=0,rr=0;
260     xl p[n];
261     LINE q[n];q[0]=l[0];
262     re(i,1,n-1)
263     {
264         while(ll<rr&&!onleft(p[rr-1],l[i]))rr--;
265         while(ll<rr&&!onleft(p[ll],l[i]))ll++;
266         q[++rr]=l[i];
267         if(!dcmp(X_(q[rr].v-q[rr].u,q[rr-1].v-q[rr-1].u)))
268         {
269             rr--;
270             if(onleft(l[i].u,q[rr]))q[rr]=l[i];
271         }
272         if(ll<rr)p[rr-1]=getjd(q[rr-1],q[rr]);
273     }
274     while(ll<rr&&!onleft(p[rr-1],q[ll]))rr--;
275     if(rr-ll<=1)return 0;
276     p[rr]=getjd(q[rr],q[ll]);
277     re(i,ll,rr)ch[ch.n++]=p[i];
278     ch[ch.n]=ch[0];return 1;
279 }
280 int main()
281 {
282      return 0;
283 }