AGC051E[Middle Point] 题解
题面链接
简要题意
给定若干整点，可以对任意两个点取中点并加入点集，问在有限次操作后点集中最多能够有多少个整点。
条件转化
我们记：
令下文向量均为二维向量，记给定点集为 \({\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}\)，则等价于
我们考虑找到一个序列 \(\{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}\) 十分棘手，不妨做个变形：
解法
参考多元Bezet定理，我们最终将能够将方程写作 \(a\vec{u}+b\vec{v}=\vec{x}\) 的形式，其中 \(u_x,u_y,v_y\in M,v_x=0\)，则有
记 \(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\) 维后回到情况一。
- \(\theta<\beta\)
注意事项
注意到在处理多元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;
}
