计算几何学库函数.txt

解决排列组合问题的通用算法,自然数得拆分,求最大公约数的算法

/********************************************* *                                          *
*        计算几何学库函数                  *
*                                          *
*                  copyright starfish      *
*                          2000/10/24      *
*                                          *
\*********************************************/

#include <stdlib.h>


#define infinity 1e20
#define EP 1e-10

struct TPoint {
    float x,y;
    };

struct TLineSeg {
    TPoint a,b;
    };



//求平面上两点之间的距离
float distance(TPoint p1,TPoint p2)
{
    return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));  
}


/******************************************************** *                                                      *
*  返回(P1-P0)*(P2-P0)的叉积。                        *
*  若结果为正，则<P0,P1>在<P0,P2>的顺时针方向；        *
*  若为0则<P0,P1><P0,P2>共线；                        *
*  若为负则<P0,P1>在<P0,P2>的在逆时针方向;            *
*  可以根据这个函数确定两条线段在交点处的转向,        *
*  比如确定p0p1和p1p2在p1处是左转还是右转，只要求      *
*  (p2-p0)*(p1-p0)，若<0则左转，>0则右转，=0则共线    *
*                                                      *
\********************************************************/

float multiply(TPoint p1,TPoint p2,TPoint p0)
{
    return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));    
}



//确定两条线段是否相交
int intersect(TLineSeg u,TLineSeg 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));
}

//判断点p是否在线段l上
int online(TLineSeg l,TPoint p)
{
    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 )) );
}


//判断两个点是否相等
int Euqal_Point(TPoint p1,TPoint p2)
{
return((abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP));
}

//一种线段相交判断函数，当且仅当u,v相交并且交点不是u,v的端点时函数为true;
int intersect_A(TLineSeg u,TLineSeg v)
{
    return((intersect(u,v))&&
          (!Euqal_Point(u.a,v.a))&&
          (!Euqal_Point(u.a,v.b))&&
          (!Euqal_Point(u.b,v.a))&&
          (!Euqal_Point(u.b,v.b)));
}


/************************************************* *                                              *
*      判断点q是否在多边形Polygon内，          *
*      其中多边形是任意的凸或凹多边形，        *
*      Polygon中存放多边形的逆时针顶点序列      *
*                                              *
\*************************************************/

int InsidePolygon(int vcount,TPoint Polygon[],TPoint q)
{
    int c=0,i,n;
    TLineSeg l1,l2;
        
    l1.a=q;
    l1.b=q;
    l1.b.x=infinity;
    n=vcount;
    for (i=0;i<vcount;i++)
    {
        l2.a=Polygon[i];
        l2.b=Polygon[(i+1)%n];
        if (  (intersect_A(l1,l2))||
              (
                (online(l1,Polygon[(i+1)%n]))&&
                (
                (!online(l1,Polygon[(i+2)%n]))&&
                (multiply(Polygon[i],Polygon[(i+1)%n],l1.a)*multiply(Polygon[(i+1)%n],Polygon[(i+2)%n],l1.a)>0)
                ||
                (online(l1,Polygon[(i+2)%n]))&&
                (multiply(Polygon[i],Polygon[(i+2)%n],l1.a)*multiply(Polygon[(i+2)%n],Polygon[(i+3)%n],l1.a)>0)        
                )
              )
          ) c++;
        }
        return(c%2!=0);
}


/******************************************** *                                          *
*      计算多边形的面积                    *
*                                          *
*    要求按照逆时针方向输入多边形顶点    *
*    可以是凸多边形或凹多边形            *
*                                          *
\********************************************/

float area_of_polygon(int vcount,float x[],float y[])
{
  int i;
  float s;
  if (vcount<3) return 0;
  s=y[0]*(x[vcount-1]-x[1]);
  for (i=1;i<vcount;i++)
    s+=y[i]*(x[(i-1)]-x[(i+1)%vcount]);
  return s/2;
  }


/*====================================================================================
  寻找凸包 graham 扫描法

  PointSet为输入的点集；
  ch为输出的凸包上的点集，按照逆时针方向排列;
  n为PointSet中的点的数目
  len为输出的凸包上的点的个数；

      设下一个扫描的点PointSet[i]=P2,当前栈顶的两个点ch[top]=P1,ch[top-1]=P0,
    如果P1P2相对于P0P1在点P1向左旋转(共线也不行)，则P0,P1一定是凸包上的点；
    否则P1一定不是凸包上的点，应该将其出栈。
    比如下面左图中点1就一定是凸包上的点，右图中点1就一定不是凸包上的点，因为
    可以连接点0,2构成凸包的边
          2
          |
          |                _____2
          1              1
          /              /
    ____0/          ____0/

====================================================================================*/

void Graham_scan(TPoint PointSet[],TPoint ch[],int n,int &len)
{
    int i,j,k=0,top=2;
    TPoint tmp;
    
    //选取PointSet中y坐标最小的点PointSet[k]，如果这样的点右多个，则取最左边的一个
    for(i=1;i<n;i++)
        if ((PointSet[i].y<PointSet[k].y)||((PointSet[i].y==PointSet[k].y)&&(PointSet[i].x<PointSet[k].x)))
        k=i;
    tmp=PointSet[0];
    PointSet[0]=PointSet[k];
    PointSet[k]=tmp; //现在PointSet中y坐标最小的点在PointSet[0]
    for (i=1;i<n-1;i++) //对顶点按照相对PointSet[0]的极角从小到大进行排序，极角相同的按照距离PointSet[0]从近到远进行排序
        {  k=i;
            for (j=i+1;j<n;j++)
                if ( (multiply(PointSet[j],PointSet[k],PointSet[0])>0)
                    ||
                    ( (multiply(PointSet[j],PointSet[k],PointSet[0])==0)
                      &&(distance(PointSet[0],PointSet[j])<distance(PointSet[0],PointSet[k])) )
                  ) k=j;
            tmp=PointSet[i];
            PointSet[i]=PointSet[k];
            PointSet[k]=tmp;
        }
    ch[0]=PointSet[0];
    ch[1]=PointSet[1];
    ch[2]=PointSet[2];  
    for (i=3;i<n;i++)
        {
            while (multiply(PointSet[i],ch[top],ch[top-1])>=0) top--;
            ch[++top]=PointSet[i];
            }
    len=top+1;
}