ABC202F - Integer Convex Hull 题解

本文同步发表在 hexo 博客

Description

平面内 $N$ 个坐标均为整数的点，其中任意三点不共线，令所有点 $P_i$ 构成一个全集 $U$，子集 $S\subseteq U$ 且 $|S|\ge 3$。问凸包面积为整数的子集 $S$ 的个数模 $10^9+7$ 的值。

$3\le N\le 80$。

Solution

比较神奇的题目。

考虑 Andrew 算法的执行流程：将所有点按照先 $x$ 轴升序后 $y$ 轴升序排序，然后从左往右弄出下凸包，再从右往左弄出上凸包。本题可以利用类似的思路。

然后我们由叉积计算面积的方法知道，整点三角形的面积乘二之后必然是整数。于是我们考虑依次用三角形弄出凸包，具体地，可以 dp：

• $O(n)$ 枚举凸包最左端的端点 $P_s$。

• ${\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}_{i,jk}$：凸包满足如下条件的点集 $S$ 的数量：

  • 从 $P_s$ 顺时针考虑每个点，最后两个分别是 $P_i$ 和 $P_j$；
  • 所有的点都在直线 $P_s{P}_j$ 以上；
  • 面积乘以二之后模二为 $k$。

• ${\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}}_{i,j,k}$：凸包满足如下条件的点集 $S$ 的数量：

  • 从 $P_s$ 逆时针考虑每个点，最后两个分别是 $P_i$ 和 $P_j$；
  • 所有的点都在直线 $P_s{P}_j$ 以下；
  • 面积乘以二之后模二为 $k$。

• 确定每个最右端的点 $P_j$，然后枚举 $i$ 和 $k$ 统计答案：
  a n s s = ∑ j = s + 1 n ∑ k ∈ { 0 , 1 } ( ∑ i u p p e r i , j , k × ∑ i l o w e r i , j , k ) \mathrm{ans}_s = \sum_{j = s + 1}^n\sum_{k \in \{0, 1\}}\left(\sum_{i}\mathrm{upper}_{i,j,k} \times \sum_{i}\mathrm{lower_{i,j,k}} \right) anss​=j=s+1∑n​k∈{0,1}∑​(i∑​upperi,j,k​×i∑​loweri,j,k​)

这样做的原理其实就是合并上下凸壳：根据上面设计的 dp 状态可知我们确定了左端点 $P_s$ 和右端点 $P_j$ 后，满足上凸壳的点集个数即为枚举倒数第二个点 $i$ 然后求和的结果 ${\sum }_i{\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}_{i,j,k}$，满足下凸壳的点集个数即为 ${\sum }_i{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}}_{i,j,k}$。

怎么计算 $\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}$ 和 $\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}$ 呢？我们可以通过 ${\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}_{i,j,k}$ 递推出 ${\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}_{j,l,k}$，其原理就是判断 $P_l$ 能否被加入凸壳中，能则更新，不能就更新 ${\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}}_{j,l,k}$（加入下凸壳）。新加入的部分即为 $△P_s{P}_j{P}_l$，需要根据这个面积的二倍的奇偶性判断更新哪个 $k$。

这个做法要求我们提前预处理每个三角形的面积的二倍的奇偶性 ${\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{y}}_{i,j,k}$ 和三角形内点的个数 ${\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{d}\mathrm{e}}_{i,j,k}$。总时间复杂度 $O(n^4)$。

#include <cstdio>
#include <cctype>
#include <cmath>
#include <algorithm>
#define il inline
#define FOR(i, a, b) for (int i = a; i <= b; ++i)
#define DEC(i, a, b) for (int i = a; i >= b; --i)

const int maxn = 85;

int read()
{
    int s = 0, x = 0;
    char c = getchar();
    while (!isdigit(c))
        x |= (c == '-'), c = getchar();
    while (isdigit(c))
        s = s * 10 + c - '0', c = getchar();
    return x ? -s : s;
}

struct Point
{
    int x, y;
    Point() {}
    Point(int _x, int _y) {x = _x, y = _y;}
    il bool operator<(const Point &b) const {return x == b.x ? y < b.y : x < b.x;}
} P[maxn];

typedef Point Vector;

struct modint
{
    typedef long long ll;
    const static ll mod = 1e9 + 7;
    ll val;
    modint(int _val = 0) {val = _val;}
    modint operator+(const modint &b) const {return modint((val + b.val) % mod);}
    modint operator-(const modint &b) const {return modint((val - b.val + mod) % mod);}
    modint operator*(const modint &b) const {return modint((val * b.val) % mod);}
    modint operator+=(const modint &b) {return *this = *this + b;}
    modint operator*=(const modint &b) {return *this = *this * b;}
    modint operator-=(const modint &b) {return *this = *this - b;}
    modint operator+(const int &b) const {return modint((val + b) % mod);}
    modint operator-(const int &b) const {return modint((val - b + mod) % mod);}
    modint operator*(const int &b) const {return modint((val * b) % mod);}
    modint operator+=(const int &b) {return *this = *this + b;}
    modint operator-=(const int &b) {return *this = *this - b;}
    modint operator*=(const int &b) {return *this = *this * b;}
};

il Vector operator-(const Point &a, const Point &b) {return Point(b.x - a.x, b.y - a.y);}
il int operator^(const Vector &a, const Vector &b) {return a.x * b.y - a.y * b.x;}//叉积

int area2(Point x, Point y, Point z)//返回面积的二倍
{
    return abs((y - x) ^ (z - x));
}

int n, parity[maxn][maxn][maxn], inside[maxn][maxn][maxn];
modint pow2[maxn], upper[maxn][maxn][2], lower[maxn][maxn][2];

int main()
{
    n = read();
    FOR(i, 1, n) P[i].x = read(), P[i].y = read();
    std::sort(P + 1, P + n + 1);//对点进行排序

    FOR(i, 1, n)
        FOR(j, 1, n)
            FOR(k, 1, n)
            {
                if (i == j || j == k || i == k) continue;
                parity[i][j][k] = area2(P[i], P[j], P[k]) % 2;//更新面积奇偶性
                FOR(l, 1, n)//枚举每个 l 点判断 l 在不在三角形内
                    if (l != i && l != j && l != k && area2(P[i], P[j], P[k]) == area2(P[l], P[j], P[k]) + area2(P[i], P[l], P[k]) + area2(P[i], P[j], P[l]))
                        inside[i][j][k] += 1;
            }
    

    pow2[0] = 1;
    FOR(i, 1, n) pow2[i] = pow2[i - 1] * 2;//预处理二的次幂

    modint ans = 0;

    DEC(must, n, 1)//枚举最左边的端点
    {
        FOR(i, must, n)
            FOR(j, must, n)
                FOR(k, 0, 1)
                    upper[i][j][k] = lower[i][j][k] = 0;//清空数组
        FOR(i, must + 1, n) upper[must][i][0] = lower[must][i][0] = 1;//边界条件
        FOR(i, must, n)
            FOR(j, must + 1, n)
                FOR(k, 0, 1)
                    FOR(l, j + 1, n)
                        if (((P[l] - P[j]) ^ (P[j] - P[i])) > 0)//加入上凸壳
                            upper[j][l][k ^ parity[must][j][l]] += upper[i][j][k] * pow2[inside[must][j][l]];//乘法原理
                        else lower[j][l][k ^ parity[must][j][l]] += lower[i][j][k] * pow2[inside[must][j][l]];
        FOR(j, must + 1, n)//枚举最右端点
            FOR(k, 0, 1)
            {
                modint up = 0, lo = 0;//上下凸壳
                FOR(i, must, j - 1)
                    up += upper[i][j][k], lo += lower[i][j][k];
                ans += up * lo;
            }
    }

    printf("%d\n", (ans - n * (n - 1) / 2).val);//最后要去除线段的贡献
    return 0;
}