2018-Reconquista-ComputationGeometry

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

#include<bits/stdc++.h>
#define ERROR 1e99
#define N 100005
using namespace std;
typedef double data;
typedef long long LL;
const data PI=acos(-1.);
const data EPS=1E-10;
int Zero(data O){return fabs(O)<=EPS;}
void Read(data &now){double IO;scanf("%lf",&IO);now=IO;}
struct Point{
    data x,y;
    void read(){Read(x);Read(y);}
    void Print(){printf("%lf %lf\n",(double)x,(double)y);}
    void init(data _x,data _y){x=_x;y=_y;}
    data operator * (const Point &b)const{
        return x*b.y-y*b.x;
    }
    Point operator * (const data &b)const{
        return (Point){x*b,y*b};
    }
    data operator ^ (const Point &b)const{
        return x*b.x+y*b.y;
    }
    Point operator - (const Point &b)const{
        return (Point){x-b.x,y-b.y};
    }
    Point operator + (const Point &b)const{
        return (Point){x+b.x,y+b.y};
    }
    int operator < (const Point &b)const{
        return x<b.x||x==b.x&&y<b.y;
    }
    int operator == (const Point &b)const{
        return Zero(x-b.x)&&Zero(y-b.y);
    }
    void Rotate(data alpha){
        data _x=x*cos(alpha)-y*sin(alpha);
        data _y=y*cos(alpha)+x*sin(alpha);
        x=_x;y=_y;
    }
}ERRORPOINT=(Point){ERROR,ERROR};
struct Line{
    Point p,q;
};
struct Circle{
    Point O;data r;
    void read(){O.read();Read(r);}
};
data sqr(data x){return x*x;}
data dist(Point A,Point B){return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y));}
int Intersection(Line A,Line B,Point &ret){
    ret=A.p;
    data u=(A.p-B.p)*(B.p-B.q);
    data v=(A.p.x-A.q.x)*(B.p.y-B.q.y)-(A.p.y-A.q.y)*(B.p.x-B.q.x);
    if (Zero(v)) return 0;u/=v;
    ret.x+=(A.q.x-A.p.x)*u;
    ret.y+=(A.q.y-A.p.y)*u;
    return 1;
}
int Intersection(Circle C,Line L,Point &A,Point &B){
    Point p=C.O;p.x+=L.p.y-L.q.y;p.y+=L.q.x-L.p.x;
    Intersection((Line){p,C.O},L,p);
    data tmp=sqr(C.r)-sqr(p.x-C.O.x)-sqr(p.y-C.O.y);
    if (tmp<-EPS) return 0;
    data t=sqrt(max(tmp,0.))/dist(L.p,L.q);
    A=p+(L.q-L.p)*t;B=p-(L.q-L.p)*t;
    return 1;
}
int Intersection(Circle C1,Circle C2,Point &A, Point &B) {
    double d = dist(C1.O,C2.O);
    if (d<fabs(C1.r-C2.r)-EPS || d>fabs(C1.r+C2.r)+EPS) return 0;
    double cosa = (sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2 * C1.r * d);
    double sina = sqrt(max(0., 1. - sqr(cosa)));
    A = B = C1.O;
    A.x += C1.r / d * ((C2.O.x - C1.O.x) * cosa + (C2.O.y - C1.O.y) * -sina);
    A.y += C1.r / d * ((C2.O.x - C1.O.x) * sina + (C2.O.y - C1.O.y) * cosa);
    B.x += C1.r / d * ((C2.O.x - C1.O.x) * cosa + (C2.O.y - C1.O.y) * sina);
    B.y += C1.r / d * ((C2.O.x - C1.O.x) *-sina + (C2.O.y - C1.O.y) * cosa);
    return 1;
}

//此处默认data是int 
int Polar_cmp(const Point &A,const Point &B){
    if (A.y==0&&B.y==0) return A.x>0&&B.x<0;
    if (A.y==0) return A.x>0||B.y<0;
    if (B.y==0) return B.x<0&&A.y>0;
    if ((A.y>0)^(B.y>0)) return A.y>0;
    return (A*B)>0;
}

vector<Point>Make_Convex(vector<Point>now){
    sort(now.begin(),now.end());
    if (now.size()<3) return now;
    for (int i=0;i<now.size();i++)
        now[i].Print();
    puts("");

    vector<Point>up;up.clear();
    up.push_back(now[0]);
    for (int i=1;i<now.size();i++){
        for (;up.size()>1&&(up[up.size()-2]-up.back())*(now[i]-up.back())<=0;up.pop_back());
        up.push_back(now[i]);
    }
    vector<Point>down;down.clear();
    down.push_back(now.back());
    for (int i=now.size()-1;i>=0;i--){
        for (;down.size()>1&&(down[down.size()-2]-down.back())*(now[i]-down.back())<=0;down.pop_back());
        down.push_back(now[i]);
    }

    vector<Point>ret;ret.clear();
    for (int i=down.size()-1;i>0;i--)
        ret.push_back(down[i]);
    for (int i=up.size()-1;i>0;i--)
        ret.push_back(up[i]);
    return ret;
}

//以下是最小圆覆盖

//两个点 
Circle make_2(Point p,Point q){
    Circle ret;
    ret.O.x=(p.x+q.x)/2.0;
    ret.O.y=(p.y+q.y)/2.0;
    ret.r=dist(ret.O,p);
    return ret;
}
//返回三角形的外心  
Circle make_3(const Point &a,const Point &b,const Point &c)  {
    data a1=b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;  
    data a2=c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;  
    data d=a1*b2-a2*b1;
    Circle ret;
    ret.O.x=a.x+(c1*b2-c2*b1)/d;
    ret.O.y=a.y+(a1*c2-a2*c1)/d;
    ret.r=dist(ret.O,a);
    return ret;
} 
bool IN(Circle now,Point a){
    return (dist(now.O,a)<=now.r+EPS);
}
Circle Smallest_Circle(Point *a,int n){
    if (n==1) return (Circle){a[1],0};
    srand(2333);
    random_shuffle(a+1,a+n+1);
    Circle cur=make_2(a[1],a[2]);
    for (int i=2;i<=n;i++)
        if (!IN(cur,a[i])){
            cur=make_2(a[1],a[i]);
            for (int j=1;j<i;j++)
                if (!IN(cur,a[j])){
                    cur=make_2(a[i],a[j]);
                    for (int k=1;k<j;k++)
                        if (!IN(cur,a[k]))
                            cur=make_3(a[i],a[j],a[k]);
                }
        }
    return cur;
}

#include

#define ERROR 1e99

#define N 100005

using namespace std;

typedef double data;

typedef long long LL;

const data PI=acos(-1.);

const data EPS=1E-10;

int Zero(data O){return fabs(O)<=EPS;}

void Read(data &now){double IO;scanf("%lf",&IO);now=IO;}

struct Point{

data x,y;

void read(){Read(x);Read(y);}

void Print(){printf("%lf %lf\n",(double)x,(double)y);}

void init(data _x,data _y){x=_x;y=_y;}

data operator * (const Point &b)const{

return x*b.y-y*b.x;

}

Point operator * (const data &b)const{

return (Point){x*b,y*b};

}

data operator ^ (const Point &b)const{

return x*b.x+y*b.y;

}

Point operator - (const Point &b)const{

return (Point){x-b.x,y-b.y};

}

Point operator + (const Point &b)const{

return (Point){x+b.x,y+b.y};

}

int operator < (const Point &b)const{

return x

}

int operator == (const Point &b)const{

return Zero(x-b.x)&&Zero(y-b.y);

}

void Rotate(data alpha){

data _x=x*cos(alpha)-y*sin(alpha);

data _y=y*cos(alpha)+x*sin(alpha);

x=_x;y=_y;

}

}ERRORPOINT=(Point){ERROR,ERROR};

struct Line{

Point p,q;

};

struct Circle{

Point O;data r;

void read(){O.read();Read(r);}

};

data sqr(data x){return x*x;}

data dist(Point A,Point B){return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y));}

int Intersection(Line A,Line B,Point &ret){

ret=A.p;

data u=(A.p-B.p)*(B.p-B.q);

data v=(A.p.x-A.q.x)*(B.p.y-B.q.y)-(A.p.y-A.q.y)*(B.p.x-B.q.x);

if (Zero(v)) return 0;u/=v;

ret.x+=(A.q.x-A.p.x)*u;

ret.y+=(A.q.y-A.p.y)*u;

return 1;

}

int Intersection(Circle C,Line L,Point &A,Point &B){

Point p=C.O;p.x+=L.p.y-L.q.y;p.y+=L.q.x-L.p.x;

Intersection((Line){p,C.O},L,p);

data tmp=sqr(C.r)-sqr(p.x-C.O.x)-sqr(p.y-C.O.y);

if (tmp<-EPS) return 0;

data t=sqrt(max(tmp,0.))/dist(L.p,L.q);

A=p+(L.q-L.p)*t;B=p-(L.q-L.p)*t;

return 1;

}

int Intersection(Circle C1,Circle C2,Point &A, Point &B) {

double d = dist(C1.O,C2.O);

if (dfabs(C1.r+C2.r)+EPS) return 0;

double cosa = (sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2 * C1.r * d);

double sina = sqrt(max(0., 1. - sqr(cosa)));

A = B = C1.O;

A.x += C1.r / d * ((C2.O.x - C1.O.x) * cosa + (C2.O.y - C1.O.y) * -sina);

A.y += C1.r / d * ((C2.O.x - C1.O.x) * sina + (C2.O.y - C1.O.y) * cosa);

B.x += C1.r / d * ((C2.O.x - C1.O.x) * cosa + (C2.O.y - C1.O.y) * sina);

B.y += C1.r / d * ((C2.O.x - C1.O.x) *-sina + (C2.O.y - C1.O.y) * cosa);

return 1;

}

//此处默认data是int

int Polar_cmp(const Point &A,const Point &B){

if (A.y==0&&B.y==0) return A.x>0&&B.x<0;

if (A.y==0) return A.x>0||B.y<0;

if (B.y==0) return B.x<0&&A.y>0;

if ((A.y>0)^(B.y>0)) return A.y>0;

return (A*B)>0;

}

vectorMake_Convex(vectornow){

sort(now.begin(),now.end());

if (now.size()<3) return now;

for (int i=0;i

now[i].Print();

puts("");

vectorup;up.clear();

up.push_back(now[0]);

for (int i=1;i

for (;up.size()>1&&(up[up.size()-2]-up.back())*(now[i]-up.back())<=0;up.pop_back());

up.push_back(now[i]);

}

vectordown;down.clear();

down.push_back(now.back());

for (int i=now.size()-1;i>=0;i--){

for (;down.size()>1&&(down[down.size()-2]-down.back())*(now[i]-down.back())<=0;down.pop_back());

down.push_back(now[i]);

}

vectorret;ret.clear();

for (int i=down.size()-1;i>0;i--)

ret.push_back(down[i]);

for (int i=up.size()-1;i>0;i--)

ret.push_back(up[i]);

return ret;

}

//以下是最小圆覆盖

//两个点

Circle make_2(Point p,Point q){

Circle ret;

ret.O.x=(p.x+q.x)/2.0;

ret.O.y=(p.y+q.y)/2.0;

ret.r=dist(ret.O,p);

return ret;

}

//返回三角形的外心

Circle make_3(const Point &a,const Point &b,const Point &c) {

data a1=b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;

data a2=c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;

data d=a1*b2-a2*b1;

Circle ret;

ret.O.x=a.x+(c1*b2-c2*b1)/d;

ret.O.y=a.y+(a1*c2-a2*c1)/d;

ret.r=dist(ret.O,a);

return ret;

}

bool IN(Circle now,Point a){

return (dist(now.O,a)<=now.r+EPS);

}

Circle Smallest_Circle(Point *a,int n){

if (n==1) return (Circle){a[1],0};

srand(2333);

random_shuffle(a+1,a+n+1);

Circle cur=make_2(a[1],a[2]);

for (int i=2;i<=n;i++)

if (!IN(cur,a[i])){

cur=make_2(a[1],a[i]);

for (int j=1;j

if (!IN(cur,a[j])){

cur=make_2(a[i],a[j]);

for (int k=1;k

if (!IN(cur,a[k]))

cur=make_3(a[i],a[j],a[k]);

}

}

return cur;

}