计算几何

KellyWLJ · 2023-07-23 11:09:04 · 个人记录

嗯，是我讨厌且不擅长的东西）

现在好像很少直接考，但是可以利用它的性质与 DP 等东西结合，eg：2023春测 T3

做题时多画图手玩，不要把叉积方向弄反了

基础知识

向量：

通常用坐标 \vec a=(x,y) 表示

运算

加减法，数乘，点积见数学必修2

主要是叉积，\vec a\times \vec b=|\vec a||\vec b|\sin\theta，其中 \theta 为 \vec a 逆时针旋转到 \vec b 的角度

写成坐标为 (x_1,y_1)\times (x_2,y_2)=x_1y_2-x_2y_1

它的几何意义是 \vec a,\vec b 张成的平行四边形的有向面积

是判断位置关系，求面积非常好用的运算

• 当 \vec a\times \vec b=0 时，\vec a,\vec b 共线，方向相同或相反

• 当 \vec a\times \vec b>0 时，\vec a 在 \vec b 的下方

• 当 \vec a\times \vec b<0 时，\vec a 在 \vec b 的上方

注意叉积没有交换律，\vec a\times \vec b=-\vec a\times \vec b

极角排序

这里是把向量按逆时针方向排序，可以理解为排成一圈

以第三象限为起点，逆时针把四个象限标号为 1\sim 4，标号小的在前面

比较时先判断两个向量是否在 x 轴同一侧，如果不在则直接排序

如果在则有两种方式继续判断：

• 叉积，判断是否在逆时针方向，适用于整数运算，对精度要求较高的题

• 三角函数 atan2(y,x)，返回 (-\pi,\pi] 之间的值，表示以 x 轴非负半轴为始边，以向量所在射线为终边的角的弧度

向量旋转

先研究以原点为起点的

设 \vec u=(x,y)，模长为 a，则 \vec u=(a\cos\alpha,a\sin\alpha)

将 u 逆时针旋转 \theta，则 \vec u=(a\cos(\alpha+\theta),a\sin(\alpha+\theta))

其实可以这么硬算，但太麻烦，且卡精度

拆开，代入 x,y，得 \vec u=(x\cos\theta-y\sin\theta,y\cos\theta+x\sin\theta)

如果不以原点为起点，则将起点看作原点，求出新向量后旋转，还原回去即可

pdd zhuan(pdd x, pdd y, ld ang) // x绕 y逆时针旋转，角度为 ang
{
    x = x - y;
    return mp(x.fi * cosl(ang) - x.se * sinl(ang) + y.fi, x.se * cosl(ang) + x.fi * sinl(ang) + y.se);
}

求两向量的夹角

可以利用点积，叉积求出三角函数值，再逆着求回去

也可以将某个看成坐标轴，整体旋转上去，用 atan2 或只要求比大小时叉积也行

Nearest vectors

ll dot(pll x, pll y)    {return x.fi * y.fi + x.se * y.se;}
ll cross(pll x, pll y)  {return x.fi * y.se - x.se * y.fi;}
ll cmp(const ll &a, const ll &b)
{
    if(p[a].se <= 0 && p[b].se > 0) return 1;
    if(p[a].se > 0 && p[b].se <= 0) return 0;
    if(p[a].se == 0 && p[b].se == 0)    return p[a].fi < p[b].fi;
    return cross(p[a], p[b]) > 0;
}
int main()
{
    n = read();
    for(ll i = 1; i <= n; ++i)  p[i].fi = read(), p[i].se = read(), id[i] = i;
    sort(id + 1, id + n + 1, cmp);
    num = 2, id[n + 1] = id[1];
    for(ll i = 3; i <= n + 1; ++i)
    { // 把后一个旋转至与 x轴正半轴重合，叉积比较另一个的位置关系 
        pll ang1 = mp(dot(p[id[i]], p[id[i - 1]]), abs(cross(p[id[i]], p[id[i - 1]])));
        pll ang2 = mp(dot(p[id[num]], p[id[num - 1]]), abs(cross(p[id[num]], p[id[num - 1]])));
        if(cross(ang1, ang2) > 0) // 叉积式子是约去分母后的 
            num = i;
    }
    printf("%lld %lld", id[num - 1], id[num]);
    return 0;
}

直线：

有向量后，可以用点向式描述一条直线

把直线上的某点 (x,y) 记录，同时记录方向 \vec v=(1,k)

还可以直接用两个直线上的点，看作两个向量记录

不推荐用直线解析式，精度损失较大，也容易出现 k=0 等边界情况，很复杂

点与直线的位置关系

取直线上两点 A,B，直接用叉积判断 (P-A)\times (B-A) 的符号

特别的，若为 0，则作为点在直线上的判断

线段的位置关系

判断两条线段是否相交，可以用快速排斥实验和跨立实验

快速排斥实验是初步判断，判断以两条线段为对角线的四边与坐标轴平行的矩形有没有交，简单说就是线段坐标覆盖的范围如果没交就直接返回

设线段 a 的端点为 A,B，线段 b 的端点为 C,D

那么如果两线段相交，则 \overrightarrow{AB} 一定在 \overrightarrow{AC},\overrightarrow{AD} 之间

用叉积判断，(\overrightarrow{AB}\times \overrightarrow{AC})\cdot (\overrightarrow{AB}\times \overrightarrow{AD})\ge 0 且另一条线段，以 C 为起点，同理 (\overrightarrow{CD}\times \overrightarrow{CA})\cdot (\overrightarrow{CD}\times \overrightarrow{CB})\ge 0，则两条线段相交

取等是为处理端点相交的情况

好像不用快速排斥也行

U130828 【模板】计算几何1

struct line
{
    pdd x, y;
}p[N]; 
pdd operator-(const pdd &a, const pdd &b)   {return mp(a.fi - b.fi, a.se - b.se);} 
ld cross(pdd a, pdd b)  {return a.fi * b.se - a.se * b.fi;}
int kl(line a, line b)
{
    if(cross(a.y - a.x, b.x - a.x) * cross(a.y - a.x, b.y - a.x) > eps) return 2;
    if(cross(b.y - b.x, a.x - b.x) * cross(b.y - b.x, a.y - b.x) > eps) return 2;
    return 1;
}

两条直线的交点

这里已经默认用叉积排除了平行的情况

设直线 m 一点为 P，方向为 \vec v，直线 n 上一点为 Q，方向为 \vec w，\vec u=\overrightarrow{QP}

则交点为 P+\frac{w\times u}{v\times w}\cdot v

可以用相似证，求出交点离 P 的长度

pdd cross(line x, line y) // 求直线（向量）x,y的交点 
{
    pdd w = y.ed - y.st, v = x.ed - x.st;
    return x.st + ((w ^ (x.st - y.st)) / (v ^ w)) * v;
}

凸包

分为上凸壳和下凸壳，下凸壳的斜率从左到右单调递增，上凸壳的斜率单调递减

凸包是凸多边形，可以在斜率优化中运用，也有较好的性质

比如，一个点集中相距最远的两点肯定是这个点集凸包的两个端点

求凸包

两种算法，Andrew 和 Graham，老师上课讲的是 Andrew，但我个人更喜欢 Graham，因为我觉得 Andrew 分别求上下凸包再合并细节比较多，不优美）

Graham 算法流程：

• 找到最靠下且最靠左的一个点，以它为基准

• 对剩下点与它组成的向量极角排序，由于相当于都在 x 轴上方，直接用 atan2

• 先将基准点入栈，参照斜率优化时维护凸包的方法，新点入栈时把前面不合法的弹出，用叉积判断下凸壳的方法

注意舍去多点共线的情况

P2742 [USACO5.1] 圈奶牛Fencing the Cows /【模板】二维凸包

int cmp(const pdd &a, const pdd &b)
{
    double anga = atan2(a.se - p[1].se, a.fi - p[1].fi);
    double angb = atan2(b.se - p[1].se, b.fi - p[1].fi);
    if(!equ(anga, angb))    return anga < angb;
    return len(a - p[1]) < len(b - p[1]);
}
int main()
{
    ios::sync_with_stdio(false), cin.tie(0);
    cin >> n;
    for(int i = 1; i <= n; ++i) cin >> p[i].fi >> p[i].se;
    for(int i = 2; i <= n; ++i)
        if((p[i].se < p[1].se) || (equ(p[i].se, p[1].se) && p[i].fi < p[1].fi)) swap(p[i], p[1]);
    sort(p + 2, p + n + 1, cmp), stk[++top] = 1;
    for(int i = 2; i <= n; ++i)
    {
        while(top > 1 && cross(p[i] - p[stk[top]], p[stk[top]] - p[stk[top - 1]]) >= 0) --top;
        stk[++top] = i;
    }
    for(int i = 2; i <= top; ++i)   ans += len(p[stk[i]] - p[stk[i - 1]]);
    printf("%.2f", ans + len(p[stk[top]] - p[1]));
    return 0;
}

P3829 [SHOI2012] 信用卡凸包

有圆弧不好处理，考虑以四角的圆心为端点

画图可以发现，所有的圆弧的角度加起来为 2\pi

因为每个顶点处对应的圆弧可以看作凸包中这个角的外角

然后就变成了求凸包模板题，注意旋转

pdd zhuan(pdd x, pdd y, ld ang) // x绕 y逆时针旋转，角度为 ang
{
    x = x - y;
    return mp(x.fi * cosl(ang) - x.se * sinl(ang) + y.fi, x.se * cosl(ang) + x.fi * sinl(ang) + y.se);
}
int cmp(const pdd &x, const pdd &y)
{
    ld ang1 = atan2(x.se - p[1].se, x.fi - p[1].fi);
    ld ang2 = atan2(y.se - p[1].se, y.fi - p[1].fi);
    if(!equ(ang1, ang2))    return ang1 < ang2;
    return len(x - p[1]) < len(y - p[1]); 
}
int main()
{
    ios::sync_with_stdio(false), cin.tie(0);
    cin >> n >> b >> a >> r;
    nx[1] = nx[4] = a / (ld)2.0 - r, nx[2] = nx[3] = -a / (ld)2.0 + r;
    ny[1] = ny[2] = b / (ld)2.0 - r, ny[3] = ny[4] = -b / (ld)2.0 + r;
    for(int i = 1; i <= n; ++i) 
    {
        cin >> xx >> yy >> theta;
        for(int j = 1; j <= 4; ++j)
        {
            pdd nw = mp(xx + nx[j], yy + ny[j]);
            p[++cnt] = zhuan(nw, mp(xx, yy), theta);
        }
    }
    for(int i = 2; i <= cnt; ++i)
        if(p[i].se < p[1].se || (equ(p[i].se, p[1].se) && p[i].fi < p[1].fi))   swap(p[i], p[1]);
    sort(p + 2, p + cnt + 1, cmp), stk[++top] = 1;
    for(int i = 2; i <= cnt; ++i)
    {
        while(top > 1 && cross(p[i] - p[stk[top]], p[stk[top]] - p[stk[top - 1]]) >= 0) --top;
        stk[++top] = i;
    }
    for(int i = 2; i <= top; ++i)   ans += len(p[stk[i]] - p[stk[i - 1]]);
    printf("%.2Lf", ans + len(p[stk[top]] - p[1]) + pi * r * (ld)2.0);
    return 0;
}

求凸包面积

随便找一个点，逆时针排列点，把答案加上这个点和枚举的两个点组成的三角形有向面积

用叉积计算，发现这样多余部分刚好相互抵消，实现上这个点一般直接取一个顶点

for(int i = 2; i < idx; ++i)    ans += S(a[1], a[i], a[i + 1]); 
// 求这个多边形的面积

点与凸包的位置关系

边界上：直接判断是否在每条向量上

方法一：

可以以这个点为起点，随机生成一条射线，并求出与凸包的交点个数

如果交点为奇数，则通常在凸包内部，如果交点个数为偶数（包括 0），则在外部

但是，如果在外部且随机的直线刚好与凸包相切，那就完了

可以多随机几次，一般出错概率较小

方法二：

同 Graham 扫描法，以凸包最靠下，左的顶点为起点，逆时针排列凸包上的点

先通过与顶点和它两边的点形成向量的位置关系排除不可能的情况

然后二分查找出这个点与顶点相连形成的向量在哪两个顶点这样形成的向量之间

用叉积判断即可，如果共线则判断距离

注意细节较多

ll cmp(const pdd &x, const pdd &y)
{
    if(sgn((x - bas) ^ (y - bas)) == 0) return len(x - bas) < len(y - bas);
    return sgn((x - bas) ^ (y - bas)) > 0;
}
ll in(vector<pdd> &p, pdd x)
{
    if(p.size() <= 1)   return 2;
    bas = p[0]; ll pos = lower_bound(p.begin() + 1, p.end(), x, cmp) - p.begin() - 1;
    if(sgn((x - bas) ^ (p.back() - bas)) < 0 || sgn((x - bas) ^ (p[1] - bas)) > 0)  return 2;
    if(sgn((p[(pos + 1) % p.size()] - p[pos]) ^ (x - p[pos])) == 0)    
        return len(x - p[pos]) <= len(p[(pos + 1) % p.size()] - p[pos]) ? 1 : 2;
    return sgn((p[(pos + 1) % p.size()] - p[pos]) ^ (x - p[pos])) > 0 ? 1 : 2;
}

旋转卡壳

求法

凸包直径：凸包边上相距最远的两个点

然而可以证这两个点一定是顶点

于是用两条平行线去截凸包，这两条线一直在旋转

枚举边，则对面的点要到这条线的距离最远，这个点叫对踵点

可以用三角形的面积大小判断距离，这样用指针扫这个点，总的移动量是 O(n)

注意判断凸包的边也可能是答案，且面积相等要继续走，但是特判成一条直线的情况，不然会死循环

P1452 [USACO03FALL] Beauty Contest G /【模板】旋转卡壳

ll S(pll x, pll y, pll z) // x,y,z 形成三角形面积的 2倍 
{
    return llabs(cross(y - x, z - x));
}
ll cmp(const pll &x, const pll &y)
{
    double ang1 = atan2(x.se - p[1].se, x.fi - p[1].fi);
    double ang2 = atan2(y.se - p[1].se, y.fi - p[1].fi);
    if(!equ(ang1, ang2))    return ang1 < ang2;
    return len(x - p[1]) < len(y - p[1]);
}
void gettb()
{
    for(int i = 2; i <= n; ++i)
        if(p[i].se < p[1].se || (p[i].se == p[1].se && p[i].fi < p[1].fi))  swap(p[1], p[i]);
    sort(p + 2, p + n + 1, cmp), stk[++top] = 1; 
    for(int i = 2; i <= n; ++i)
    {
        while(top > 1 && cross(p[i] - p[stk[top]], p[stk[top]] - p[stk[top - 1]]) >= 0) --top;
        stk[++top] = i;
    }
    stk[++top] = 1;
    for(int i = 1; i <= top; ++i)   a[i] = p[stk[i]];
}
void xzqk() // 旋转卡壳 
{
    ll num = 3;
    if(top == 3)
    {
        ans = len(p[stk[top - 1]] - p[stk[top - 2]]);
        return;
    }    
    for(ll i = 1; i < top; ++i) // 枚举边 
    {
        ans = max(ans, len(a[i + 1] - a[i])); // 用凸包上的边更新 
        while(S(a[i], a[i + 1], a[num]) <= S(a[i], a[i + 1], a[num % (top - 1) + 1]))   num = num % (top - 1) + 1;
        ans = max({ans, len(a[num] - a[i]), len(a[num] - a[i + 1])}); // 距离更大就继续走，更新答案 
    }
}

应用

• 平面最远点对

• 最小矩形覆盖 P3187 [HNOI2007] 最小矩形覆盖

可证这个矩形一定至少一条边与凸包的一条边重合

类似旋转卡壳的思路，找到对踵点，则矩形的上边卡住这个点，下边与枚举的边重合

同时还要维护矩形的宽，用点积判断边的端点与那个点形成的向量在下边上的投影长度，也是用指针扫 l,r

先移动 l，注意初始时 r 应在对踵点 num 的位置，因为 l,r 分别在对踵点两侧，r 如果在 l 那侧则会不满足单调性

求出下边，先用投影求出下面两个点，然后利用旋转求出上面两个点

pdd rotate(pdd x, ld ang)
{
    return mp(-x.se, x.fi); 
}
ld S(pdd x, pdd y, pdd z)   {return fabs((x - y) ^ (z - y));}
int cmp(const pdd &x, const pdd &y)
{
    ld ang1 = atan2l(x.se - p[1].se, x.fi - p[1].fi);
    ld ang2 = atan2l(y.se - p[1].se, y.fi - p[1].fi);
    if(!equ(ang1, ang2))    return ang1 < ang2;
    return len(x - p[1]) < len(y - p[1]);
}
void graham()
{
    for(int i = 2; i <= n; ++i)
        if(p[i].se < p[1].se || (equ(p[i].se, p[1].se) && p[i].fi < p[1].fi))   swap(p[i], p[1]);
    sort(p + 2, p + n + 1, cmp), stk[++top] = 1;
    for(int i = 2; i <= n; ++i)
    {
        while(top > 1 && ((p[i] - p[stk[top]]) ^ (p[stk[top]] - p[stk[top - 1]])) >= 0) --top;
        stk[++top] = i;
    }
    stk[++top] = 1;
    for(int i = 1; i <= top; ++i)   lsh[i] = p[stk[i]];
    memcpy(p, lsh, sizeof(p)), --top;
}
void getans(pdd a, pdd b, pdd c, pdd d, pdd x)
{
    ld h = S(a, b, x) / len(b - a), lin = len(b - a); // 计算矩形的宽 
    ld L = fabsl((c - a) * (b - a)) / lin, R = fabsl((d - b) * (a - b)) / lin; // 投影长度 
    ld sum = (L + R - lin) * h; // 计算面积 
    if(sum >= mj)   return;
    mj = sum;
    ans[0] = b + R / lin * (a - b); // 顶点
    ans[1] = a + L / lin * (b - a);
    ans[2] = ans[1] + h / lin * rotate(b - a, pi / 2); // 加上 b-a 逆时针旋转 90°得到的单位向量 * h 
    ans[3] = ans[0] + h / lin * rotate(b - a, pi / 2);
}
void work()
{
    if(top == 2)
    {
        mj = 0, ans[0] = ans[2] = p[1], ans[1] = ans[3] = p[2];
        return;
    }
    num = 3, l = r = 3;
    for(int i = 1; i <= top; ++i)
    {
        int nxt = num % top + 1;
        while(S(p[num], p[i], p[i + 1]) <= S(p[nxt], p[i], p[i + 1]))   num = nxt, nxt = num % top + 1;
        nxt = l % top + 1;
        while((p[l] - p[i]) * (p[i + 1] - p[i]) <= (p[nxt] - p[i]) * (p[i + 1] - p[i])) 
            l = nxt, nxt = l % top + 1;  // 不加绝对值！正负的判断是需要的！必须转到投影为正 
        if(i == 1)   r = num; // 注意 r初始要在另一边，否则单调性不满足 
        nxt = r % top + 1;
        while((p[r] - p[i + 1]) * (p[i] - p[i + 1]) <= (p[nxt] - p[i + 1]) * (p[i] - p[i + 1]))
            r = nxt, nxt = r % top + 1;
        getans(p[i], p[i + 1], p[l], p[r], p[num]);
    }
}
int main()
{
    cin >> n;
    for(int i = 1; i <= n; ++i) cin >> p[i].fi >> p[i].se;
    graham();
    work(), id = 0;
    for(int i = 1; i < 4; ++i)
        if(ans[i].se < ans[id].se || (equ(ans[i].se, ans[id].se) && ans[i].fi < ans[id].fi))    id = i;
    printf("%.5Lf\n", mj); 
    for(int i = id; i < 4; ++i) printf("%.5Lf %.5Lf\n", sgn(ans[i].fi), sgn(ans[i].se));
    for(int i = 0; i < id; ++i) printf("%.5Lf %.5Lf\n", sgn(ans[i].fi), sgn(ans[i].se));
    return 0;
}

半平面交

半平面：一条直线和直线的一侧

半平面是一个一条直线和直线的一侧构成的点集，当包含直线时，称为闭半平面；当不包含直线时，称为开半平面

在计算几何中用向量表示，整个题统一以向量的左侧或右侧为半平面

半平面交：多个半平面的点集的交集

求法

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

规定向量左侧为目标半平面

先将划分的向量按极角逆时针排序，如果极角相同，在左侧的优先，用叉积判断

因为在左侧的半平面被右侧的完全包含，因此右侧的不需要

维护队列表示在边界上的向量，它们的交点应该形成了一个凸包

先弹出队尾的，如果原来的交点在新向量的右侧，则队尾不需要了

队头同理

最后再用队头弹出队尾

如果队列中的向量少于 3 个，则无解，因为题目中的面积不可能无限大

求出交点后计算凸包面积即可

struct line
{
    pdd st, ed; double ang;
    void getang()   {ang = atan2(ed.se - st.se, ed.fi - st.fi);}
    bool operator<(const line &x)const
    {
        if(!equ(ang, x.ang))    return ang < x.ang;
        return ((ed - x.st) ^ (x.ed - x.st)) < 0; // 极角相同，在左边的排前面 
    }
}p[N];
pdd cross(line x, line y) // 求直线（向量）x,y的交点 
{
    pdd w = y.ed - y.st, v = x.ed - x.st;
    return x.st + ((w ^ (x.st - y.st)) / (v ^ w)) * v;
}
int onright(line x, pdd y) // y 在 x右边，返回 1 
{
    return ((y - x.st) ^ (x.ed - x.st)) >= 0;
}
void halfplane()
{
    sort(p + 1, p + cnt + 1);
    for(int i = 1; i <= cnt; ++i) // 注意，不能把最开头 2个点直接加入，否则它们可能平行！ 
    {
        if(i > 1 && equ(p[i - 1].ang, p[i].ang))    continue;
        while(hh < tt && onright(p[i], cross(p[q[tt]], p[q[tt - 1]])))  --tt; // 注意先弹出队尾 
        while(hh < tt && onright(p[i], cross(p[q[hh]], p[q[hh + 1]])))  ++hh;
        q[++tt] = i;
    }
    while(hh < tt && onright(p[q[hh]], cross(p[q[tt]], p[q[tt - 1]])))  --tt; // 队头还可以更新队尾 
    if(tt - hh < 2) return; // 至少有 3个点 
    for(int i = hh; i < tt; ++i)    a[++idx] = cross(p[q[i]], p[q[i + 1]]);
    a[++idx] = cross(p[q[tt]], p[q[hh]]);
    for(int i = 2; i < idx; ++i)    ans += S(a[1], a[i], a[i + 1]); // 求这个多边形的面积 
}
int main()
{
    cin >> n;
    for(int i = 1; i <= n; ++i)
    {
        cin >> m[i];
        cin >> xi >> yi, stx = xi, sty = yi;
        for(int j = 1; j < m[i]; ++j)
        {
            cin >> xj >> yj;
            p[++cnt].st = mp(xi, yi), p[cnt].ed = mp(xj, yj), p[cnt].getang();
            xi = xj, yi = yj;
        }
        p[++cnt].st = mp(xi, yi), p[cnt].ed = mp(stx, sty), p[cnt].getang();
    }   
    halfplane();
    printf("%.3f", ans / 2);
    return 0;
}

应用

P3256 [JLOI2013] 赛车

画出 x-t 图像，发现符合要求的车代表的向量在半平面交的边界上

注意半平面强制规定在第一象限，要加入 y 轴以限定范围

此时出现某个向量刚好过前两个的交点时，就不能弹出，因为并列第一也算

注意去重，以及巨卡精度，eps=10^{-18}……和 0 没有本质区别

好像用单调栈做没有这些问题

struct line
{
    pdd st, ed; ld ang; ll id, sd;
    void getang()   {ang = atan2l(ed.se - st.se, ed.fi - st.fi);}
    bool operator<(const line &x)const
    {
        if(sd != x.sd)  return sd < x.sd;
        return ((ed - x.st) ^ (x.ed - x.st)) < -eps; 
    }
}p[N];
pdd cross(line x, line y)
{
    pdd v = x.ed - x.st, w = y.ed - y.st; 
    return x.st + ((w ^ (x.st - y.st)) / (v ^ w)) * v;
}
ll onright(line x, pdd y)   
{
    ld res = (y - x.st) ^ (x.ed - x.st);
    return res > eps; // 可能有的线只出现在交点处，为了判断这里就不删除 
}
void halfplane()
{
    sort(p + 1, p + cnt + 1);
    for(ll i = 1; i <= cnt; ++i)
    {
        if(i > 1 && equ(p[i].ang, p[i - 1].ang))    continue;
        while(hh < tt && onright(p[i], cross(p[q[tt]], p[q[tt - 1]])))  
            --tt;
        while(hh < tt && onright(p[i], cross(p[q[hh]], p[q[hh + 1]])))  
            ++hh;
        q[++tt] = i;
    }
    while(hh < tt && onright(p[q[hh]], cross(p[q[tt]], p[q[tt - 1]])))  
        --tt;
    for(ll i = hh; i <= tt; ++i)    book[p[q[i]].id] = 1;
}
int main()
{
    n = read();
    for(ll i = 1; i <= n; ++i)  ps[i] = read(), fa[i] = i;
    for(ll i = 1; i <= n; ++i)  k[i] = read(), qc[mp(ps[i], k[i])].pb(i);
    p[++cnt].st = mp(0, 0), p[cnt].ed = mp((ld)1, 0), p[cnt].getang(), p[cnt].id = 0, p[cnt].sd = 0;
    p[++cnt].st = mp(0, (ld)1), p[cnt].ed = mp(0, 0), p[cnt].getang(), p[cnt].id = 0, p[cnt].sd = -inf; 
    for(auto it = qc.begin(); it != qc.end(); ++it)
    {
        vector<ll> nn = it -> se;   ll i = nn[0];
        p[++cnt].st = mp(0, (ld)ps[i]), p[cnt].ed = mp(1, (ld)ps[i] + (ld)k[i]), p[cnt].getang(), p[cnt].id = i, p[cnt].sd = k[i];
        for(int j : nn) rel[j] = i;
    }
    halfplane();
    for(ll i = 1; i <= n; ++i)  ans += book[rel[i]];
    print(ans), putchar('\n');
    for(ll i = 1; i <= n; ++i)
        if(book[rel[i]])    print(i), putchar(' ');
    return 0; 
}

闵可夫斯基和

定义

两个点集 A,B 的闵可夫斯基和是 C=\{a+b,a\in A,b\in B\}

或从原点向点集 A 的每一个点做向量，将整个点集 B 沿每个向量移动，所有点的最终位置的并便是闵可夫斯基和

性质

• 两个大小分别为 n,m 的凸包的闵可夫斯基和仍是凸包，边是由原凸包上的边构成的

• 两个大小分别为 n,m 的点集的闵可夫斯基和，只取凸包的话大小为 \le n+m

• 如果 f 为凸函数，且 f_x=\max_{i+j=x}\{g_i+h_j\}，则 f 可以看作 g,h 对应点集的闵可夫斯基和

  因为对凸包差分，则差分数组， 也就是斜率单调，那么 f_x 可看作 g,h 各取某个前缀相加且总长度为 x，让总和最大

  贪心的选总共的前 x 大就好了，卷积后也就是对差分数组归并排序

  写成这样后，把 wqs 二分的朴素方程列出来，j 看作自变量，DP 值为因变量，则写成 f_k=\max_{j+u=k}\{f'_j+f''_u\} 的形式

  然后把 DP 写成区间 DP 的形式，可以分治，每层用闵可夫斯基和优化，求出对于每个 m=0\sim k f_{n,m} 的值

求法

求出两个点集的凸包，这里用的 Graham，已经以最下左的点为起点，按极角排序了

然后把它们差分，差分即为向量相减，注意凸包是环

加入两个点集第一个点的和，对差分数组归并排序，极角小的排在前面，用叉积判断

做一遍前缀和就得到了闵可夫斯基和，最后可以再求一遍凸包，因为要排除多点共线的情况

vector<pdd> minkowski(vector<pdd> p, vector<pdd> q)
{
    cnt1.clear(), cnt2.clear(), tb.clear();
    if(q.empty())   return p;
    if(p.empty())   return q;
    for(int i = 1; i < p.size(); ++i)   cnt1.pb(p[i] - p[i - 1]);   if(cnt1.size() > 1) cnt1.pb(p[0] - p.back());
    for(int i = 1; i < q.size(); ++i)   cnt2.pb(q[i] - q[i - 1]);   if(cnt2.size() > 1) cnt2.pb(q[0] - q.back());
    tb.pb(p[0] + q[0]);
    for(int i = 0, j = 0; ; )
    {
        if(i >= cnt1.size() && j >= cnt2.size())    break;
        if(j >= cnt2.size() || (i < cnt1.size() && sgn(cnt1[i] ^ cnt2[j]) >= 0))   tb.pb(tb.back() + cnt1[i]), ++i;
        else if(j < cnt2.size())    tb.pb(tb.back() + cnt2[j]), ++j;
    }
    return tb;
}
a = graham(a), b = graham(b);
ans = minkowski(a, b), ans = graham(ans);

应用

UVA10256 The Great Divide

对红蓝点集求凸包，问题转化为判断凸包是否有交

可以用旋转卡壳做，但是也可以转化为闵可夫斯基和