#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <vector>
#include <complex>
#include <queue>
#define eps 1e-10
#define pi acos(-1.0)
#define INF 0x7ffffff
#define maxn 50005
#define maxm 50005
#define next(i) ((i+1)%n)
using namespace std;
int cmp(double x)
{
if(fabs(x)<eps)
return 0;
if(x>0)
return 1;
return -1;
}
inline double sqr(double x)
{
return x*x;
}
struct point//点
{
double x,y;
point() {}
point(double a,double b):x(a),y(b) {}
void input()
{
scanf("%lf%lf",&x,&y);
}
friend point operator + (const point &a,const point &b)
{
return point(a.x+b.x,a.y+b.y);
}
friend point operator - (const point &a,const point &b)
{
return point(a.x-b.x,a.y-b.y);
}
friend bool operator == (const point &a,const point &b)
{
return cmp(a.x-b.x)==0 &&cmp(a.y-b.y)==0;
}
friend point operator * (const point &a,const double &b)
{
return point(a.x*b,a.y*b);
}
friend point operator * (const double &a,const point &b)
{
return point(a*b.x,a*b.y);
}
friend point operator / (const point &a,const double &b)
{
return point(a.x/b,a.y/b);
}
double norm()
{
return sqrt(sqr(x)+sqr(y));
}
};
double det (const point &a,const point &b)
{
return a.x*b.y-a.y*b.x;
}//叉积
double dot (const point &a,const point &b)
{
return a.x*b.x+a.y*b.y;
}//点乘
double dist (const point &a,const point &b)
{
return (a-b).norm();
}//距离
point rotate_point(const point &p,double A)
{
double tx=p.x,ty=p.y;
return point (tx*cos(A)-ty*sin(A),tx*sin(A)+ty*cos(A));
}//旋转
struct line
{
point a,b;
line() {}
line(point x,point y):a(x),b(y) {}
};
line point_make_line(const point a,const point b)
{
return line(a,b);
}//两个点a,b生成的一个线段
double dis_point_segment(const point p,const point s,const point t)
{
if(cmp(dot(p-s,t-s))<0)
return (p-s).norm();
if(cmp(dot(p-t,s-t))<0)
return (p-t).norm();
return fabs(det(s-p,t-p)/dist(s,t));
}//p到线段st的距离
void PointProLine(const point p,const point s,const point t,point &cp)
{
double r=dot((t-s),(p-s))/dot(t-s,t-s);
cp=s+r*(t-s);
}//p到线段st的垂足,存在cp
bool PointOnSegment (point p,point s,point t)
{
return cmp(det(p-s,t-s))==0&&cmp(dot(p-s,p-t))<=0;
}//点p是否在线段st上
bool parallel(line a,line b)
{
return !cmp(det(a.a-a.b,b.a-b.b));
}//平行
bool line_make_point (line a,line b,point &res)
{
if(parallel(a,b))
return false;
double s1=det(a.a-b.a,b.b-b.a);
double s2=det(a.b-b.a,b.b-b.a);
res=(s1*a.b-s2*a.a)/(s1-s2);
return true;
}//交点
line move_d(line a,const double &len)
{
point d=a.b-a.a;
d=d/d.norm();
d=rotate_point(d,pi/2);
return line(a.a+d*len,a.b+d*len);
}
struct polygon//多边形
{
int n;
point a[maxn];
polygon() {}
double perimeter()
{
double sum=0;
a[n]=a[0];
for(int i=0; i<n; i++)
sum+=(a[i+1]-a[i]).norm();
return sum;
}
double area()
{
double sum=0;
a[n]=a[0];
for(int i=0; i<n; i++)
sum+=det(a[i+1],a[i]);
return sum/2.0;
}
int Point_In(point t)
{
int num=0,i,d1,d2,k;
a[n]=a[0];
for(int i=0; i<n; i++)
{
if(PointOnSegment(t,a[i],a[i+1]))
return 2;
k=cmp(det(a[i+1]-a[i],t-a[i]));
d1=cmp(a[i].y=t.y);
d2=cmp(a[i+1].y-t.y);
if(k>0&&d1<=0&&d2>0)
num++;
if(k<0&&d2<=0&&d1>0)
num--;
}
return num!=0;
}//点是否在多边形内部
point MassCenter()
{
point ans=point (0,0);
if(cmp(area())==0)
return ans;
a[n]=a[0];
for(int i=0; i<n; i++)
ans=ans+(a[i]+a[i+1])*det(a[i+1],a[i]);
return ans/area()/6.0;
}//重心
int Border_Int_Point_Num()
{
int num=0;
a[n]=a[0];
for(int i=0; i<n; i++)
num+=__gcd (abs(int(a[i+1].x-a[i].x)),abs(int(a[i+1].y-a[i].y)));
return num;
}
int Inside_Int_Point_Num()
{
return int(area())+1-Border_Int_Point_Num()/2;
}//多边形内格点数
};
struct polygon_convex
{
vector <point> P;
polygon_convex(int size=0)
{
P.resize(size);
}
};
bool comp_less(const point &a,const point &b)
{
return cmp(a.x-b.x)<0||cmp(a.x-b.x)==0&&cmp(a.y-b.y)<0;
}
polygon_convex convex_hull(vector <point> a)
{
polygon_convex res(2*a.size()+5);
sort(a.begin(),a.end(),comp_less);
a.erase(unique(a.begin(),a.end()),a.end());
int m=0;
for(int i=0; i<a.size(); i++)
{
while(m>1 && cmp(det(res.P[m-1]-res.P[m-2],a[i]-res.P[m-2]))<=0)
m--;
res.P[m++]=a[i];
}
int k=m;
for(int i=int(a.size())-2; i>=0; i--)
{
while(m>k && cmp(det(res.P[m-1]-res.P[m-2],a[i]-res.P[m-2]))<=0)
m--;
res.P[m++]=a[i];
}
res.P.resize(m);
if(a.size()>1)
res.P.resize(m-1);
return res;
}//凸包
bool contain(const polygon_convex &a,const point &b)
{
int n=a.P.size();
point g=(a.P[0]+a.P[n/3]+a.P[2*n/3])/3.0;
int l=0,r=n;
//二分凸包 g-a.P[a]-a.P[b]
while(l+1<r)
{
int mid=(l+r)/2;
if(cmp(det(a.P[l]-g,a.P[mid]-g))>0)
{
if(cmp(det(a.P[l]-g,b-g))>=0&&cmp(det(a.P[mid]-g,b-g))<0)
r=mid;
else l=mid;
}
else
{
if(cmp(det(a.P[l]-g,b-g))<0&&cmp(det(a.P[mid]-g,b-g))>=0)
l=mid;
else r=mid;
}
}
r%=n;
int z=cmp(det(a.P[r]-b,a.P[l]-b))-1;
if(z==-2)
return 1;
return z;
}//O(longn)的判断点是否在凸包里
double convex_diameter (polygon_convex &a,int First,int Second)
{
vector <point> p=a.P;
int n=p.size();
double maxd=0;
if(n==1)
{
First=Second=0;
return maxd;
}
for(int i=0,j=1; i<n; ++i)
{
while(cmp(det(p[next(i)]-p[i],p[j]-p[i])-det(p[next(i)]-p[i],p[next(j)]-p[i]))<0)
{
j=next(j);
}
double d=dist(p[i],p[j]);
if(d>maxd)
{
maxd=d;
First=i;
Second=j;
}
d=dist(p[next(i)],p[next(j)]);
if(d>maxd)
{
maxd=d;
First=i;
Second=j;
}
}
return maxd;
}//凸包直径
inline bool operator < (const point &a,const point &b)
{
return a.y+eps<b.y ||fabs(a.y-b.y)<eps &&a.x+eps<b.x;
}
inline double getA(const point &a)
{
double res=atan2(a.y,a.x);
if(res<0)
res+=2*pi;
return res;
}
point p[maxm],hull[maxm];
int n;
double w[maxm],sum[maxm];
inline void getHull()//逆时针凸包
{
sort(p+1,p+n+1);
int N=0;
hull[++N]=p[1];
for(int i=2; i<=n; i++)
{
while(N>1 &&det(hull[N]-hull[N-1],p[i]-hull[N-1])<=0)
N--;
hull[++N]=p[i];
}
int bak=N;
for(int i=n-1; i>=1; i--)
{
while(N>bak&&det(hull[N]-hull[N-1],p[i]-hull[N-1])<=0)
N--;
hull[++N]=p[i];
}
n=N-1;
for(int i=1; i<=n; i++)
p[i+n]=p[i]=hull[i];
p[n+n+1]=p[1];
for(int i=1; i<=n; i++)
w[i+n]=w[i]=getA(p[i+1]-p[i]);
sum[0]=0;
for(int i=1; i<=n*2; i++)
sum[i]=sum[i-1]+det(p[i],p[i+1]);//预处理有向面积前缀和
}
inline int Find(double x)
{
if(x<=w[1]||x>=w[n])
return 1;
return (upper_bound(w+1,w+n+1,x)-(w+1))+1;
}//找到第一个角度>x的边
point P,Q;
inline int getInter(int l,int r)
{
int sign;
if(det(Q-P,p[l]-P)<0)
sign=-1;
else sign=1;
while(l+1<r)
{
int mid=(l+r)/2;
if(det(Q-P,p[mid]-P)*sign>0)
l=mid;
else r=mid;
}
return r;
}//找到第一个和p[l]不在P->Q向量的同侧的点
inline point Intersect (const point &a,const point &b,const point &c,const point &d)
{
double s1=det(c-a,b-a);
double s2=det(d-a,b-a);
return (c*s2-d*s1)/(s2-s1);
}
inline bool solve(point P,point Q)
{
int i=Find (getA(Q-P));
int j=Find (getA(P-Q));
//两侧各找一点
if(det(Q-P,p[i]-P)*det(Q-P,p[j]-P)>=0)//无交点
return false;
return true;
}
struct halfPlane
{
double a,b,c;//ax+by+c<=0
halfPlane(point p,point q)
{
a=p.y-q.y;
b=q.x-p.x;
c=det(p,q);
}
halfPlane(double aa,double bb,double cc)
{
a=aa;
b=bb;
c=cc;
}
};
//计算点a和b连线与半平面L的交点
double calc(halfPlane &L,point &a)
{
return a.x*L.a+a.y*L.b+L.c;
}
point Intersect(point &a,point &b,halfPlane &L)
{
point res;
double t1=calc(L,a),t2=calc(L,b);
res.x=(t2*a.x-t1*b.x)/(t2-t1);
res.y=(t2*a.y-t1*b.y)/(t2-t1);
return res;
}
//将一个凸多边形和一个半平面交
polygon_convex cut(polygon_convex &a,halfPlane &L)
{
int n=a.P.size();
polygon_convex res;
for(int i=0; i<n; i++)
{
if(calc(L,a.P[i])<-eps)
res.P.push_back(a.P[i]);
else
{
int j;
j=i-1;
if(j<0)
j=n-1;
if(calc(L,a.P[j])<-eps)
res.P.push_back(Intersect(a.P[j],a.P[i],L));
j=i+1;
if(j==n)
j=0;
if(calc(L,a.P[j])<-eps)
res.P.push_back(Intersect(a.P[i],a.P[j],L));
}
}
return res;
}
polygon_convex core(polygon &a)
{
polygon_convex res;
res.P.push_back(point (-INF,-INF));
res.P.push_back(point(INF,-INF));
res.P.push_back(point(INF,INF));
res.P.push_back(point(-INF,INF));
int n=a.n;
for(int i=0;i<n;i++)
{
halfPlane L(a.a[i],a.a[(i+1)%n]);
res=cut(res,L);
}
return res;
}
typedef complex < double > Point;
typedef pair < Point , Point > Halfplane;
typedef vector< Point > Convex;
inline double dist(Halfplane hal)
{
return sqrt((hal.second.imag()-hal.first.imag())*(hal.second.imag()-hal.first.imag())+(hal.second.real()-hal.first.real())*(hal.second.real()-hal.first.real()));
}
inline int sgn(double n)
{
return fabs(n)<eps?0:(n<0?-1:1);
}
inline double cross(Point a,Point b)
{
return (conj(a)*b).imag();
}
inline double dot(Point a,Point b)
{
return(conj(a)*b).real();
}
inline double satisfy(Point a,Halfplane p)
{
return sgn(cross(a-p.first,p.second-p.first))<=0;
}
Point crosspoint (const Halfplane &a,const Halfplane &b)
{
double k=cross(b.first-b.second,a.first-b.second);
k=k/(k-cross(b.first-b.second,a.second-b.second));
return a.first+(a.second-a.first)*k;
}
bool cmp2(const Halfplane &a,const Halfplane &b)
{
int res= sgn(arg(a.second-a.first)-arg(b.second-b.first));
return res==0?satisfy(a.first,b):res<0;
}
vector < Point > halfplaneIntersection(vector<Halfplane>v)
{
sort(v.begin(),v.end(),cmp2);
deque < Halfplane > q;
deque < Point > ans;
q.push_back(v[0]);
for(int i=1; i<int(v.size()); i++)
{
if(sgn(arg(v[i].second - v[i].first)-arg(v[i-1].second - v[i-1].first))==0)
continue;
while(ans.size()>0 && !satisfy(ans.back(),v[i]))
{
ans.pop_back();
q.pop_back();
}
while(ans.size()>0 && !satisfy(ans.front(),v[i]))
{
ans.pop_front();
q.pop_front();
}
ans.push_back(crosspoint (q.back(),v[i]));
q.push_back(v[i]);
}
while(ans.size() > 0 && !satisfy(ans.back(),q.front()))
{
ans.pop_back();
q.pop_back();
}
while(ans.size() > 0 && !satisfy(ans.front(),q.back()))
{
ans.pop_front();
q.pop_front();
}
ans.push_back(crosspoint(q.back(),q.front()));
return vector < Point > (ans.begin(),ans.end());
}
Convex convexIntersection(Convex v1,Convex v2)
{
vector <Halfplane> h;
for(int i=0; i<int(v1.size()); i++)
h.push_back(Halfplane(v1[i],v1[(i+1) % v1.size()]));
for(int i=0; i<int(v2.size()); i++)
h.push_back(Halfplane(v2[i],v2[(i+1) % v2.size()]));
return halfplaneIntersection(h);
}//凸多边形交
int dcmp(double k)
{
return k<-eps?-1:k>eps?1:0;
}
double abs(const point &o)
{
return sqrt(dot(o,o));
}
point crosspt(const point &a,const point &b,const point &p,const point &q)
{
double a1=det(b-a,p-a);
double a2=det(b-a,q-a);
return (p*a2-q*a1)/(a2-a1);
}
double mysqrt(double n)
{
return sqrt(max(0.0,n));
}
void circle_cross_line(point a,point b,point o,double r,point ret[],int &num)
{
double x0=o.x,y0=o.y;
double x1=a.x,y1=a.y;
double x2=b.x,y2=b.y;
double dx=x2-x1,dy=y2-y1;
double A=dx*dx+dy*dy;
double B=2*dx*(x1-x0) +2*dy*(y1-y0);
double C=sqr(x1-x0)+sqr(y1-y0)-sqr(r);
double delta =B*B-4*A*C;
num=0;
if(dcmp(delta) >= 0)
{
double t1=(-B-mysqrt(delta))/(2*A);
double t2=(-B+mysqrt(delta))/(2*A);
if(dcmp(t1-1)<=0&&dcmp(t1)>=0)//如果是直线,去掉这行if
ret[num++]=point(x1+t1*dx,y1+t1*dy);
if(dcmp(t2-1)<=0&&dcmp(t2)>=0)//如果是直线,去掉这行if
ret[num++]=point(x1+t2*dx,y1+t2*dy);
}
}//圆和直线交,ab(线段或直线),o圆心,r半径,ret交点,num交点数
point res[maxn];
double r;
double sector_area(const point &a , const point &b)
{
double theta =atan2(a.y,a.x) -atan2(b.y,b.x);
while(theta<=0)theta+=2*pi;
while(theta>2*pi)theta-=2*pi;
theta=min(theta,2*pi-theta);
return r*r*theta/2;
}
double calc1(const point &a,const point &b)
{
point p[2];
int num=0;
int ina=dcmp(abs(a)-r)<0;
int inb=dcmp(abs(b)-r)<0;
if(ina)
{
if(inb)
return fabs(det(a,b))/2.0;
else
{
circle_cross_line(a,b,point(0,0),r,p,num);
return sector_area(b,p[0])+fabs(det(a,p[0]))/2.0;
}
}
else
{
if(inb)
{
circle_cross_line(a,b,point(0,0),r,p,num);
return sector_area(p[0],a)+fabs(det(p[0],b))/2.0;
}
else
{
circle_cross_line(a,b,point(0,0),r,p,num);
if(num==2)
{
return sector_area(a,p[0])+sector_area(p[1],b)+fabs(det(p[0],p[1]))/2.0;
}
else return sector_area(a,b);
}
}
}
double area()
{
double ret=0;
for(int i=0;i<n;i++)
{
int sgn=dcmp(det(res[i],res[i+1]));
if(sgn!=0)
ret+=sgn*calc1(res[i],res[i+1]);
}
return ret;
}//res[n]=res[0],n全局变量,多边形点数,res,全局变量,逆时针存入多边形的所有点,r,全局变量,圆的半径
几何模板
最新推荐文章于 2020-09-26 17:39:35 发布