AGC051E[Middle Point] 题解

题面链接#

简要题意#

给定若干整点，可以对任意两个点取中点并加入点集，问在有限次操作后点集中最多能够有多少个整点。

条件转化#

我们记：

$$M=\{x|x=\frac{a}{2^b},a,b\in \mathbb{Z}\}\\ M^*=\{x|x\in M,x\ge0\}$$

令下文向量均为二维向量，记给定点集为 ${\vec{p_n}}$

那么原题即为求满足 $\exists\{a_n\}$ 使得 $a_i\in M^*,\displaystyle{\sum a_i=1,\begin{pmatrix}u\\v\end{pmatrix}}=\sum_{i=1}^n \begin{pmatrix}x_i\\ y_i\end{pmatrix}a_i$ 的整数点 $(u,v)$ 的个数。

注意到这类似于多元Bezet定理的形式，但是斐蜀定理需要对系数进行辗转相减法，这意味着系数需要满足对于加法和减法的封闭性，而 $M^*$ 并不满足。

注意到满足条件的点必然在所有点构成的凸包上或凸包内，对于凸包上的点，我们可以快速计算，接下来考虑凸包内的点对答案的贡献。

我们不妨猜想，对于凸包内的点，我们可以把条件中的 $M^*$ 放宽到对加减法都封闭的 $M$。
对于点在凸包内，将其平移至原点，假如我们已经有数列 $\{a'_n\}$ 满足 $a'_i\in M,\displaystyle{\sum a'_i=1},\vec{0}=\sum_{i=1}^n \begin{pmatrix}x_i\\ y_i\end{pmatrix}a'_i$，我们试图构造一个数列 $\{c'_n\}$，满足 $c'_i\in M^*,\displaystyle{\sum c'_i=1},\vec{0}=\sum_{i=1}^n \begin{pmatrix}x_i\\ y_i\end{pmatrix}c'_i$，取足够大的 $\alpha,\beta\in \mathbb{N}$，则等价于

$$\begin{align*} \exists\{a_n\}\ s.t.a_i\in \mathbb{Z},\displaystyle{\sum a_i=2^\alpha},\vec{0}=\sum_{i=1}^n \begin{pmatrix}x_i\\ y_i\end{pmatrix}a_i\\\Rightarrow \exists\{c_n\}\ s.t.c_i\in\mathbb{N},\sum c_i=2^\beta,\vec{0}=\sum_{i=1}^n \begin{pmatrix}x_i\\ y_i\end{pmatrix}c_i \end{align*}$$

我们考虑找到一个序列 $\{b_n\}$ 满足 $b_i\in \mathbb{N}^*,\displaystyle{\sum b_i=X,\sum_{i=1}^n \begin{pmatrix}x_i\\ y_i\end{pmatrix}c_i=\vec{0}}$，将 $\{a_n\}$ 加上若干倍 $\{b_n\}$ 构造出 $\{c_n\}$。显然有 $(u,v)\in ConvexHull \Leftrightarrow \exists \{b_n\}$，令 $c_i=a_i+kb_i\ge0$，则有 $\displaystyle\sum_{i=1}^n c_i=(2^\alpha+kX)$，则有 $kX=2^\beta-2^\alpha$，令 $k=2^\alpha g$，则有 $Xg+1=2^{\beta-\alpha}$，而 $\beta$ 可以任意大，故 $g$ 和 $k$ 可以任意大，必然存在一个充分大的 $k$，使得 $c_i\ge 0$，得证，平移回去则得到凸包内任意点的证明。

关于凸包的证明，利用凸包的性质即得，这里不再赘述。

又注意到 $\displaystyle{\sum c_i=2^\beta}$ 十分棘手，不妨做个变形：

$$\begin{align*} &\sum \vec{p_i}c_i=2^\beta\vec{x}\\ \Leftrightarrow &\sum_{i\ne 1} \vec{p_i}c_i=2^\beta\vec{x}-\vec{p_1}(2^\beta-\sum_{i\ne1}c_i)\\ \Leftrightarrow &\sum_{i\ne 1} (\vec{p_i}-\vec{p_1})c_i=2^\beta(\vec{x}-\vec{p_1})\\ \Leftrightarrow &\sum_{i\ne 1} \vec{p'_i}c_i=2^\beta\vec{x'} \end{align*}$$

解法#

参考多元Bezet定理，我们最终将能够将方程写作 $a\vec{u}+b\vec{v}=\vec{x}$ 的形式，其中 $u_x,u_y,v_y\in M,v_x=0$，则有

$$\begin{align*} & x_x=au_x\\ & x_y=au_y+bv_y \end{align*}$$

记 $f(x)$ 表示 $x$ 的因子 $2$ 的个数，$g(x)$ 表示对于 $x\in M$，分母中 $2$ 的幂次。
记 $\alpha=f(u_x),\beta=f(u_y),\theta=f(v_y),u_x=2^\alpha m,u_y=2^\beta n,v_y=2^\theta l$

• $\alpha\le \beta$
  因 $x_x$ 为整数，故 $g(a)\le \alpha\le \beta$，故 $au_y\in \mathbb{Z}$，因此 $bv_y$ 必须为整数，则有 $g(b)\le \theta$，将 $x$ 轴拉长 $2^{\alpha}$，$y$ 轴拉长 $2^{\theta}$，计算将基向量换成 ${\vec{u},\vec{v}}$ 的凸包中的整点个数即可，利用毕克定理即可解决。
• $\alpha> \beta$
  考虑将其转化到第一种情况：
  • $\theta<\beta$
    考虑将 $u_y$ 加上 $kv_y$，使 $f(u_y+kv_y)\ge\alpha$，即 $2^\beta n+2^\theta km=2^\alpha t\Leftrightarrow 2^{\beta-\theta}n+km=2^{\alpha-\theta}t \Leftrightarrow -2^{\beta-\theta}n\equiv km(\operatorname{mod}\ 2^{\alpha-\theta})$，取 $k=-2^{\beta-\theta}nm^{-1}$ 即可。
  • $\theta\ge\beta$
    考虑对 $y$ 维做辗转相减法，最终会得到
    $$\begin{align*} &\vec{u'}=\begin{pmatrix}k_1u_x\\ \gcd(u_y,v_y)\end{pmatrix}\\ &\vec{v'}=\begin{pmatrix}k_2u_x\\ 0\end{pmatrix} \end{align*}$$
    有 $f(u'_y)=\beta,f(u'_x)\ge \alpha>\beta$，即 $f(u'_y)>f(u'_x)$，交换 $x$ 维和 $y$ 维后回到情况一。

注意事项#

注意到在处理多元Bezet定理时，如果使用欧几里得算法，在缩减一维的情况下容易导致另一维的值域指数爆炸，而欧几里得算法的值域相当玄学，我们考虑换一个不那么玄学的算法（thanks https://codeforces.com/blog/entry/86100 @ ecnerwala），我们每次将三个向量的其中一个变成 $\vec{0}$，在这过程中我们考虑沿着 $max{|\vec{v_i}|}$ 最小的贪心策略消减，或者因为欧几里得距离不好实现，我们也可以把欧几里得距离换成曼哈顿距离。

代码#

#include <bits/stdc++.h>
using namespace std;
using i128 = __int128;
using i64 = long long;
const int M = 110;
inline i128 aabs(i128 x) { return x < 0 ? -x : x; }
struct vec2 {
    i128 x, y;
    i128 norm() { return aabs(x) + aabs(y); }
} a[M], v[3], vbef[3], vnow[3], _;
inline vec2 operator*(vec2 &a, i128 b) { return {a.x * b, a.y * b}; }
inline vec2 operator+(vec2 a, vec2 b) { return {a.x + b.x, a.y + b.y}; }
int n;
i64 ans = 0, __, flag;
mt19937 rnd(time(0));
vector<i64> met;
i128 x, y, t;
inline vec2 operator-(vec2 &a, vec2 &b) { return {a.x - b.x, a.y - b.y}; }
inline i128 operator*(vec2 a, vec2 b) { return a.x * b.y - a.y * b.x; }
vector<vec2> convex_hull(vector<vec2> a) {
    sort(a.begin(), a.end(), [](vec2 a, vec2 b) { return a.x == b.x ? a.y < b.y : a.x < b.x; });
    vec2 stk[M];
    vector<vec2> res;
    int t = -1;
    for (vec2 v : a) {
        while (t >= 1 && (stk[t] - stk[t - 1]) * (v - stk[t - 1]) < 0) t--;
        stk[++t] = v;
    }
    for (int i = 0; i <= t; i++) res.push_back(stk[i]);
    reverse(a.begin(), a.end()), t = -1;
    for (vec2 v : a) {
        while (t >= 1 && (stk[t] - stk[t - 1]) * (v - stk[t - 1]) < 0) t--;
        stk[++t] = v;
    }
    for (int i = 1; i <= t; i++) res.push_back(stk[i]);
    return res;
}
inline int ctz(i64 x) { return x = abs(x), x ? __builtin_ctzll(x) : 1010000000; }
void reduce_by_euclidean(vec2 &a, vec2 &b) {
    if (a.x < b.x) swap(a, b);
    while (a.x && b.x) t = a.x / b.x, a.x -= b.x * t, a.y -= b.y * t, swap(a, b);
}
struct update {
    int i, j;
    i64 k;
};
inline vec2 calupd(update u) { return v[u.i] + v[u.j] * u.k; }
inline void choose_update(update &bef, update now) {
    if (!now.k) return;
    vec2 vbef = calupd(bef), vnow = calupd(now);
    i128 sum = v[0].norm() + v[1].norm() + v[2].norm(), sbef = sum - v[bef.i].norm() + vbef.norm(), snow = sum - v[now.i].norm() + vnow.norm();
    if (sbef > snow) bef = now;
    else if (sbef == snow && (rnd() % 2)) bef = now;
}
inline int hh(int x) { return x == 1010000000 ? 0 : x; }
void reduce_by_greedy(vec2 &a, vec2 &b, vec2 &c) { // span(a, b, c) = span(a', b')
    v[0] = a, v[1] = b, v[2] = c;
    while ((v[0].x || v[0].y) && (v[1].x || v[1].y) && (v[2].x || v[2].y)) {
        update upd = {0, 0, 0};
        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
                if (i != j) {
                    met.clear();
                    if (v[j].x) __ = -v[i].x / v[j].x, met.push_back(__), met.push_back(__ - 1), met.push_back(__ + 1);
                    if (v[j].y) __ = -v[i].y / v[j].y, met.push_back(__), met.push_back(__ - 1), met.push_back(__ + 1);
                    for (int o : met) choose_update(upd, {i, j, o});
                }
        v[upd.i] = calupd(upd);
    }
    for (int o = 0; o < 2; o++)
        if (!v[o].x && !v[o].y) swap(v[o], v[2]);
    a = v[0], b = v[1], c = v[2];
}
void ex__gcd(i128 a, i128 b, i128 &x, i128 &y) {
    if (!b) return x = 1, y = 0, void();
    ex__gcd(b, a % b, x, y);
    t = y, y = x - (a / b) * y, x = t;
}
i128 inv(i128 a, i128 b) { return a %= b, ex__gcd(b, a, x, y), y; }
int main() {
    scanf("%d", &n);
    for (int i = 0, x, y; i < n; i++) scanf("%d%d", &x, &y), a[i].x = x, a[i].y = y; // ast::p[i] = {a[i].x, a[i].y};
    for (int i = 1; i < n; i++) a[i].x -= a[0].x, a[i].y -= a[0].y;
    a[0] = {0, 0};
    vector<vec2> t = convex_hull(vector<vec2>(a, a + n)), s;
    for (int i = 1; i < t.size(); i++) _ = t[i] - t[i - 1], ans += 1 << min(ctz(_.x), ctz(_.y));
    for (int i = 3; i < n; i++) reduce_by_greedy(a[1], a[2], a[i]);
    vec2 u = a[1], v = a[2], _;
    reduce_by_euclidean(u, v);
    if (!u.x || !u.y) swap(u, v);
    if (!v.y) swap(u.x, u.y), swap(v.x, v.y), flag = 1;
    v.y = aabs(v.y);
    int alpha = ctz(u.x), beta = ctz(u.y), theta = ctz(v.y);
    i128 a = u.x >> hh(alpha), b = u.y >> hh(beta), c = v.y >> hh(theta);
    if (alpha > beta)
        if (theta > beta) swap(u.x, u.y), swap(v.x, v.y), reduce_by_euclidean(u, v), swap(u.x, u.y), swap(v.x, v.y);
        else u.y = (b * inv((i128) 1 << (alpha - beta), c)) << alpha;
    int delta = min(ctz(u.x), ctz(u.y)), epsilon = min(ctz(v.x), ctz(v.y));
    u.x >>= hh(delta), u.y >>= hh(delta), v.x >>= hh(epsilon), v.y >>= hh(epsilon);
    if (flag) swap(u.x, u.y), swap(v.x, v.y);
    for (vec2 g : t)
        s.push_back({(g.x * v.y - g.y * v.x) / (u.x * v.y - u.y * v.x), (g.x * u.y - g.y * u.x) / (v.x * u.y - v.y * u.x)}); // g.pr(), printf("\n");
    i128 area = 0, border = 0, in;
    for (int i = 1; i < s.size(); i++) _ = s[i] - s[i - 1], area += s[i - 1] * s[i], border += aabs(__gcd(_.x, _.y));
    area = aabs(area), in = (area - border + 2) / 2;
    printf("%lld\n", (i64) (ans + in));
    return 0;
}