[模板] 计算几何2: 自适应Simpson/凸包/半平面交/旋转卡壳/闵可夫斯基和

痛定思痛。 2023-08-17 17:16 65阅读 0赞

一些基本的定义在这里: [\[模板\] 计算几何1(基础): 点/向量/线/圆/多边形/其他运算][_1_]

## 自适应Simpson ##

Simpson's Rule:

\\\[ \\int ^b\_a f(x)dx\\approx \\frac\{b-a\}6(f(a)+4f(\\frac\{a+b\}\{2\})+f(b)) \\\]

这是对二次函数的积分估值, 对于一, 二次函数来说都是准确的.

但是对于其他函数来说, 这只是利用二次函数进行近似.

可以采用自适应精度的手段, 使得估值接近真实结果. 详见代码.

然后这是误差估计, 详见 [adaptive.pdf][]:

\\\[ \\int\_\{a\}^b f(x) \\mathrm\{d\}x = S(a, c) + S(c, b) + \\frac\{1\}\{15\}\[S(a, c) + S(b, c) - S(a, b)\] + O((b - a)^6) \\\]

一种实现:

db f(db x){
        //returns f(x)
    }
    
    db simp(db l,db r){
        db mid=(l+r)/2.0;
        return (r-l)*(f(l)+4*f(mid)+f(r))/6.0;
    }
    db asr(db l,db r,db ans){
        db mid=(l+r)/2.0;
        db vl=simp(l,mid),vr=simp(mid,r),tmp=vl+vr-ans;
        if(fabs(tmp)<=eps)return ans;
        else return asr(l,mid,vl)+asr(mid,r,vr);
    }

## 凸包 ##

Andrew 算法, 即分别求上, 下凸包. 时间复杂度 \\(O(n \\log n)\\).

struct tvec{db x,y;};
    
    il int dcmp(db a){return fabs(a)<=eps?0:(a>0?1:-1);}
    il db p2(db a){return a*a;}
    il db gougu1(db a,db b){return sqrt(p2(a)+p2(b));}
    il tvec operator+(tvec a,tvec b){return (tvec){a.x+b.x,a.y+b.y};}
    il tvec operator-(tvec a,tvec b){return (tvec){a.x-b.x,a.y-b.y};}
    il tvec operator*(tvec a,db b){return (tvec){a.x*b,a.y*b};}
    il tvec operator*(db a,tvec b){return b*a;}
    il db operator*(tvec a,tvec b){return a.x*b.y-b.x*a.y;}
    il db operator^(tvec a,tvec b){return a.x*b.x+a.y*b.y;}
    il db len(tvec a){return gougu1(a.x,a.y);}
    bool cmp(tvec a,tvec b){int tmp=dcmp(a.x-b.x);return tmp?tmp<0:dcmp(a.y-b.y)<0;}
    
    tvec li[nsz],conv[nsz];
    int pc=0;
    void getconv(){
        sort(li+1,li+n+1,cmp);
        rep(i,1,n){
            while(pc>1&&dcmp((conv[pc]-conv[pc-1])*(li[i]-conv[pc]))<=0)--pc;
            conv[++pc]=li[i];
        }
        int tmp=pc;
        repdo(i,n-1,1){
            while(pc>tmp&&dcmp((conv[pc]-conv[pc-1])*(li[i]-conv[pc]))<=0)--pc;
            conv[++pc]=li[i];
        }
        if(n>1)--pc;
    }

## 半平面交 ##

增量法, 时间复杂度 \\(O(n \\log n)\\) (排序的复杂度).

需要保证不是开放的半平面. 否则加上四个 \\(\\pm \\infty\\) 的平面即可.

细节较多. 详见代码...

const int psz=550;
    const db eps=1e-9;
    int n,m;
    
    db dcmp(db v){return fabs(v)<=eps?0:(v>0?1:-1);}
    db p2(db v){return v*v;}
    
    struct tvec{db x,y;};
    tvec operator+(tvec a,tvec b){return (tvec){a.x+b.x,a.y+b.y};}
    tvec operator-(tvec a,tvec b){return (tvec){a.x-b.x,a.y-b.y};}
    tvec operator*(tvec a,db b){return (tvec){a.x*b,a.y*b};}
    tvec operator*(db a,tvec b){return b*a;}
    db operator*(tvec a,tvec b){return a.x*b.y-a.y*b.x;}
    db operator^(tvec a,tvec b){return a.x*b.x+a.y*b.y;}
    db len(tvec a){return sqrt(p2(a.x)+p2(a.y));}
    
    struct tl{
        tvec p,v;
        db d;
        tl(){}
        tl(tvec a,tvec b):p(a),v(b-a){d=atan2(v.y,v.x);}
    }li[psz];
    int pl=0;
    bool operator<(tl a,tl b){return a.d<b.d;}
    bool isleft(tl a,tvec b){return dcmp(a.v*(b-a.p))>0;}
    il tvec inters(tl a,tl b){db v=(b.v*(a.p-b.p))/(a.v*b.v);return a.p+a.v*v;}
    
    tvec poly[psz];
    int ppo=0;
    
    tl que[psz]; //queue
    tvec qp[psz]; //intersect points
    int qh=1,qt=0;
    int hplane(){//0 fail, 1 success
        sort(li+1,li+pl+1);
        int pl1=1;//suppose that pl>=1
        rep(i,2,pl){
            if(li[i].d>li[pl1].d)li[++pl1]=li[i];
            else if(isleft(li[pl1],li[i].p))li[pl1]=li[i];
        }
        pl=pl1;
        qh=1,qt=0;
        rep(i,1,pl){
            while(qh<qt&&!isleft(li[i],qp[qt-1]))--qt;
            while(qh<qt&&!isleft(li[i],qp[qh]))++qh;
            que[++qt]=li[i];
            if(qh<qt)qp[qt-1]=inters(que[qt-1],que[qt]);
        }
        while(qh<qt&&!isleft(que[qh],qp[qt-1]))--qt; //**
        ppo=0;
        if(qt-qh<=1)return 0; //no sol
        qp[qt]=inters(que[qh],que[qt]);
        rep(i,qh,qt)poly[++ppo]=qp[i];
        return 1;
    }

## 旋转卡壳 ##

这是一种拥有 \\(4\\) 个多音字, \\(2^4 = 16\\) 种读音的优秀算法.

## 闵可夫斯基和 ##

转载于:https://www.cnblogs.com/ubospica/p/10828219.html

[_1_]: https://www.cnblogs.com/ubospica/p/10121673.html
[adaptive.pdf]: http://www2.math.umd.edu/~mariakc/teaching/adaptive.pdf