嘘~ 正在从服务器偷取页面 . . .

OI模板-计算几何


OI模板-计算几何

二维凸包

Andrew

时间复杂度 \(O(n \log n)\)

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
#define N (int)(1e5+15)
int n;
struct vct
{
    double x,y;
    vct operator-(const vct &o)const
    {
        return (vct){x-o.x,y-o.y};
    }
}p[N],ans[N];
int stk[N],top,used[N];
int cmp(vct a,vct b)
{
    return a.x==b.x?a.y<b.y:a.x<b.x;
}
double cross(vct a,vct b)
{
    return a.x*b.y-a.y*b.x;
}
double dis(vct a,vct b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cout << fixed << setprecision(2);
    cin >> n;
    for(int i=1; i<=n; i++)
        cin >> p[i].x >> p[i].y;
    sort(p+1,p+1+n,cmp);
    stk[++top]=1;
    for(int i=2; i<=n; i++)
    {
        while(top>1&&cross(p[stk[top]]-p[stk[top-1]],p[i]-p[stk[top]])<=0)
            used[stk[top--]]=0;
        used[i]=1;
        stk[++top]=i;
    }
    int tmp=top;
    for(int i=n-1; i>=1; i--)
    {
        if(!used[i])
        {
            while(top>tmp&&cross(p[stk[top]]-p[stk[top-1]],p[i]-p[stk[top]])<=0)
                used[stk[top--]]=0;
            used[i]=1;
            stk[++top]=i;
        }
    }
    for(int i=1; i<=top; i++)
        ans[i]=p[stk[i]];
    double res=0;
    for(int i=2; i<=top; i++)
        res+=dis(ans[i-1],ans[i]);
    cout << res << endl;
    return 0;
}

Graham

时间复杂度 \(O(n\log n)\)

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define N (int)(1e5+15)
int n;
struct node
{
    double x,y;
}p[N],stk[N];
double cross(node a1,node a2,node b1,node b2)
{
    return (a2.x-a1.x)*(b2.y-b1.y)-(b2.x-b1.x)*(a2.y-a1.y);
}
double dis(node a,node b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
bool cmp(node a,node b)
{
    double tmp=cross(p[1],a,p[1],b);
    if(tmp>0)return 1;
    if(tmp==0&&dis(p[1],a)<=dis(p[1],b))
        return 1;
    return 0;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cout << fixed << setprecision(2);
    cin >> n;
    for(int i=1; i<=n; i++)
    {
        cin >> p[i].x >> p[i].y;
        if(i!=1&&p[i].y<=p[1].y&&p[i].x<=p[1].x)
            swap(p[1],p[i]);
    }
    sort(p+2,p+1+n,cmp);
    stk[1]=p[1];
    int top=1;
    for(int i=2; i<=n; i++)
    {
        while(top>1&&cross(stk[top-1],stk[top],stk[top],p[i])<=0)
            --top;
        stk[++top]=p[i];
    }
    double ans=0;
    stk[top+1]=p[1];
    for(int i=1; i<=top; i++)
        ans+=dis(stk[i],stk[i+1]);
    cout << ans << endl;
    return 0;
}

平面最近点对

P7883 平面最近点对(加强加强版)

时间复杂度 \(O(n\log n)\)

注意下面的代码输出的是平方

改为求sqrt一定要改成fabs !!!!!!!!!!!

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF ((int)0x3f3f3f3f3f3f3f3f)
#define inf ((int)0xc0c0c0c0c0c0c0c0)
#define N (int)(4e5+15)
int n;
struct node
{
    int x,y;
}p[N],tmp[N];
#define pf(x) ((x)*(x))
int dis(node a,node b)
{
    return pf(a.x-b.x)+pf(a.y-b.y);
}
int cmpx(node a,node b){return a.x<b.x;}
int cmpy(node a,node b){return a.y<b.y;}
int cdq(int l,int r)
{
    if(l==r)return INF;
    if(l+1==r)
    {
        if(p[l].y>p[r].y)swap(p[l],p[r]);
        return dis(p[l],p[r]);
    }
    int mid=(l+r)>>1,t=p[mid].x;
    int d=min(cdq(l,mid),cdq(mid+1,r));
    inplace_merge(p+l,p+mid+1,p+r+1,cmpy);
    int cnt=0;
    for(int i=l; i<=r; i++)
        if(pf(p[i].x-t)<d)
            tmp[++cnt]=p[i];
    for(int i=1; i<=cnt; i++)
        for(int j=i+1; j<=cnt&&pf(tmp[i].y-tmp[j].y)<d; j++)
            d=min(d,dis(tmp[i],tmp[j]));
    return d;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cout << fixed << setprecision(4);
    cin >> n;
    for(int i=1; i<=n; i++)
        cin >> p[i].x >> p[i].y;
    sort(p+1,p+1+n,cmpx);
    cout << cdq(1,n) << endl;
    return 0;
}

半平面交

P4196 [CQOI2006]凸多边形 /【模板】半平面交

  1. 给出了凸多边形,逆时针给出各个顶点的坐标

时间复杂度 \(O(n\log n)\) (不是题目里的 \(n\)

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF ((int)0x3f3f3f3f3f3f3f3f)
#define inf ((int)0xc0c0c0c0c0c0c0c0)
#define N (int)(1e5+15)
const double eps=1e-10;
struct vct{double x,y;}p[N],in[N];
#define pf(x) ((x)*(x))
vct operator+(vct a,vct b){return {a.x+b.x,a.y+b.y};}
vct operator-(vct a,vct b){return {a.x-b.x,a.y-b.y};}
vct operator*(vct a,double b){return {a.x*b,a.y*b};}
vct operator/(vct a,double b){return {a.x/b,a.y/b};}
double cross(vct a,vct b){return a.x*b.y-a.y*b.x;}
double dot(vct a,vct b){return a.x*b.x+a.y*b.y;}
double len(vct a){return sqrt(pf(a.x)+pf(a.y));}
struct line
{
    vct p,way;
    double k;
    void mkline(vct a,vct b)
    {
        p=a;way=b;
        k=atan2(b.y,b.x);
    }
}a[N];
double dcmp(double x)
{
    if(fabs(x)<=eps)return 0;
    if(x>eps)return 1;
    return -1;
}
bool operator<(line a,line b)
{
    int d=dcmp(a.k-b.k);
    if(d==0)return dcmp(cross(b.p-a.p,b.way))>0;
    return d<0;
}
bool onright(vct a,line b){return dcmp(cross(b.way,a-b.p))<0;}
vct intersect(line a,line b)
{
    double t=cross(b.way,a.p-b.p)/cross(a.way,b.way);
    return a.p+a.way*x;
}
int st,en;
line que[N];
bool halfplane(int n)
{
    que[st=en=1]=a[1];
    for(int i=2; i<=n; i++)
    {
        while(st<en&&onright(p[en],a[i]))--en;
        while(st<en&&onright(p[st+1],a[i]))++st;
        que[++en]=a[i];
        if(st<en)p[en]=intersect(que[en-1],que[en]);
    }
    while(st<en&&onright(p[en],que[st]))--en;
    if(en-st<=1)return 0;
    p[st]=intersect(que[st],que[en]);
    return 1;
}
double calc(vct *a,int n)
{
    double res=0;
    for(int i=1; i<n; i++)
        res+=cross(a[i],a[i+1]);
    res+=cross(a[n],a[1]);
    if(fabs(res/=2)<=eps)res=0;
    return res;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cout << fixed << setprecision(3);
    int n,o=0,oo=0;
    cin >> n;
    while(n--)
    {
        int m;cin >> m;
        for(int i=1; i<=m; i++)
            cin >> in[i].x >> in[i].y;
        for(int i=1; i<m; i++)
            a[++oo].mkline(in[i],in[i+1]-in[i]);
        a[++oo].mkline(in[m],in[1]-in[m]);
    }
    sort(a+1,a+1+oo);
    for(int i=1; i<=oo; i++)
        if(i==1||dcmp(a[i].k-a[i-1].k))
            a[++o]=a[i];
    if(!halfplane(o))cout << 0.0 << endl;
    else cout << calc(p+st-1,en-st+1) << endl;
    return 0;
}

poj2451

  1. 给出了向量

时间复杂度同上

注意POJ上提交要用printf() + %f ,所以下面这个代码交上去的时候要改一下

#include <cstdio>
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <iomanip>
using namespace std;
#define int long long
// #define double long double
#define INF ((int)0x3f3f3f3f3f3f3f3f)
#define inf ((int)0xc0c0c0c0c0c0c0c0)
#define N (int)(1e5+15)
const double eps=1e-10;
struct vct{double x,y;}p[N];
#define pf(x) ((x)*(x))
vct operator+(vct a,vct b){return {a.x+b.x,a.y+b.y};}
vct operator-(vct a,vct b){return {a.x-b.x,a.y-b.y};}
vct operator*(vct a,double b){return {a.x*b,a.y*b};}
vct operator/(vct a,double b){return {a.x/b,a.y/b};}
double cross(vct a,vct b){return a.x*b.y-a.y*b.x;}
double dot(vct a,vct b){return a.x*b.x+a.y*b.y;}
double len(vct a){return sqrt(pf(a.x)+pf(a.y));}
struct line
{
    vct p,way;
    double k;
    void mkline(vct a,vct b)
    {
        p=a;way=b;
        k=atan2(b.y,b.x);
    }
    void mk(double x1,double y1,double x2,double y2)
    {
        mkline({x1,y1},{x2-x1,y2-y1});
    }
}a[N];
int dcmp(double x)
{
    if(fabs(x)<=eps)return 0;
    return x>eps?1:-1;
}
bool operator<(line a,line b)
{
    int d=dcmp(a.k-b.k);
    if(d==0)return dcmp(cross(b.p-a.p,b.way))>0;
    return d<0;
}
bool onright(vct a,line b){return dcmp(cross(b.way,a-b.p))<0;}
vct intersect(line a,line b)
{
    double x=cross(b.way,a.p-b.p)/cross(a.way,b.way);
    return a.p+a.way*x;
}
int st,en;
line que[N];
bool halfplane(int n)
{
    que[st=en=1]=a[1];
    for(int i=2; i<=n; i++)
    {
        while(st<en&&onright(p[en],a[i]))--en;
        while(st<en&&onright(p[st+1],a[i]))++st;
        que[++en]=a[i];
        if(st<en)p[en]=intersect(que[en-1],que[en]);
    }
    while(st<en&&onright(p[en],que[st]))--en;
    if(en-st<=1)return 0;
    p[st]=intersect(que[st],que[en]);
    return 1;
}
double calc(vct *a,int n)
{
    double res=0;
    for(int i=1; i<n; i++)
        res+=cross(a[i],a[i+1]);
    res+=cross(a[n],a[1]);
    if(fabs(res/=2)<=eps)res=0;
    return res;
}
signed main()
{
    int n,o=0,oo=0;
    scanf("%lld",&n);
    while(n--)
    {
        double x1,x2,y1,y2;
        scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
        if(x1==x2&&y1==y2)continue;
        a[++oo].mk(x1,y1,x2,y2);
    }
    a[++oo].mk(0,0,10000,0);
    a[++oo].mk(10000,0,10000,10000);
    a[++oo].mk(10000,10000,0,10000);
    a[++oo].mk(0,10000,0,0);
    sort(a+1,a+1+oo);
    a[0].k=a[1].k-1;
    for(int i=1; i<=oo; i++)
        if(i==1||dcmp(a[i-1].k-a[i].k))
            a[++o]=a[i];
    if(!halfplane(o))puts("0.0");
    else printf("%.1lf\n",calc(p+st-1,en-st+1));
    return 0;
}

旋转卡壳

时间复杂度 \(O(n\log n)\)

其实是二维凸包的复杂度,查找只有 \(O(n)\)

注意下面这份代码写的是距离的平方

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF ((int)0x3f3f3f3f3f3f3f3f)
#define inf ((int)0xc0c0c0c0c0c0c0c0)
#define N (int)(5e4+15)
struct vct
{
    int x,y;
    
}p[N],s[N];
vct operator-(vct a,vct b){return {a.x-b.x,a.y-b.y};}
int cross(vct a,vct b)
{
    return a.x*b.y-a.y*b.x;
}
#define pf(x) ((x)*(x))
int dis(vct a,vct b)
{
    return pf(a.x-b.x)+pf(a.y-b.y);
}
int sq(int a,int b,int c)
{
    return abs(cross(p[b]-p[a],p[c]-p[b]));
}
int n,used[N],stk[N],top;
int getmx()
{
    int a=2,mx=-INF;
    if(top<4)
        return dis(s[1],s[2]);
    for(int i=2; i<=top; i++)
    {
        while(cross(s[a%top+1]-s[i],s[i-1]-s[i])>=cross(s[a]-s[i],s[i-1]-s[i]))
            a=a%top+1;
        mx=max(mx,max(dis(s[i-1],s[a]),dis(s[i],s[a])));
    }
    return mx;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin >> n;
    for(int i=1; i<=n; i++)
        cin >> p[i].x >> p[i].y;
    sort(p+1,p+1+n,[](vct a,vct b)
    {
        return a.x==b.x?a.y<b.y:a.x<b.x;
    });
    stk[++top]=1;
    for(int i=2; i<=n; i++)
    {
        while(top>1&&cross(p[stk[top]]-p[stk[top-1]],p[i]-p[stk[top]])<=0)
            used[stk[top--]]=0;
        used[i]=1;
        stk[++top]=i;
    }
    int tmp=top;
    for(int i=n-1; i>=1; i--)
    {
        if(used[i])continue;
        while(top>tmp&&cross(p[stk[top]]-p[stk[top-1]],p[i]-p[stk[top]])<=0)
            used[stk[top--]]=0;
        used[i]=1;
        stk[++top]=i;
    }
    for(int i=1; i<=top; i++)
        s[i]=p[stk[i]];
    cout << getmx() << endl; 
    return 0;
}

扫描线

扫描线 面积并

P5490 【模板】扫描线

时间复杂度 \(O(n\log n)\)

#include <bits/stdc++.h>
using namespace std;
#define int long long
// #define double long double
#define INF ((int)0x3f3f3f3f3f3f3f3f)
#define inf ((int)0xc0c0c0c0c0c0c0c0)
#define lb lower_bound
#define N (int)(2e5+15)
int n,x[N];
struct sline
{
    int l,r,h,mark;
}s[N];
bool operator<(sline a,sline b){return a.h<b.h;}
struct node
{
    int l,r,sum,len;
//  sum: 被完全覆盖的次数
//  len: 区间内被截的长度
}t[N<<2];
#define ls(at) (at<<1)
#define rs(at) (at<<1|1)
void build(int l,int r,int at)
{
    t[at].l=l;t[at].r=r;
    if(l==r)return;
    int mid=(l+r)>>1;
    build(l,mid,ls(at));
    build(mid+1,r,rs(at));
}
void push_up(int at)
{
    int l=t[at].l,r=t[at].r;
    if(t[at].sum)
        t[at].len=x[r+1]-x[l];
    else if(l==r)t[at].len=0;
    else t[at].len=t[ls(at)].len+t[rs(at)].len;
}
void update(int nl,int nr,int c,int at)
{
    int l=t[at].l,r=t[at].r;
    if(nl<=l&&r<=nr)
    {
        t[at].sum+=c;
        push_up(at);
        return;
    }
    int mid=(l+r)>>1;
    if(nl<=mid)update(nl,nr,c,ls(at));
    if(nr>mid)update(nl,nr,c,rs(at));
    push_up(at);
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin >> n;
    for(int i=1,x1,x2,y1,y2; i<=n; i++)
    {
        cin >> x1 >> y1 >> x2 >> y2;
        x[2*i-1]=x1;x[2*i]=x2;
        s[2*i-1]={x1,x2,y1,1};
        s[2*i]={x1,x2,y2,-1};
    }
    n<<=1;
    sort(s+1,s+1+n);
    sort(x+1,x+1+n);
    int pos=unique(x+1,x+1+n)-x-1;
    for(int i=1; i<=n; i++)
    {
        s[i].l=lb(x+1,x+1+pos,s[i].l)-x;
        s[i].r=lb(x+1,x+1+pos,s[i].r)-x-1;
    }
    build(1,pos,1);
    int res=0;
    for(int i=1; i<=n; i++)
    {
        update(s[i].l,s[i].r,s[i].mark,1);
        res+=t[1].len*(s[i+1].h-s[i].h);
    }
    cout << res << endl;
    return 0;
}

扫描线 周长并(并的周长)

P1856 [IOI1998] [USACO5.5] 矩形周长Picture

【旧代码,懒得维护】

跑了两边扫描线(纵、横)

注意不要忘记把最后一条线搞一下

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF ((int)0x3f3f3f3f3f3f3f3f)
#define inf ((int)0xc0c0c0c0c0c0c0c0)
#define lb lower_bound
#define N (int)(1e4+15)
int n,ans;
struct process
{
    struct sline
    {
        int l,r,h,mark;
        bool operator<(const sline &o)const
        {
            if(h==o.h)
                return mark>o.mark;
            return h<o.h;
        }
    }s[N];
    struct node
    {
        int l,r,sum,len;
    }t[N<<2];
    int x[N],pre;
    #define ls(at) (at<<1)
    #define rs(at) (at<<1|1)
    void build(int l,int r,int at)
    {
        t[at].l=l;t[at].r=r;
        if(l==r)return;
        int mid=(l+r)>>1;
        build(l,mid,ls(at));
        build(mid+1,r,rs(at));
    }
    void push_up(int at)
    {
        int l=t[at].l,r=t[at].r;
        if(t[at].sum)
            t[at].len=x[r+1]-x[l];
        else if(l==r)t[at].len=0;
        else t[at].len=t[ls(at)].len+t[rs(at)].len;
    }
    void update(int nl,int nr,int c,int at)
    {
        int l=t[at].l,r=t[at].r;
        if(nl<=l&&r<=nr)
        {
            t[at].sum+=c;
            push_up(at);
            return;
        }
        int mid=(l+r)>>1;
        if(nl<=mid)update(nl,nr,c,ls(at));
        if(nr>mid)update(nl,nr,c,rs(at));
        push_up(at);
    }
    void solve()
    {
        sort(s+1,s+1+n);
        sort(x+1,x+1+n);
        int pos=unique(x+1,x+1+n)-x-1;
        for(int i=1; i<=n; i++)
        {
            s[i].l=lb(x+1,x+1+pos,s[i].l)-x;
            s[i].r=lb(x+1,x+1+pos,s[i].r)-x-1;
        }
        build(1,pos,1);
        for(int i=1; i<n; i++)
        {
            update(s[i].l,s[i].r,s[i].mark,1);
            ans+=abs(t[1].len-pre);
            pre=t[1].len;
        }
        ans+=x[s[n].r+1]-x[s[n].l];
    }
}p1,p2;
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin >> n;
    for(int i=1,x1,x2,y1,y2; i<=n; i++)
    {
        cin >> x1 >> y1 >> x2 >> y2;
        p1.x[2*i-1]=x1;p1.x[2*i]=x2;
        p1.s[2*i-1]={x1,x2,y1,1};
        p1.s[2*i]={x1,x2,y2,-1};

        p2.x[2*i-1]=y1;p2.x[2*i]=y2;
        p2.s[2*i-1]={y1,y2,x1,1};
        p2.s[2*i]={y1,y2,x2,-1};
    }
    n<<=1;
    p1.solve();p2.solve();
    cout << ans << endl;
    return 0;
}

随机增量法

时间复杂度 \(O(n)\)

三点确定一个圆,故三层循环

证: \[ \begin{aligned} T_1(n) &= O(n) + \sum\limits_{i=1}^{n}\dfrac{3}{i}T_2(i) \\ T_2(n) &= O(n) + \sum\limits_{i=1}^{n}\dfrac{3}{i}T_3(i) \\ T_3(n) &= O(n)\\\\ \therefore T_1(n)&=T_2(n)=T_3(n)=O(n) \end{aligned} \] 例题:P1742 最小圆覆盖

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF ((int)0x3f3f3f3f3f3f3f3f)
#define inf ((int)0xc0c0c0c0c0c0c0c0)
#define N (int)(1e5+15)
#define double long double
int n;
const double eps=1e-10;
#define pf(x) ((x)*(x))
struct vct{double x,y;}p[N],C;
vct operator+(vct a,vct b){return {a.x+b.x,a.y+b.y};}
vct operator-(vct a,vct b){return {a.x-b.x,a.y-b.y};}
vct operator*(vct a,double b){return {a.x*b,a.y*b};}
vct operator/(vct a,double b){return {a.x/b,a.y/b};}
double cross(vct a,vct b){return a.x*b.y-a.y*b.x;}
double dot(vct a,vct b){return a.x*b.x+a.y*b.y;}
double dis(vct a,vct b){return sqrt(pf(a.x-b.x)+pf(a.y-b.y));}
double len(vct a){return sqrt(pf(a.x)+pf(a.y));}
vct rot(vct a){return {-a.y,a.x};}
struct line
{
    vct p,way;
    void mkline(vct a,vct b){p=a;way=b;}
    void mk(double x1,double y1,double x2,double y2)
    {mkline({x1,y1},{x2-x1,y2-y1});}
};
vct intersect(line a,line b)
{
    double t=cross(b.way,a.p-b.p)/cross(a.way,b.way);
    return a.p+a.way*t;
}
vct getcircle(vct a,vct b,vct c)
{
    line p,q;
    p.mkline((a+b)/2,rot(b-a));
    q.mkline((a+c)/2,rot(c-a));
    return intersect(p,q);
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);vct p=(a.p+b.p)/2;
    cout << fixed << setprecision(10);
    mt19937_64 rd(time(0));
    cin >> n;
    for(int i=1; i<=n; i++)
        cin >> p[i].x >> p[i].y;
    shuffle(p+1,p+1+n,rd);
    double r=0;
    for(int i=1; i<=n; i++) if(len(C-p[i])>r)
    {
        C=p[i];r=0;
        
        for(int j=1; j<i; j++) if(len(C-p[j])>r)
        {
            C=(p[i]+p[j])/2;
            r=len(C-p[j]);
            
            for(int k=1; k<j; k++) if(len(C-p[k])>r)
            {
                C=getcircle(p[i],p[j],p[k]);
                r=len(C-p[k]);
            }
        }
    }
    cout << r << endl;
    cout << C.x << " " << C.y << endl;
    return 0;
}

文章作者: q779
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-ND 4.0 许可协议。转载请注明来源 q779 !
评论
  目录