数学库

const   double   INF=1e100;  
  const   double   ZERO=1e-6;  
  const   double   PI=2*asin(1.0);  
   
  struct   XYpoint{     //(x,y)  
        double   x;  
        double   y;  
  };  
   
  struct   XYline{       //   Ax+By+C=0;  
        double   A;  
        double   B;  
        double   C;  
  };  
   
  struct   XYsegment{  
        XYpoint   a,b;  
  };  
   
  struct   XYround{     //   圆  
        XYpoint   center;  
        double   r;  
  };  
   
  inline   double   min(double   a,double   b)  
  {  
        return   a<b?a:b;  
  }  
  inline   double   max(double   a,double   b)  
  {  
        return   a>b?a:b;  
  }  
   
   
  /********************************************/  
  /*                         两点确定一条直线                             */  
  /********************************************/  
  XYline   makeLine(double   x1,double   y1,double   x2,double   y2)  
  {  
        XYline   line;  
        line.A=(y2-y1);  
        line.B=(x1-x2);  
        line.C=y1*(x2-x1)+x1*(y1-y2);  
        return   line;  
  }  
   
  /********************************************/  
  /*                               两点间距离                                   */  
  /********************************************/  
  inline   double   dis_PP(double   x1,double   y1,double   x2,double   y2)  
  {  
        return   sqrt(   (x1-x2)*(x1-x2)   +   (y1-y2)*(y1-y2)   );  
  }  
   
  /********************************************/  
  /*                               点直线距离                                   */  
  /********************************************/  
  double   dis_PL(double   x,double   y,XYline   line)  
  {  
        return   fabs(line.A*x   +   line.B*y   +   line.C)   /   sqrt(line.A*line.A   +   line.B*line.B);  
  }  
   
  /********************************************/  
  /*                     海伦公式求三角形面积                         */  
  /********************************************/  
  double   triangleArea(double   x1,double   y1,double   x2,double   y2,double   x3,double   y3)  
  {  
        double   a=sqrt(   (x1-x2)*(x1-x2)   +   (y1-y2)*(y1-y2)   );  
        double   b=sqrt(   (x1-x3)*(x1-x3)   +   (y1-y3)*(y1-y3)   );  
        double   c=sqrt(   (x3-x2)*(x3-x2)   +   (y3-y2)*(y3-y2)   );  
        double   s=(a+b+c)/2;  
        return   sqrt(   s   *   (s-a)   *   (s-b)   *   (s-c)   );  
  }  
  inline   double   tArea(XYpoint   p1,XYpoint   p2,XYpoint   p3)  
  {  
          return   fabs(0.5*((p3.y-p1.y)*(p2.x-p1.x)-(p2.y-p1.y)*(p3.x-p1.x)));  
  }  
   
  /********************************************/  
  /*                             判断直线相交                                 */  
  /*         1   -   有交点     0   -   无交点     -1   -   重合           */  
  /********************************************/  
  int   inter_LL(   XYline   line1   ,   XYline   line2   ,   XYpoint   &point   )  
  {  
        if(line1.A*line2.B   ==   line1.B*line2.A)  
        {  
              if(line1.A*line2.C   ==   line1.C*line2.A   &&   line1.B*line2.C   ==   line1.C*line2.B)  
                    return   -1;//   重合  
              else  
                    return   0;//   平行  
        }  
        else  
        {  
              point.x=   -1   *   (line1.C*line2.B-line2.C*line1.B)   /   (line1.A*line2.B-line2.A*line1.B);  
              point.y=   -1   *   (line1.C*line2.A-line2.C*line1.A)   /   (line1.B*line2.A-line2.B*line1.A);  
              //   交点  
              return   1;  
        }  
  }  
   
  /********************************************/  
  /*                         两直线夹角   （弧度）                       */  
  /********************************************/  
  double   angle_LL(XYline   line1   ,   XYline   line2)  
  {  
        double   t;  
        if(   (line1.A*line2.A   +   line1.B*line2.B)==0   )  
              t=INF;  
        else  
              t=fabs(   (line1.A*line2.B   -   line2.A*line1.B)   /   (line1.A*line2.A   +   line1.B*line2.B)   );  
        return   atan(t);  
  }  
   
  /********************************************/  
  /*                     过某点和直线垂直的直线                     */  
  /********************************************/  
  XYline   ver_PL(   XYpoint   p   ,   XYline   line   )  
  {  
        XYline   temp;  
        temp.A=line.B;  
        temp.B=-1   *   line.A;  
        temp.C=-temp.A*p.x-temp.B*p.y;  
        return   temp;  
  }  
   
  /********************************************/  
  /*                       多边形面积   (   v[0]=v[n]   )               */  
  /*                     n                                                               */  
  /*     a=0.5*sigma{   X[i-1]*Y[i]-X[i]*Y[i-1]   }     */  
  /*                   i=1                                                             */  
  /********************************************/  
  double   polygonArea(   int   vcount   ,   XYpoint   ver[]   )  
  {  
        int   i;  
        double   s=0;  
        XYpoint   t;  
        if(   vcount<3   )  
              return   0;  
        t=ver[vcount];  
        ver[vcount]=ver[0];  
        for(i=1;i<=vcount;i++)  
              s+=ver[i-1].x*ver[i].y   -   ver[i].x*ver[i-1].y   ;  
        ver[vcount]=t;  
        return   s/2   ;  
  }  
   
  /********************************************/  
  /*   返回(P1-P0)*(P2-P0)的叉积。                             */  
  /*   若为正，则<P0,P1>在<P0,P2>的顺时针方向       */  
  /*   若为零，则<P0,P1><P0,P2>共线；                       */  
  /*   若为负，则<P0,P1>在<P0,P2>的在逆时针方向   */  
  /*   这个函数可以确定两条线段在交点处的转向       */  
  /*   比如确定p0p1和p1p2在p1处是左转还是右转，   */  
  /*   只要求   (p2-p0)*(p1-p0)，                                   */  
  /*   若<0则左转，>0则右转，=0则共线                       */  
  /********************************************/  
  double   multiply(   XYpoint   p1   ,   XYpoint   p2   ,   XYpoint   p0   )  
  {  
          /*  
          double   x0=p0.x;  
          double   x1=p1.x;  
          double   y0=p0.y;  
          double   y0=p1.y;                          
          double   r;  
          r=((x1-x0)*(y2-y0)-(x2-x0)*(y1-y0));  
          if(   r<0   )   return   -1;  
          if(   r==0   )   return   0;  
          if(   r>0   )   return   1;  
          */  
  return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));  
   
  }

/********************************************/  
  /*                           判断直线段相交                               */  
  /********************************************/  
  bool   inter_SS(   XYsegment   u   ,   XYsegment   v   )  
  {  
  return(   (max(u.a.x,u.b.x)>=min(v.a.x,v.b.x))&&  
                  (max(v.a.x,v.b.x)>=min(u.a.x,u.b.x))&&  
                  (max(u.a.y,u.b.y)>=min(v.a.y,v.b.y))&&  
                  (max(v.a.y,v.b.y)>=min(u.a.y,u.b.y))&&  
                  (multiply(v.a,u.b,u.a)*multiply(u.b,v.b,u.a)>=0)&&  
                  (multiply(u.a,v.b,v.a)*multiply(v.b,u.b,v.a)>=0));  
  }  
   
   
  /********************************************/  
  /*                           判断两点是否共点                           */  
  /********************************************/  
  bool   equal_PP(XYpoint   a,XYpoint   b)  
  {  
        if(   fabs(a.x-b.x)<ZERO   &&   fabs(a.y-b.y)<ZERO   )  
              return   true;  
        else  
              return   false;  
  }  
   
  /********************************************/  
  /*                       判断点p是否在线段l上                       */  
  /********************************************/  
  bool   online(   XYsegment   l   ,   XYpoint   p   )  
  {  
        if(   equal_PP(p,l.a)   ||   equal_PP(p,l.b)   )  
              return   true;  
  return(   (multiply(l.b,p,l.a)==0)&&(   ((p.x-l.a.x)*(p.x-l.b.x)<0   )||(   (p.y-l.a.y)*(p.y-l.b.y)<0   ))   );  
  }  
   
  /********************************************/  
  /*                     判断点是否在多边形内                         */  
  /********************************************/  
  bool   insidePolygon(   int   vcount   ,   XYpoint   ver[]   ,   XYpoint   point   )  
  {  
        int   i,c=0;  
        XYpoint   t;  
        XYsegment   line,edge;  
        t=ver[vcount];  
        ver[vcount]=ver[0];  
        line.a=line.b=point;  
        line.b.x=INF;  
        for(i=0;i<vcount;i++)  
        {  
              edge.a=ver[i];   edge.b=ver[i+1];  
              if(   online(edge,point)   )  
                    return   true;  
              if(   ver[i].y==ver[i+1].y   )  
                    continue;  
              if(   min(ver[i].y,ver[i].y)!=point.y   &&   inter_SS(line,edge)   )  
                    c++;  
        }  
        ver[vcount]=t;  
        return   c%2   ;  
  }  
   
  /********************************************/  
  /*                           得到线段的中点                               */  
  /********************************************/  
  XYpoint   middle(   XYpoint   p1   ,   XYpoint   p2   )  
  {  
        XYpoint   point;  
        point.x=   (p1.x+p2.x)/2;  
        point.y=   (p1.y+p2.y)/2;  
        return   point;  
  }  
   
  /********************************************/  
  /*                     判断线段是否在多边形内                     */  
  /*         1   -   在内部     0   -   在外部     -1   -   相交           */  
  /********************************************/  
  int   segmentPolygon(   int   vcount   ,   XYpoint   ver[]   ,   XYsegment   seg   )  
  {  
        int   i;  
        XYpoint   point,t;  
        XYsegment   edge;  
        t=ver[vcount];  
        ver[vcount]=ver[0];  
        for(i=0;i<vcount;i++)  
        {  
              edge.a=ver[i];   edge.b=ver[i+1];  
              if(   inter_SS(seg,edge)   )//&&   !online(edge,seg.a)   &&   !online(edge,seg.b)   )  
                    return   -1;  
        }  
        return   insidePolygon(   vcount   ,   ver   ,   middle(seg.a,seg.b)   )   ;  
  }  
   
  /********************************************/  
  /*             计算线段（p0->p1)的幅角     [0,2*PI)       */  
  /********************************************/  
  double   angle_PP(   XYpoint   p0   ,   XYpoint   p1   )  
  {  
        double   x,y;  
        x=p1.x-p0.x;  
        y=p1.y-p0.y;  
        if(   equal_PP(   p0   ,   p1   )   )  
              return   -1;  
        if(   x==0   )   {  
              if(   y>0   )  
                    return   PI/2;  
              else  
                    return   PI+PI/2;  
        }  
        if(   x>0   &&   y>=0   )   return   atan(y/x);  
        if(   x>0   &&   y<0   )   return   2*PI+atan(y/x);  
        if(   x<0   &&   y>=0   )   return   PI+atan(y/x);  
        if(   x<0   &&   y<0   )   return   PI+atan(y/x);  
  }  
   
  /********************************************/  
  /*                             卷包裹法寻找凸包                         */  
  /********************************************/  
  int   convexHull(   int   pointNum   ,   XYpoint   pointSet[]   ,   XYpoint   convex[]   )  
  {  
        int   i;  
        int   verNum;  
        int   lastPoint,nowPoint;  
        double   lastAngle,nowAngle,tempAngle,t1,t2;  
        //   找到   最右   最下的点  
        for(i=0,lastPoint=0;i<pointNum;i++)  
        {  
              if(   pointSet[i].y<pointSet[lastPoint].y   )  
                    lastPoint=i;  
              else   if(   pointSet[i].y==pointSet[lastPoint].y   &&  
                                pointSet[i].x<pointSet[lastPoint].x   )  
                    lastPoint=i;  
        }  
        //   寻找凸包  
        lastAngle=PI;  
        convex[0]=pointSet[lastPoint];  
        verNum=1;  
        do   {  
              nowAngle=lastAngle-2*PI;  
              for(i=0;i<pointNum;i++)  
              {  
                    tempAngle=angle_PP(   convex[verNum-1]   ,   pointSet[i]   );  
                    if(   tempAngle!=-1   )  
                    {  
                          t1=lastAngle-tempAngle;   t2=lastAngle-nowAngle;  
                          if(t1<0)   t1+=2*PI;  
                          if(t2<0)   t2+=2*PI;  
                          if(   t1<t2   )   {  
                                nowPoint=i;   nowAngle=tempAngle;  
                          }  
                          else   if(   t1==t2   )   {  
                                if(   dis_PP(convex[verNum-1].x,convex[verNum-1].y,pointSet[i].x,pointSet[i].y)  
                                      <dis_PP(convex[verNum-1].x,convex[verNum-1].y,pointSet[nowPoint].x,pointSet[nowPoint].y)   )  
                                        nowPoint=i;   nowAngle=tempAngle;  
                          }  
                    }  
              }  
              lastAngle=nowAngle;  
              convex[verNum]=pointSet[nowPoint];  
              verNum++;  
        }  
        while(   !equal_PP(convex[0],convex[verNum-1])   )   ;  
        return   verNum-1;  
  }  
   
  /********************************************/  
  /*                             两点确定一个圆                             */  
  /********************************************/  
  XYround   makeRound2(XYpoint   p1,XYpoint   p2)  
  {  
          XYround   r;  
          r.center=middle(p1,p2);  
          r.r=dis_PP(p1.x,p1.y,r.center.x,r.center.y);  
          return   r;  
  }  
   
  /********************************************/  
  /*                             三点确定一个圆                             */  
  /********************************************/  
  XYround   makeRound3(XYpoint   p1,XYpoint   p2,XYpoint   p3)  
  {  
          XYround   r;  
          XYline   l1,l2;  
          if(   fabs(   (p1.x-p2.x)/(p1.y-p2.y)-(p1.x-p3.x)/(p1.y-p3.y)   )<ZERO   )  
          {  
                  r.r=INF;  
                  return   r;  
          }  
          l1=ver_PL(   middle(p1,p2),   makeLine(p1.x,p1.y,p2.x,p2.y)   );  
          l2=ver_PL(   middle(p1,p3),   makeLine(p1.x,p1.y,p3.x,p3.y)   );        
          inter_LL(l1,l2,r.center);  
          r.r=dis_PP(p1.x,p1.y,r.center.x,r.center.y);  
          return   r;  
  }  
   
  /********************************************/  
  /*                             判断点是否在圆内                         */  
  /********************************************/  
  inline   bool   insideRound(   XYround   r,XYpoint   p   )  
  {  
          return   (   dis_PP(p.x,p.y,r.center.x,r.center.y)<=r.r   );  
  }