计算几何模板

/*///////////////////////////////////
***************Content***************
范数（模的平方）
向量的模(求线段距离/两点距离)
点乘
叉乘
判断向量垂直
判断a1-a2与b1-b2垂直
判断线段垂直
判断向量平行
判断a1-a2与b1-b2平行
判断线段平行
求点在线段上的垂足
求点关于直线的对称点
判断P0,P1,P2三点位置关系,Vector p0-p2 在p0-p1的相对位置
判断p1-p2与p3-p4是否相交
判断线段是否相交
a,b两点间的距离
点到直线的距离
点到线段的距离
两线段间的距离(相交为0)
求两线段的交点
求直线交点
中垂线
向量a,b的夹角，范围[0,180]
向量（点）极角，范围[-180,180]
角度排序(从x正半轴起逆时针一圈)范围为[0,180)
判断点在多边形内部
凸包(CCW/CW)
求向量A,向量B构成三角形的面积
(旋转卡壳)求平面中任意两点的最大距离
三点所成的外接圆
三点所成的内切圆
过点p求过该点的两条切线与圆的两个切点
极角(直线与x轴的角度)
点与圆的位置关系
已知点与切线求圆心
已知两直线和半径,求夹在两直线间的圆
求与给定两圆相切的圆
多边形(存储方式:点->线)
半平面交
求最近点对的距离(分治)
*////////////////////////////////////
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<vector>
using namespace std;
#define EPS (1e-10)
#define equals(a, b) (fabs(a-b)<EPS)
const double PI = acos(-1);
const double INF = 1e20;

static const int CCW=1;//逆时针
static const int CW=-1;//顺时针
static const int BACK=-2;//后面
static const int FRONT=2;//前面
static const int ON=0;//线段上

struct Point
{
    double x, y;
    Point(double x=0, double y=0):x(x), y(y) {}
};//点
typedef Point Vector;//向量

struct Segment
{
    Point p1, p2;
    Segment(Point p1=Point(), Point p2=Point()):p1(p1), p2(p2) {}
    double angle;
};//线段
typedef Segment Line;//直线
typedef vector<Point> Polygon;//多边形(存点)
typedef vector<Segment> Pol;//多边形（存边）

class Circle
{
public:
    Point c;
    double r;
    Circle(Point c=Point(), double r=0.0):c(c), r(r) {}
    Point point(double a)
    {
        return Point(c.x + cos(a)*r, c.y + sin(a)*r);
    }
};//圆

Point operator + (Point a, Point b)
{
    return Point(a.x+b.x, a.y+b.y);
}
Point operator - (Point a, Point b)
{
    return Point(a.x-b.x, a.y-b.y);
}
Point operator * (Point a, double p)
{
    return  Point(a.x*p, a.y*p);
}
Point operator / (Point a, double p)
{
    return Point(a.x/p, a.y/p);
}
//排序左下到右上
bool operator < (const Point &a,const Point &b)
{
    return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (Point a, Point b)
{
    return fabs(a.x-b.x)<EPS && fabs(a.y-b.y)<EPS;
}
//范数（模的平方）
double norm(Vector a)
{
    return a.x*a.x+a.y*a.y;
}
//向量的模(求线段距离/两点距离)
double abs(Vector a)
{
    return sqrt(norm(a));
}
//点乘
double dot(Vector a, Vector b)
{
    return a.x*b.x+a.y*b.y;
}
//叉乘
double cross(Vector a, Vector b)
{
    return a.x*b.y-a.y*b.x;
}
//判断向量垂直
bool isOrthgonal(Vector a, Vector b)
{
    return equals(dot(a, b), 0.0);
}
//判断a1-a2与b1-b2垂直
bool isOrthgonal(Point a1, Point a2, Point b1, Point b2)
{
    return isOrthgonal(a1-a2, b1-b2);
}
//判断线段垂直
bool isOrthgonal(Segment s1, Segment s2)
{
    return equals(dot(s1.p2-s1.p1, s2.p2-s2.p1), 0.0);
}
//判断向量平行
bool isParallel(Vector a, Vector b)
{
    return equals(cross(a, b), 0.0);
}
//判断a1-a2与b1-b2平行
bool isParallel(Point a1, Point a2, Point b1, Point b2)
{
    return isParallel(a1-a2, b1-b2);
}
//判断线段平行
bool isParallel(Segment s1, Segment s2)
{
    return equals(cross(s1.p2-s1.p1, s2.p2-s2.p1), 0.0);
}
//求点在线段上的垂足
Point project(Segment s, Point p)
{
    Vector base=s.p2-s.p1;
    double r=dot(p-s.p1, base)/norm(base);
    return s.p1+base*r;
}
//求点关于直线的对称点
Point reflect(Segment s, Point p)
{
    return p+(project(s, p)-p)*2.0;
}
//判断P0,P1,P2三点位置关系,Vector p0-p2 在p0-p1的相对位置
int ccw(Point p0, Point p1, Point p2)
{
    Vector a=p1-p0;
    Vector b=p2-p0;
    if( cross(a, b)>EPS ) return CCW;
    if( cross(a, b)<-EPS ) return CW;
    if( dot(a, b)<-EPS ) return BACK;
    if( norm(a)<norm(b) ) return FRONT;
    return ON;
}
//判断p1-p2与p3-p4是否相交
bool intersect(Point p1, Point p2, Point p3, Point p4)
{
    return (ccw(p1, p2, p3)*ccw(p1, p2, p4)<=0 &&ccw(p3, p4, p1)*ccw(p3, p4, p2)<=0);
}
//判断线段是否相交
bool intersect(Segment s1, Segment s2)
{
    return intersect(s1.p1, s1.p2, s2.p1, s2.p2);
}
//a,b两点间的距离
double getDistance(Point a, Point b)
{
    return abs(a-b);
}
//点到直线的距离
double getDistanceLP(Line l, Point p)
{
    return abs(cross(l.p2-l.p1, p-l.p1)/abs(l.p2-l.p1));
}
//点到线段的距离
double getDistanceSP(Segment s, Point p)
{
    if(dot(s.p2-s.p1, p-s.p1)<0.0) return abs(p-s.p1);
    if(dot(s.p1-s.p2, p-s.p2)<0.0) return abs(p-s.p2);
    return getDistanceLP(s, p);
}
//两线段间的距离(相交为0)
double getDistanceSS(Segment s1, Segment s2)
{
    if(intersect(s1, s2)) return 0.0;
    return min( min(getDistanceSP(s1,s2.p1), getDistanceSP(s1,s2.p2)),
                min(getDistanceSP(s2,s1.p1), getDistanceSP(s2,s1.p2)) );
}
//求两线段的交点
Point getCrossPoint(Segment s1, Segment s2)
{
    Vector base=s2.p2-s2.p1;
    double d1=abs(cross(base, s1.p1-s2.p1));
    double d2=abs(cross(base, s1.p2-s2.p1));
    double t=d1/(d1+d2);
    return s1.p1+(s1.p2-s1.p1)*t;
}
//求直线交点
Point intersectL(Segment a, Segment b)
{
    double x1=a.p1.x,y1=a.p1.y,x2=a.p2.x,y2=a.p2.y;
    double x3=b.p1.x,y3=b.p1.y,x4=b.p2.x,y4=b.p2.y;
    double k1=(x4-x3)*(y2-y1),k2=(x2-x1)*(y4-y3);
    double ans_x=(k1*x1-k2*x3+(y3-y1)*(x2-x1)*(x4-x3))/(k1-k2);
    double ans_y=(k2*y1-k1*y3+(x3-x1)*(y2-y1)*(y4-y3))/(k2-k1);
    return Point(ans_x,ans_y);
}
//中垂线
Line mid_vert(Line l)
{
    double x1=l.p1.x,y1=l.p1.y;
    double x2=l.p2.x,y2=l.p2.y;
    double xm=(x1+x2)/2,ym=(y1+y2)/2;
    Line s;
    s.p1.x=xm+ym-y1;
    s.p1.y=ym-xm+x1;
    s.p2.x=xm-ym+y1;
    s.p2.y=ym+xm-x1;
    return s;
}
//向量a,b的夹角，范围[0,180]
double Angle(Vector a, Vector b)
{
    return acos(dot(a, b)/(abs(a)*abs(b)));
}
//向量（点）极角，范围[-180,180]
double angle(Vector v)
{
    return atan2(v.y, v.x);
}
//角度排序(从x正半轴起逆时针一圈)范围为[0,180)
double SortAngle(Vector a)
{
    Point p0(0.0, 0.0);
    Point p1(a.x, a.y);
    Point p2(-1.0, 0.0);
    Vector b=p2;
    if(ccw(p0, p1, p2)==CW) return acos(dot(a, b)/(abs(a)*abs(b)));
    if(ccw(p0, p1, p2)==CCW) return 2*PI-acos(dot(a, b)/(abs(a)*abs(b)));
    if(ccw(p0, p1, p2)==BACK) return PI;
    else return 0;
}
/*
判断点在多边形内部
IN 2
ON 1
OUT 0
*/
int contains(Polygon g, Point p)
{
    int n=g.size();
    bool x=false;
    for(int i=0; i<n; i++)
    {
        Point a=g[i]-p;
        Point b=g[(i+1)%n]-p;
        if(abs(cross(a, b))<EPS && dot(a, b)<EPS) return 1;
        if(a.y>b.y) swap(a,b);
        if(a.y<EPS && EPS<b.y && cross(a,b)>EPS) x=!x;
    }
    return (x? 2 : 0);
}
//凸包(CCW/CW)
Polygon andrewScan(Polygon s)
{
    Polygon u, l;
    if(s.size()<3) return s;
    sort(s.begin(), s.end());
    u.push_back(s[0]);
    u.push_back(s[1]);
    l.push_back(s[s.size()-1]);
    l.push_back(s[s.size()-2]);
    for(int i=2; i<s.size(); i++)
    {
        for(int n=u.size(); n>=2 && ccw(u[n-2], u[n-1], s[i])!=CW; n--)
            u.pop_back();
        u.push_back(s[i]);
    }
    for(int i=s.size()-3; i>=0; i--)
    {
        for(int n=l.size(); n>=2 && ccw(l[n-2], l[n-1], s[i])!=CW; n--)
            l.pop_back();
        l.push_back(s[i]);
    }
    reverse(l.begin(), l.end());
    for(int i=u.size()-2; i>=1; i--) l.push_back(u[i]);
    return l;
}
//求向量A,向量B构成三角形的面积
double TriArea(Vector a, Vector b)
{
    return 0.5*abs(cross(a,b));
}
//求平面中任意两点的最大距离(旋转卡壳)
double RotatingCalipers(const Polygon& s)
{
    Polygon l;
    double dis, maxn=0.0;
    int len, i, k;
    l=andrewScan(s);
    len=l.size();
    if(len>=3)
    {
        for(i=0, k=2; i<len; i++)
        {
            while(cross(l[(k+1)%len]-l[i], l[(k+1)%len]-l[(i+1)%len])>=cross(l[k%len]-l[i], l[k%len]-l[(i+1)%len]))
                k++;
            dis=max(norm(l[k%len]-l[i]), norm(l[k%len]-l[(i+1)%len]));
            if(dis>maxn) maxn=dis;
        }
    }
    else maxn=norm(l[1]-l[0]);
    return maxn;
}
//三点所成的外接圆
Circle CircumscribedCircle(Point a, Point b, Point c)
{
    double x=0.5*(norm(b)*c.y+norm(c)*a.y+norm(a)*b.y-norm(b)*a.y-norm(c)*b.y-norm(a)*c.y)
             /(b.x*c.y+c.x*a.y+a.x*b.y-b.x*a.y-c.x*b.y-a.x*c.y);
    double y=0.5*(norm(b)*a.x+norm(c)*b.x+norm(a)*c.x-norm(b)*c.x-norm(c)*a.x-norm(a)*b.x)
             /(b.x*c.y+c.x*a.y+a.x*b.y-b.x*a.y-c.x*b.y-a.x*c.y);
    Point O(x, y);
    double r=abs(O-a);
    Circle m(O, r);
    return m;
}
//三点所成的内切圆
Circle InscribedCircle(Point a, Point b, Point c)
{
    double A=abs(b-c), B=abs(a-c), C=abs(a-b);
    double x=(A*a.x+B*b.x+C*c.x)/(A+B+C);
    double y=(A*a.y+B*b.y+C*c.y)/(A+B+C);
    Point O(x, y);
    Line l(a, b);
    double r=getDistanceLP(l, O);
    Circle m(O, r);
    return m;
}
//过点p求过该点的两条切线与圆的两个切点
Segment TangentLineThroughPoint(Circle m, Point p)
{
    Point c=m.c;
    double l=abs(c-p);
    double r=m.r;
    double k=(2*r*r-l*l+norm(p)-norm(c)-2*p.y*c.y+2*c.y*c.y)/(2*(p.y-c.y));
    double A=1+(p.x-c.x)*(p.x-c.x)/((p.y-c.y)*(p.y-c.y));
    double B=-(2*k*(p.x-c.x)/(p.y-c.y)+2*c.x);
    double C=c.x*c.x+k*k-r*r;
    double x1, x2, y1, y2;

    x1=(-B-sqrt(B*B-4*A*C))/(2*A);
    x2=(-B+sqrt(B*B-4*A*C))/(2*A);
    y1=(2*r*r-l*l+norm(p)-norm(c)-2*(p.x-c.x)*x1)/(2*(p.y-c.y));
    y2=(2*r*r-l*l+norm(p)-norm(c)-2*(p.x-c.x)*x2)/(2*(p.y-c.y));
    Point p1(x1, y1), p2(x2, y2);
    Segment L(p1, p2);
    return L;
}
//极角(直线与x轴的角度)
double PolarAngle(Vector a)
{
    Point p0(0.0, 0.0);
    Point p1(1.0, 0.0);
    Point p2(a.x, a.y);
    Vector b=p1;
    double ans=0;
    if(ccw(p0, p1, p2)==CW) ans=180-acos(dot(a, b)/(abs(a)*abs(b)))*180/acos(-1);
    else if(ccw(p0, p1, p2)==CCW) ans=acos(dot(a, b)/(abs(a)*abs(b)))*180/acos(-1);
    else ans=0;
    if(ans>=180) ans-=180;
    if(ans<0) ans+=180;
    return ans;
}
//点与圆的位置关系
int CircleContain(Circle m, Point p)
{
    double r=m.r;
    double l=abs(p-m.c);
    if(r>l) return 2;
    if(r==l) return 1;
    if(r<l) return 0;
}
//已知点与切线求圆心
void CircleThroughAPointAndTangentToALineWithRadius(Point p, Line l, double r)
{
    Point m=project(l, p);
    if(abs(p-m)>2*r)
    {
        printf("[]\n");
    }
    else if(abs(p-m)==2*r)
    {
        Circle c((p+m)/2, r);
        printf("[(%.6f,%.6f)]\n", c.c.x, c.c.y);
    }
    else if(abs(p-m)<EPS)
    {
        Point m0(m.x+10, m.y);
        if(abs(m0-project(l, m0))<EPS) m0.y+=20;
        Point m1=project(l, m0);
        Circle c1(m-(m0-m1)/abs(m0-m1)*r, r);
        Circle c2(m+(m0-m1)/abs(m0-m1)*r, r);
        if(c1.c.x>c2.c.x) swap(c1, c2);
        else if(c1.c.x==c2.c.x && c1.c.y>c2.c.y) swap(c1, c2);
        printf("[(%.6f,%.6f),(%.6f,%.6f)]\n", c1.c.x, c1.c.y, c2.c.x, c2.c.y);
    }
    else if(abs(p-m)<2*r)
    {
        double s=abs(p-m);
        double d=sqrt(r*r-(r-s)*(r-s));
        Point m1, m2;
        m1=(m+(l.p1-l.p2)/abs(l.p1-l.p2)*d);
        m2=(m-(l.p1-l.p2)/abs(l.p1-l.p2)*d);
        Circle c1(m1+(p-m)/abs(p-m)*r, r);
        Circle c2(m2+(p-m)/abs(p-m)*r, r);
        if(c1.c.x>c2.c.x) swap(c1, c2);
        else if(c1.c.x==c2.c.x && c1.c.y>c2.c.y) swap(c1, c2);
        printf("[(%.6f,%.6f),(%.6f,%.6f)]\n", c1.c.x, c1.c.y, c2.c.x, c2.c.y);
    }
    return ;
}

bool cmp_CircleTangentToTwoLinesWithRadius(Circle x, Circle y)
{
    if(x.c.x==y.c.x) return x.c.y<y.c.y;
    else return x.c.x<y.c.x;
}
//已知两直线和半径,求夹在两直线间的圆
void CircleTangentToTwoLinesWithRadius(Line l1, Line l2, double r)
{
    Point p=intersectL(l1, l2);
    Vector a, b;
    l1.p2=p+p-l1.p1;
    l2.p2=p+p-l2.p1;

    a=(l1.p1-p)/abs(l1.p1-p), b=(l2.p2-p)/abs(l2.p2-p);
    Circle c1(p+(a+b)/abs(a+b)*(r/sin(Angle(a,a+b))),r);

    a=(l1.p1-p)/abs(l1.p1-p), b=(l2.p1-p)/abs(l2.p1-p);
    Circle c2(p+(a+b)/abs(a+b)*r/sin(Angle(a,a+b)),r);

    a=(l1.p2-p)/abs(l1.p2-p), b=(l2.p1-p)/abs(l2.p1-p);
    Circle c3(p+(a+b)/abs(a+b)*(r/sin(Angle(a,a+b))),r);

    a=(l1.p2-p)/abs(l1.p2-p), b=(l2.p2-p)/abs(l2.p2-p);
    Circle c4(p+(a+b)/abs(a+b)*(r/sin(Angle(a,a+b))),r);

    vector<Circle> T;
    T.push_back(c1);
    T.push_back(c2);
    T.push_back(c3);
    T.push_back(c4);
    sort(T.begin(), T.end(), cmp_CircleTangentToTwoLinesWithRadius);
    printf("[(%.6f,%.6f),(%.6f,%.6f),(%.6f,%.6f),(%.6f,%.6f)]\n", T[0].c.x, T[0].c.y, T[1].c.x, T[1].c.y, T[2].c.x, T[2].c.y, T[3].c.x, T[3].c.y);
}
//求与给定两圆相切的圆
void CircleTangentToTwoDisjointCirclesWithRadius(Circle C1, Circle C2)
{
    double d = abs(C1.c-C2.c);
    if(d<EPS) printf("[]\n");
    else if(fabs(C1.r+C2.r)<d) printf("[]\n");
    else if(fabs(C1.r-C2.r)>d) printf("[]\n");
    else
    {
        double sita = angle(C2.c - C1.c);
        double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2 * C1.r*d));
        Point p1 = C1.point(sita - da), p2 = C1.point(sita + da);
        if(p1.x>p2.x) swap(p1, p2);
        else if(p1.x==p2.x && p1.y>p2.y) swap(p1, p2);
        if (p1 == p2) printf("[(%.6f,%.6f)]\n", p1.x, p1.y);
        else printf("[(%.6f,%.6f),(%.6f,%.6f)]\n", p1.x, p1.y, p2.x, p2.y);
    }
}
//多边形(存储方式:点->线)
Pol Polygon_to_Pol(Polygon L)
{
    Pol S;
    Segment l;
    for(int i=0; i+1<L.size(); i++)
    {
        l.p1=L[i];
        l.p2=L[i+1];
        swap(l.p1, l.p2);//注意顺序
        l.angle=angle(l.p2-l.p1);
        S.push_back(l);
    }
    l.p1=L[L.size()-1];
    l.p2=L[0];
    swap(l.p1, l.p2);//注意顺序
    l.angle=angle(l.p2-l.p1);
    S.push_back(l);
    return S;
}
//半平面交
bool SortAnglecmp(Segment a, Segment b)
{
    if(fabs(a.angle-b.angle)>EPS)
        return a.angle>b.angle;
    return ccw(b.p1, b.p2, a.p1)!=CW;
}

int intersection_of_half_planes(Pol s)
{
    Segment deq[1505];
    Segment l[1505];
    Point p[1505];
    memset(deq, 0, sizeof(deq));
    memset(l, 0, sizeof(l));
    memset(p, 0, sizeof(p));

    sort(s.begin(), s.end(), SortAnglecmp);
    int cnt=0;
    for(int i=0; i<s.size(); i++)
        if(fabs(s[i].angle-l[cnt].angle)>EPS)
            l[++cnt]=s[i];
    int le=1,ri=1;
    for(int i=1; i<=cnt; i++)
    {
        while(ri>le+1 && ccw(l[i].p1, l[i].p2, intersectL(deq[ri-1],deq[ri-2]))==CW) ri--;
        while(ri>le+1 && ccw(l[i].p1, l[i].p2, intersectL(deq[le],deq[le+1]))==CW) le++;
        deq[ri++]=l[i];
    }
    while(ri>le+2 && ccw(deq[le].p1, deq[le].p2, intersectL(deq[ri-1],deq[ri-2]))==CW) ri--;
    while(ri>le+2 && ccw(deq[ri-1].p1, deq[ri-1].p2, intersectL(deq[le],deq[le+1]))==CW) le++;
    //***************getArea******************
    /*
    if(ri<=le+2)
    {
        printf("0.00\n");
        return 0;
    }
    deq[ri]=deq[le];
    cnt=0;
    for(int i=le; i<ri; i++)
        p[++cnt]=intersectL(deq[i],deq[i+1]);
    double ans=0.0;
    for(int i=2; i<cnt; i++)
        ans+=fabs(TriArea(p[i]-p[1],p[i+1]-p[1]));
    printf("%.2f\n", ans);
    return 0;
    */
    //****************************************
    if(ri>le+2) return 1;
    return 0;
}
//求最近点对的距离(分治)
//************************************************************************
bool cmpxy(Point a, Point b)
{
    if(a.x != b.x) return a.x < b.x;
    else return a.y < b.y;
}
bool cmpy(Point a, Point b)
{
    return a.y < b.y;
}
int n;
Point closest_p[100010];//将点都存入closest_p中！
Point closest_tmpt[100010];
double Closest_Pair(int left,int right)
{
    double d = INF;
    if(left == right) return d;
    if(left + 1 == right)
        return getDistance(closest_p[left],closest_p[right]);
    int mid = (left+right)/2;
    double d1 = Closest_Pair(left,mid);
    double d2 = Closest_Pair(mid+1,right);
    d = min(d1,d2);
    int k = 0;
    for(int i = left; i <= right; i++)
        if(fabs(closest_p[mid].x - closest_p[i].x) <= d)
            closest_tmpt[k++] = closest_p[i];
    sort(closest_tmpt,closest_tmpt+k,cmpy);
    for(int i = 0; i <k; i++)
        for(int j = i+1; j < k && closest_tmpt[j].y - closest_tmpt[i].y < d; j++)
            d = min(d,getDistance(closest_tmpt[i],closest_tmpt[j]));
    return d;
}
double Closest()//直接调用此函数即可
{
    sort(closest_p,closest_p+n,cmpxy);
    return Closest_Pair(0,n-1)/2;
}
//************************************************************************
int main()
{
    while(scanf("%d",&n)==1 && n)
    {
        for(int i = 0; i < n; i++)
            scanf("%lf%lf",&closest_p[i].x,&closest_p[i].y);
        printf("%.2lf\n",Closest());
    }
    return 0;
}