几何模板

花了些时间准备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);
 }