几何模板

花了些时间准备acm的比赛，现在比赛结束有一段时间了，忙着上自习、做实验。有空就看看linux的书，至于算法，很久没有碰了。昨天刚想回toj做几个题，想不到传来toj系统维护的消息。

回顾比赛，计算几何的东西基本没有用到，这个不过这个模板上的线段树还是发挥了点作用。把模板传上来，算是一个阶段性的总结。也希望大牛多给我指出些不对的地方。

/*

*常量定义

*/

const double INF=1e200;

const double EPS=1e-7;

const double PI =acos(-1);

/*

*点的定义

*/

struct PT

{

double x,y;

PT(){}

PT(double xx,double yy)//构造函数

{

x=xx,y=yy;

return;

}

};

/*

*线段的定义

*/

struct SEG

{

PT st,end;

int len;

};

/*

*圆的定义

*/

struct CIR

{

double x,y,r;

CIR(){}

CIR(double xx,double yy,double rr)

{

x=xx,y=yy,r=rr;

return;

}

};

/*

*符号函数

*/

inline int sign(double a)

{

if(fabs(a)<EPS)

return 0;

return a>EPS?1:-1;

}

/*

*叉积

*/

inline int crossP(PT &a,PT &b,PT &c)

//cross product of ab and ac

{

return (c.x-a.x)*(b.y-a.y)-(b.x-a.x)*(c.y-a.y);

}

inline double crossP2(const PT &a,const PT &b,const PT &c,const PT &d)

//cross product of ab and cd

{

return (b.x-a.x)*(d.y-c.y)-(b.y-a.y)*(d.x-c.x);

}

/*

*点积

*/

inline double dotP(PT &a,PT &b,PT &c)

{

return (c.x-a.x)*(b.x-a.x)+(c.y-a.y)*(b.y-a.y);

}

/*

*两点间距离

*/

inline int dis(PT a,PT b)

{

return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));

}

/*

*点到线段的距离

*/

double dis2seg(PT pt,PT st,PT end)

{

if(dotP(st,end,pt)*dotP(end,st,pt)<0)

return min(dis(pt,st),dis(pt,end));

return fabs(crossP(st,end,pt)/dis(st,end));

}

/*

*已知圆上三点求该圆

*pt1,pt2,pt3是圆上三点

*cen是圆心的引用，返回半径

*/

double getRadiusBy3PT(PT pt1,PT pt2,PT pt3,PT &cen)

{

double a1,a2,b1,b2,c1,c2;

a1=2*(pt2.x-pt1.x);

a2=2*(pt3.x-pt1.x);

b1=2*(pt2.y-pt1.y);

b2=2*(pt3.y-pt1.y);

c1=pt2.x*pt2.x+pt2.y*pt2.y-pt1.x*pt1.x-pt1.y*pt1.y;

c2=pt3.x*pt3.x+pt3.y*pt3.y-pt1.x*pt1.x-pt1.y*pt1.y;

cen.y=(a2*c1-a1*c2)/(a2*b1-a1*b2);

cen.x=(b2*c1-b1*c2)/(b2*a1-b1*a2);

return dis(pt1,c);

}

/*

*两圆相交部分的面积

*不相交返回0

*/

double CircleIntersectionArea(CIR a,CIR b)

{

double len=dis(a,b);

if(sign(len-a.r-b.r)>=0 || a.r<EPS || b.r< EPS)

return 0;

else if(sign(a.r+len-b.r)<=0)

return PI*a.r*a.r;

else if(sign(b.r+len-a.r)<=0)

return PI*b.r*b.r;

double ang1=2*acos((a.r*a.r+len*len-b.r*b.r)/(2*a.r*len));

double ang2=2*acos((b.r*b.r+len*len-a.r*a.r)/(2*b.r*len));

double re=a.r*a.r*ang1*0.5-0.5*a.r*a.r*sin(ang1);

re+=b.r*b.r*ang2*0.5-0.5*b.r*b.r*sin(ang2);

return re;

}

/*

*判断线段是否相交

*端点相交返回false

*如果需要，把<改为<=

*那么，共线也会返回true

*/

inline bool intersect(PT s1,PT e1,PT s2,PT e2)

{

return((max(s1.x,e1.x)>=min(s2.x,e2.x)) &&

(max(s2.x,e2.x)>=min(s1.x,e1.x)) &&

(max(s1.y,e1.y)>=min(s2.y,e2.y)) &&

(max(s2.y,e2.y)>=min(s1.y,e1.y)) &&

(crossP(s1,e1,s2)*crossP(s1,e1,e2)<0) &&

(crossP(s2,e2,s1)*crossP(s2,e2,e1)<0));

}

/*

判断线段是否严格相交

*/

int strictcross(point p1, point p2, point p3, point p4) {

T a, b, c, d;

a=(p3.y-p2.y)*(p1.x-p2.x)-(p1.y-p2.y)*(p3.x-p2.x);

b=(p4.y-p2.y)*(p1.x-p2.x)-(p1.y-p2.y)*(p4.x-p2.x);

c=(p1.y-p4.y)*(p3.x-p4.x)-(p3.y-p4.y)*(p1.x-p4.x);

d=(p2.y-p4.y)*(p3.x-p4.x)-(p3.y-p4.y)*(p2.x-p4.x);

return a*b<0 && c*d<0;

}

/*

*判断简单多边形是否是凸多边形

*vpt[]:所有点

*n:点的数量

*/

bool isprotude(PT vpt[],int n)

{

int t,tt,pre;

for(t=0;t<n;t++)

{

tt=0;

while(pre=sign(crossP(vpt[t],vpt[(t+1)%n],vpt[(t+2+tt)%n])),pre==0)

tt++;

for(;tt<n-2;tt++)

{

double d=crossP(vpt[t],vpt[(t+1)%n],vpt[(t+2+tt)%n]);

if(sign(d)==0)

continue;

if(sign(d)!=pre)

return false;

}

}

return true;

}

/*

*判断点是否在凸多边形内部

*在边上为false

*for convex polygen only!

*/

bool inPoly(PT vpt[],int n,PT &pt)

{

int t;

PT st,end;

int dr=sign(crossP(vpt[0],vpt[1],pt));

for(t=1;t<n;t++)

{

st=vpt[t];

end=vpt[(t+1)%n];

if(sign(crossP(st,end,pt))!=dr)

return false;

}

return true;

}

/*

*Graham法求凸包

*双链法扫描

*pt[]:所有点

*n:pt[]中点的个数

*ans[]:用于返回凸包中的点

*返回凸包中点的个数

*默认逆时针时针方向

*/

int graham(int n,PT pt[],PT ans[])

{

int stk[MAX];

bool vis[MAX];

int t,p,ps,pans=0;

sort(pt,pt+n,cmp);

memset(vis,false,sizeof(vis));

stk[0]=0;

stk[1]=1;

vis[1]=true;

ps=p=2;

while(p<n)

{

if(ps<2 || crossP(pt[stk[ps-2]],pt[stk[ps-1]],pt[p])>=0)

//or crossP()<=0 to be clockwise

//if points in stk is not enough

//or pt[p] is on the left side or on segmemt

//push pt[p] into stack

//and mark pt[p] as visited

{

stk[ps++]=p++;

vis[p]=true;

}

else

//pop the top point

//and mark it as unvisited

vis[stk[--ps]]=false;

}

for(t=0;t<ps;t++)

ans[pans++]=pt[stk[t]];

vis[n-1]=false;//second travle start here

ps=0;//clear the stack to start new travle

for(p=n-1;p>=0;p--)

if(!vis[p])

{

stk[ps++]=p;

if(ps==2)//push first 2 points

break;

}

p--;//next point

while(p>=0)//mark is not needed any more

{

if(vis[p])

{

p--;

continue;

}

if(ps<2 || crossP(pt[stk[ps-2]],pt[stk[ps-1]],pt[p])>=0)

//or crossP()<=0 to be clockwise

stk[ps++]=p--;

else

--ps;

}

for(t=1;t<ps-1;t++)

ans[pans++]=pt[stk[t]];

return pans;

}

/*

旋转卡壳求最远点距离

*/

double RotatingCaliper(const int &n,PT pt[])

{

double re=0;

int t=0,tt=1;

while(sign(crossP2(pt[t],pt[t+1],pt[tt],pt[tt+1]))==1)

tt++;

while);<=n)

{

int tmp=sign(crossP2(pt[t%n],pt[(t+1)%n],pt[tt%n],pt[(tt+1)%n]));

if(tmp==1)

tt=(tt+1)%n;

else if(tmp==-1)

t++;

else

{

t++;

continue;

}

double d=dis(pt[t],pt[tt]);

if(d>re)

re=d;

}

return re;

}

/*

*线段树，离散化

*/

struct Line

{

double x,tp,bt;

bool  left;

}line[2005];

double y[2005];

int num,ny;//离散化前后y[]的大小

inline int find(double a)//二分查找对应的离散值

{

int low=0,high=ny-1,mid;

while(low<=high)

{

mid=(low+high)/2;

if(a==y[mid])

return mid;

else if(a<y[mid])

high=mid;

else

low=mid+1;

}

return -1;

}

struct node

{

int st,end,c;//c是区间被覆盖次数

double m;//m是测度

node *ls,*rs;

}* root;

node * build(int st,int end)

{

node *p=new node;

p->st=st;

p->end=end;

p->c=0;

p->m=0;

if(st+1<end)//注意一下这里

{

int mid=(st+end)/2;

p->ls=build(st,mid);

p->rs=build(mid,end);

}

else

p->ls=p->rs=NULL;

return p;

}

void init(node *p)//初始化

{

p->c=0;

p->m=0;

if(p->ls!=NULL)

init(p->ls);

if(p->rs!=NULL)

init(p->rs);

return;

}

inline void update(node *p)//更新测度

{

if(p->c>0)

p->m=y[p->end]-y[p->st];

else if(p->st+1==p->end)

p->m=0;

else

p->m=p->ls->m+p->rs->m;

return;

}

void ins(node *p,int st,int end)

{

if(st<=p->st && p->end<=end)

{

p->c++;

update(p);

return;

}

int mid=(p->st + p->end)/2;

if(p->ls!=NULL)

{

ins(p->ls,st,end);

update(p);

}

if(p->rs!=NULL)

{

ins(p->rs,st,end);

update(p);

}

return;

}

void del(node *p,int st,int end)

{

if(st<=p->st && p->end<=end)

{

p->c--;

update(p);

return;

}

if(p->ls!=NULL)

{

del(p->ls,st,end);

update(p);

}

if(p->rs!=NULL )

{

del(p->rs,st,end);

update(p);

}

return;

}

/*

*树状数组，有待完善

*/

const int MAX=32005*3;//大一点好

struct node

{

int st,end,c;

}T[MAX];

inline int ls(int k)

{

return 2*k;

}

inline int rs(int k)

{

return 2*k+1;

}

void build(int p,int st,int end)//建树

{

T[p].st=st;

T[p].end=end;

T[p].c=0;

if(st<end)

{

int mid=(st+end)/2;

build(ls(p),st,mid);

build(rs(p),mid+1,end);

}

return;

}

void ins(int p,int k)//插入一个点

{

if(T[p].st<=k && k<=T[p].end)

{

T[p].c++;

if(T[p].st<T[p].end)

{

ins(ls(p),k);

ins(rs(p),k);

}

}

return;

}

int ask(int p,int st,int end)//询问一个区间内有几个点

{

if(st<=T[p].st && T[p].end<=end)

return T[p].c;

else if(st>T[p].end || end<T[p].st)

return 0;

else

return ask(ls(p),st,end)+ask(rs(p),st,end);

}