计算几何小结

好像又鸽子了几天博客

poj真是神奇的网站……

基础知识

点积

\(a·b=a.x*b.x+a.y*b.y\)

\(a\)在\(b\)上的投影乘以\(b\)的模长

叉积

\(a×b=a.x*b.y-a.y*b.x\)

\(a,b\)围成的平行四边形的有向面积

判断线段相交（跨立实验）

由一条线段的端点向另一条线段两端点做叉积，判断符号是否一致

向量旋转

逆时针旋转\(k\)度：

inline vector turn (vector a,double k)
{
    return vector(a.x*cos(k)-a.y*sin(k),a.x*sin(k)+a.y*cos(k));
}

判断点在多边形内部

由一点引一条射线，若与多边形有奇数个交点则在多边形内部，否则在多边形外

二维凸包

求凸包

二维平面上给定\(n\)个点求凸包

显然\(y\)最小的点应该在凸包上

我们把\(y\)最小的点先选出来，将其他点极角排序

按照极角序将点入栈，利用斜率单调性判断是否弹栈

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
namespace red{
#define y1 qwq
#define int long long
#define eps (1e-10)
    inline int read()
    {
        int x=0;char ch,f=1;
        for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
        if(ch=='-') f=0,ch=getchar();
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
        return f?x:-x;
    }
    const int N=1e5+10;
    int n;
    int st[N],top;
    double ret;
    struct node
    {
        double x,y;
        node (double tx=0,double ty=0){x=tx,y=ty;}
        inline double operator ^ (const node &t) const
        {
            return x*t.y-y*t.x;
        }
        inline node operator - (const node &t) const
        {
            return node(x-t.x,y-t.y);
        }
    }a[N];
    inline double sqr(double x){return x*x;}
    inline double dis(node a,node b)
    {
        return sqrt(sqr(b.x-a.x)+sqr(b.y-a.y));
    }
    inline bool cmp1(node a,node b)
    {
        return a.x==b.x?a.y<b.y:a.x<b.x;
    }
    inline bool cmp2(node q,node b)
    {
        double tmp=(q-a[1])^(b-a[1]);
        if(!tmp) return q.x==b.x?q.y>b.y:q.x<b.x;
        return tmp>0;
    }
    inline void main()
    {
        n=read();
        for(int i=1;i<=n;++i)
        {
            scanf("%lf%lf",&a[i].x,&a[i].y);
        }
        sort(a+1,a+n+1,cmp1);
        int tot=0;
        a[0].x=-1e9+7;
        for(int i=1;i<=n;++i)
        {
            if(a[i].x!=a[tot].x||a[i].y!=a[tot].y) a[++tot]=a[i];
        }
        sort(a+2,a+tot+1,cmp2);
        st[++top]=1;st[++top]=2;
        for(int i=3;i<=tot;++i)
        {
            while(top>2&&((a[st[top-1]]-a[st[top]])^(a[i]-a[st[top]]))>=0) --top;
            st[++top]=i;
        }
        for(int i=1;i<top;++i)
        {
            ret+=dis(a[st[i]],a[st[i+1]]);
        }
        ret+=dis(a[st[top]],a[st[1]]);
        printf("%.2lf",ret);
    }
}
signed main()
{
    red::main();
    return 0;
}

二维凸包面积

利用叉积

inline double cross(node a,node b,node c)
{
    return (b-a)^(c-a);//叉积
}
inline double are()
{
    double ret=0;
    for(int i=1;i<n;++i) ret+=fabs(cross(a[1],a[i],a[i+1]));
    return ret*2;
}

多次判断点在凸包内部

将凸包分成\(n-2\)个三角形，每次取中间的三角形，用叉积判断当前点在三角形的内部还是左侧或右侧，二分求解

没写过

动态二维凸包

初始有三个点，每次加入一个点，维护凸包

由于凸包只会扩大，所以以初始的三点围成的三角形的重心为基准点求其他点的极角，将点插入合适的位置，再判断两侧的点是否删除

用平衡树/\(set\)维护

注意每次插入一个点不一定只弹出相邻的两个点……

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
using namespace std;
namespace red{
#define y1 qwq
#define int long long
#define double long double
#define iter set<node>::iterator
#define eps (1e-8)
    inline int read()
    {
        int x=0;char ch,f=1;
        for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
        if(ch=='-') f=0,ch=getchar();
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
        return f?x:-x;
    }
    const int N=1e5+10;
    int n,opt;
    struct node
    {
        double x,y,ang;
        node(double tx=0,double ty=0,double tang=0){x=tx,y=ty,ang=tang;}
        friend bool operator < (node a,node b)
        {
            return (a.ang<b.ang)||(a.ang==b.ang&&a.x<b.x);
        }
        inline node operator + (const node b) const
        {
            return node(x+b.x,y+b.y,0);
        }
        inline double operator ^ (const node &b) const
        {
            return x*b.y-y*b.x;
        }
        inline node operator - (const node &b) const
        {
            return node(x-b.x,y-b.y,0);
        }
    }a[N],zx;
    set<node> q;
    iter it1,it2;
    iter pre(iter x)
    {
        return x==q.begin()?--q.end():--x;
    }
    iter nxt(iter x)
    {
        return ++x==q.end()?q.begin():x;
    }
    inline void main()
    {
        n=read();
        for(int i=1;i<=3;++i)
        {
            opt=read();
            a[i].x=read(),a[i].y=read();
            zx=zx+a[i];
            a[i].x*=3,a[i].y*=3;
        }
        for(int i=1;i<=3;++i)
        {
            a[i].ang=atan2(a[i].y-zx.y,a[i].x-zx.x);
            q.insert(a[i]);
        }
        for(int i=4;i<=n;++i)
        {
            opt=read();
            a[i].x=read()*3,a[i].y=read()*3;
            a[i].ang=atan2(a[i].y-zx.y,a[i].x-zx.x);
            it2=q.lower_bound(a[i]);
            if(it2==q.end()) it2=q.begin();
            it1=pre(it2);
            if(opt==1)
            {
                if(((*it2-a[i])^(*it1-a[i]))<=0) continue;
                q.insert(a[i]);
                iter now=q.find(a[i]);
                iter tx=pre(now),ty=pre(tx);
                
                while(((*ty-*now)^(*tx-*now))<=0) q.erase(tx),tx=ty,ty=pre(tx);
                tx=nxt(now),ty=nxt(tx);
                while(((*ty-*now)^(*tx-*now))>=0) q.erase(tx),tx=ty,ty=nxt(tx);
            }
            else
            {
                if(((*it2-a[i])^(*it1-a[i]))<=0) puts("YES");
                else puts("NO");
            }
        }
    }
}
signed main()
{
    red::main();
    return 0;
}

旋转卡壳

组合数学入门题

众所周知，旋转卡壳有\(2*3*2*2=24\)种读法（雾

其实就是单调性优化的枚举？

给定凸包求最大直径

然后我们发现枚举一个点的话另一个点有单调性，所以这个过程是\(O(n)\)的……

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
namespace red{
#define y1 qwq
#define int long long
#define eps (1e-10)
    inline int read()
    {
        int x=0;char ch,f=1;
        for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
        if(ch=='-') f=0,ch=getchar();
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
        return f?x:-x;
    }
    const int N=1e5+10;
    int n;
    int st[N],top;
    double ret;
    struct node
    {
        double x,y;
        node (double tx=0,double ty=0){x=tx,y=ty;}
        inline double operator ^ (const node &t) const
        {
            return x*t.y-y*t.x;
        }
        inline node operator - (const node &t) const
        {
            return node(x-t.x,y-t.y);
        }
    }a[N];
    inline double sqr(double x){return x*x;}
    inline double dis(node a,node b)
    {
        return sqr(b.x-a.x)+sqr(b.y-a.y);
    }
    inline bool cmp1(node a,node b)
    {
        return a.x==b.x?a.y<b.y:a.x<b.x;
    }
    inline bool cmp2(node q,node b)
    {
        double tmp=(q-a[1])^(b-a[1]);
        if(!tmp) return q.x==b.x?q.y>b.y:q.x<b.x;
        return tmp>0;
    }
    inline double are(node a,node b,node c)
    {
        return fabs((a.x*b.y+b.x*c.y+c.x*a.y-a.x*c.y-b.x*a.y-c.x*b.y)/2);
    }
    inline void main()
    {
        n=read();
        for(int i=1;i<=n;++i)
        {
            scanf("%lf%lf",&a[i].x,&a[i].y);
        }
        sort(a+1,a+n+1,cmp1);
        int tot=0;
        a[0].x=-1e9+7;
        for(int i=1;i<=n;++i)
        {
            if(a[i].x!=a[tot].x||a[i].y!=a[tot].y) a[++tot]=a[i];
        }
        sort(a+2,a+tot+1,cmp2);
        st[++top]=1;st[++top]=2;
        for(int i=3;i<=tot;++i)
        {
            while(top>2&&((a[st[top-1]]-a[st[top]])^(a[i]-a[st[top]]))>=0) --top;
            st[++top]=i;
        }
        tot=top;
        for(int i=1;i<=tot;++i) st[++top]=st[i];
        int tmp=3;
        for(int i=1;i<=tot;++i)
        {
            while(are(a[st[i]],a[st[i+1]],a[st[tmp+1]])>are(a[st[i]],a[st[i+1]],a[st[tmp]])) ++tmp;
            ret=max(ret,max(dis(a[st[i]],a[st[tmp]]),dis(a[st[i+1]],a[st[tmp]])));
        }
        printf("%lld",(int)ret);
    }
}
signed main()
{
    red::main();
    return 0;
}

另一个经典应用：给定两凸包求两凸包上点的最近距离

取其中一个凸包的\(y\)最小的点，另一凸包\(y\)最大的点开始枚举，方向相反，这样就有单调性了

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
namespace red{
#define int long long
#define ls(p) (p<<1)
#define rs(p) (p<<1|1)
#define eps (1e-9)
    inline int read()
    {
        int x=0;char ch,f=1;
        for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
        if(ch=='-') f=0,ch=getchar();
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
        return f?x:-x;
    }
    const int N=5e4+10;
    int n,m;
    struct node
    {
        double x,y;
        node() {};
        node (double _x,double _y):x(_x),y(_y){}
        inline double operator ^ (const node &t) const {return (x*t.y-y*t.x);}
        inline double operator * (const node &t) const {return (x*t.x+y*t.y);}
        inline node operator - (const node &t) const {return node(x-t.x,y-t.y);}
        inline bool operator < (const node &b) const {return atan2(y,x)<atan2(b.y,b.x);}
    }a[N],b[N];
    inline double sqr(double x) {return x*x;}
    inline double dist(node a,node b)
    {
        return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
    }
    double multi(node A,node B,node C)//点积
    {
        return (B-A)*(C-A);
    }
    double cross(node A,node B,node C)//叉积
    {
        return (B-A)^(C-A);
    }
    inline double disline(node A,node B,node C)
    {
        if(dist(A,B)<eps) return dist(B,C);
        if(multi(A,B,C)<-eps) return dist(A,C);//用点积判断是否在线段内部
        if(multi(B,A,C)<-eps) return dist(B,C);
        return fabs(cross(A,B,C)/dist(A,B));//点到线段距离
    }
    double mindist(node A,node B,node C,node D)
    {
        return min(min(disline(A,B,C),disline(A,B,D)),min(disline(C,D,A),disline(C,D,B)));
    }
    inline void check(node a[],int n)
    {
        for(int i=1;i<n-1;++i)
        {
            double tmp=cross(a[i],a[i+1],a[i+2]);
            if(tmp<-eps) return;
            else if(tmp>eps)
            {
                reverse(a+1,a+n+1);
                return;
            }
        }
    }
    inline double solve(node a[],node b[],int n,int m)
    {
        int ymina=1,ymaxb=1;
        for(int i=1;i<=n;++i) if(a[i].y<a[ymina].y) ymina=i;
        for(int i=1;i<=m;++i) if(a[i].y>a[ymaxb].y) ymaxb=i;
        double tmp,ret=1e9+7;
        a[n+1]=a[1],b[m+1]=b[1];
        for(int i=1;i<=n;++i)
        {
            while(tmp=cross(a[ymina+1],b[ymaxb+1],a[ymina])-cross(a[ymina+1],b[ymaxb],a[ymina])<eps)
            {
                ++ymaxb;
                if(ymaxb>m) ymaxb=1;
            }
            ret=min(ret,mindist(a[ymina],a[ymina+1],b[ymaxb],b[ymaxb+1]));
            ++ymina;
            if(ymina>n) ymina=1;
        }
        return ret;
    }
    inline void main()
    {
        while("haku")
        {
            n=read(),m=read();
            if(!(n|m)) return ;
            for(int i=1;i<=n;++i) scanf("%lf%lf",&a[i].x,&a[i].y);
            for(int i=1;i<=m;++i) scanf("%lf%lf",&b[i].x,&b[i].y);
            check(a,n);check(b,m);
            printf("%.5lf ",min(solve(a,b,n,m),solve(b,a,m,n)));
        }
    }
}
signed main()
{
    red::main();
    return 0;
}

半平面交

有点恶心的东西

给定\(n\)条直线，每条直线保留左侧部分，求最后留下的面积

首先我们有\(n^2\)的暴力算法：每加入一条直线，就枚举以前的直线

如果在左侧则保留

如果在右侧把它扔掉

如果相交保留左侧端点，右侧改为交点

当然一般来说\(n^2\)是不够用的

我们可以先将所有直线按照与\(x\)轴的夹角排序，夹角相同的保留靠左的

用双端队列维护半平面交

维护方式：先保证队列中至少有两个元素

如果队尾的两个元素的交点在当前直线右侧，把队尾弹出

队首同理

但是记住一定要先弹出队尾，因为新加入的元素与一开始加入的元素围成的面积更小

如在这张图中，如果标号就是加入队列的顺序，那么当\(3\)加入的时候，先判断\(1\)就会踢\(1\)，先判断\(2\)就会踢\(2\)，但是新加入的元素和队首围成的面积会更小，所以应该先踢队尾

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define ls(p) (p<<1)
#define rs(p) (p<<1|1)
#define eps (1e-8)
    inline int read()
    {
        int x=0;char ch,f=1;
        for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
        if(ch=='-') f=0,ch=getchar();
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
        return f?x:-x;
    }
    const int N=5000;
    int n,m,tot,sum,head,tail;
    struct node
    {
        double x,y;
        inline double operator ^ (const node &t) const{return x*t.y-y*t.x;}
        inline node operator - (const node &t) const{return (node){x-t.x,y-t.y};}
        inline node operator + (const node &t) const{return (node){x+t.x,y+t.y};}
        inline node operator * (const double &t) const{return (node){x*t,y*t};}
    }a[N],c[N];
    struct seg
    {
        node a,b;
        double k;
        seg(){}
        seg(const node &aa,const node &bb):a(aa),b(bb){k=atan2(b.y,b.x);}
        inline bool operator < (const seg &t) const
        {
            return k<t.k;
        }
    }q[N],que[N];
    inline double cross(node a,node b,node c)
    {
        return (b-a)^(c-a);
    }
    node get_node(seg A,seg B)//求交点
    {
        node C=A.a-B.a;
        double t=(B.b^C)/(A.b^B.b);
        return A.a+A.b*t; 
    }
    inline int dcmp(double x)
    {
        if(fabs(x)<eps) return 0;
        return x>0?1:-1;
    }
    inline double are()
    {
        double are=0;
        for(int i=head;i<tail;++i){ are+=fabs(cross(c[head],c[i],c[i+1])); }
        return are/2;
    }
    inline void work()
    {
        head=tail=1;
        que[tail]=q[1];
        for(int i=2;i<=sum;++i)
        {
            while(head<tail&&(q[i].b^(c[tail-1]-q[i].a))<=eps) --tail;
            while(head<tail&&(q[i].b^(c[head]-q[i].a))<=eps) ++head;
            que[++tail]=q[i];
            if(fabs(que[tail].b^que[tail-1].b)<=eps) //平行保留较左的
            {
                --tail;
                if((que[tail].b^(q[i].a-que[tail].a))>eps) que[tail]=q[i];
            }
            if(head<tail) c[tail-1]=get_node(que[tail-1],que[tail]);//c[i]表示c[i]和c[i+1]的交点
        }
        while(head<tail&&(que[head].b^(c[tail-1]-que[head].a))<=eps) --tail;
        if(tail-head<=1) return;
        c[tail]=get_node(que[head],que[tail]);
    }
    inline void main()
    {
        n=read();
        for(int i=1;i<=n;++i)
        {
            m=read();
            for(int j=1;j<=m;++j)
            {
                a[j].x=read(),a[j].y=read();
            }
            for(int j=1;j<=m;++j)
            {
                ++sum;
                int t=j+1;
                if(j==m) t=1;
                q[sum]=seg(a[j],a[t]-a[j]);
            }
        }
        sort(q+1,q+sum+1);
        work();
        printf("%.3lf",are());
    }
}
signed main()
{
    red::main();
    return 0;
}

自适应辛普森积分法

\(simple\)积分法

求积分\(\int_{l}^{r}\frac{cx+d}{ax+b}dx\)

先来说什么是辛普森积分法

假设有函数\(g(x)=Ax^2+Bx+C\)

\[\int_{a}^{b}g(x)dx\]

微积分基本定理：

\[=\int_{a}^{b}\frac{A}{3}(b^3-a^3)+\frac{B}{2}(b^2-a^2)+C(b-a)\]

大力化简

\[=\frac{A}{3}(b-a)(a^2+ab+b^2)+\frac{B}{2}(b+a)(b-a)+C(b-a)\]

\[=\frac{b-a}{6}(2A(a^2+ab+b^2)+3B(b+a)+6C)\]

\[=\frac{b-a}{6}(2Aa^2+2Aab+2Ab^2+3Ba+3Bb+6C)\]

\[=\frac{b-a}{6}(Aa^2+Ba+C+Ab^2+Bb+C+4A(\frac{a+b}{2})^2+4B(\frac{a+b}{2})+4C)\]

\[=\frac{b-a}{6}(g(a)+g(b)+4g(\frac{a+b}{2}))\]

得到辛普森积分公式

\[\int_{a}^{b}g(x)dx=\frac{b-a}{6}(g(a)+g(b)+4g(\frac{a+b}{2}))\]

不过这是二次函数的积分公式，和我们上面那个式子有啥关系？

我们知道，一个复杂的函数可以近似拟合为一个二次函数

如果\(f(x)\)的拟合二次函数为\(g(x)\)，那么

\[\int_{a}^{b}f(x)dx≈\frac{b-a}{6}(f(a)+f(b)+4f(\frac{a+b}{2}))\]

当然拟合肯定存在误差，不过当定义域越小的时候，这个误差也就越小

我们可以把原函数分段，对于每一段求拟合函数\(g(x)\)的值再加起来，只要分得的段足够小就可以得到正确答案

但是分的太小会导致答案不正确， 太多会降低程序效率

这个时候我们可以让程序自适应

二分递归答案的段数和区间长度，精度达到要求就停止递归返回答案

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<set>
using namespace std;
namespace red{
#define y1 qwq
    inline int read()
    {
        int x=0;char ch,f=1;
        for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
        if(ch=='-') f=0,ch=getchar();
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
        return f?x:-x;
    }
    double a,b,c,d,l,r;
    inline double f(double x)
    {
        return (c*x+d)/(a*x+b);
    }
    inline double simpson(double l,double r)
    {
        double mid=(l+r)/2;
        return (f(l)+4*f(mid)+f(r))*(r-l)/6;
    }
    inline double asr(double l,double r,double eps,double ans)
    {
        double mid=(l+r)/2;
        double tl=simpson(l,mid),tr=simpson(mid,r);
        if(fabs(tl+tr-ans)<=15*eps) return tl+tr+(tl+tr-ans)/15;
        return asr(l,mid,eps/2,tl)+asr(mid,r,eps/2,tr);
    }
    inline double asr(double l,double r,double eps)
    {
        return asr(l,r,eps,simpson(l,r));
    }
    inline void main()
    {
        scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
        printf("%.6f\n",asr(l,r,1e-6));
    }
}
signed main()
{
    red::main();
    return 0;
}

计算几何小结

原文：https://www.cnblogs.com/knife-rose/p/12179026.html