Submission #35203333

---

Source Code Expand

// O(N^3 + MAX^2) 時間 / 全部整数 (MAX^7 * 3^7 の精度が必要)
// O(N^3) 凸包 (2 点とって辺を凸包に採用できるか)
// 半空間の積に分解
// z 座標で切って半平面の積
// y 座標で切ってしゃくとり
// 2D は gcd をうまく使うと単位格子の面積を求められる
#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;
using ll = long long;
using Modint = atcoder::modint998244353;
bool chmin(ll& a, ll b){ if(a <= b) return 0; a = b; return 1; }
bool chmax(ll& a, ll b){ if(a >= b) return 0; a = b; return 1; }
template<class T, class F = equal_to<>> void uniq(T& A, F f = equal_to<>{}){
    A.erase(unique(A.begin(), A.end(), f), A.end());
}
template<class T, class F> void remove_if(T& A, F f){
    A.erase(remove_if(A.begin(), A.end(), f), A.end());
}
ll div_floor(ll a, ll b){
    assert(b > 0);
    ll ans = a / b;
    if(a % b < 0) ans--;
    return ans;
}
ll div_ceil(ll a, ll b){
    assert(b > 0);
    ll ans = a / b;
    if(a % b > 0) ans++;
    return ans;
}
ll gcd(ll a, ll b, ll c){
    return gcd(a, gcd(b, c));
}


struct Point{
    ll x = 0, y = 0, z = 0;
    Point& operator-=(Point a){
        x -= a.x; y -= a.y; z -= a.z;
        return *this;
    }
    Point& operator/=(ll a){
        x /= a; y /= a; z /= a;
        return *this;
    }
    Point operator-() const {
        return {-x, -y, -z};
    }
    friend bool operator==(Point a, Point b){
        return a.x == b.x && a.y == b.y && a.z == b.z;
    }
    friend bool operator!=(Point a, Point b){
        return a.x != b.x || a.y != b.y || a.z != b.z;
    }
    friend bool operator<(Point a, Point b){
        return tie(a.x, a.y, a.z) < tie(b.x, b.y, b.z);
    }
    friend bool operator>(Point a, Point b){
        return tie(a.x, a.y, a.z) > tie(b.x, b.y, b.z);
    }
    friend Point operator+(Point a, Point b){
        return {a.x + b.x, a.y + b.y, a.z + b.z};
    }
    friend Point operator-(Point a, Point b){
        return {a.x - b.x, a.y - b.y, a.z - b.z};
    }
    friend Point operator*(ll a, Point b){
        return {b.x * a, b.y * a, b.z * a};
    }
    friend ll dot(Point a, Point b){
        return a.x * b.x + a.y * b.y + a.z * b.z;
    }
    friend Point cross(Point a, Point b){
        return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
    }
};
bool one_line(Point a, Point b){
    return cross(a, b) == Point{};
}
bool in_segment(Point a, Point b, Point x){
    a -= x;
    b -= x;
    return one_line(a, b) && dot(a, b) <= 0;
}

Modint solve1D(Point a, Modint K){
    // 全ての点が 1 直線上にある場合
    return K * gcd(a.x, a.y, a.z) + 1;
}
Modint solve2D(vector<Point> poly, Modint K){
    // 全ての点が 1 平面上にある場合
    // 辺上の格子点の数、面積 / 単位格子の面積 を求めてエルハート多項式 (ピックの定理)
    const ll N = poly.size();
    vector<Modint> ehrhart = {1, 0, 0};
    // 面積 / 単位格子の面積
    for(ll i = 1; i + 1 < N; i++){
        Modint area = 1;
        Point u = poly[i] - poly[0], v = poly[i + 1] - poly[0];
        {
            const ll g = gcd(u.x, u.y, u.z);
            u /= g;
            area *= g;
        }
        {
            const ll g = gcd(v.x, v.y, v.z);
            v /= g;
            area *= g;
        }
        u = cross(u, v);
        area *= gcd(u.x, u.y, u.z);
        ehrhart[2] += area;
    }
    // 辺上の格子点の数
    poly.push_back(poly[0]);
    for(ll i = 0; i < N; i++){
        const Point v = poly[i + 1] - poly[i];
        ehrhart[1] += gcd(v.x, v.y, v.z);
    }
    ehrhart[2] /= 2;
    ehrhart[1] /= 2;
    return ((ehrhart[2] * K) + ehrhart[1]) * K + ehrhart[0];
}
Modint solve(vector<Point> A, Modint K){
    // 同一直線上に 3 点あれば内側の点を削除
    for(ll i = A.size(); i--; ) if([&]() -> bool {
        for(ll j = A.size(); j--; ) if(j != i){
            for(ll k = j; k--; ) if(k != i){
                // A[j] - A[i] - A[k] の順に並んでいるか？
                if(in_segment(A[j], A[k], A[i])) return 1;
            }
        }
        return 0;
    }()) A.erase(A.begin() + i);

    // 全ての点が 1 直線上にある場合
    if(A.size() == 2) return solve1D(A[1] - A[0], K);

    // 凸包の辺を列挙
    vector<pair<ll, ll>> edge;
    vector g(A.size(), vector<ll>{});
    for(ll i = A.size(); i--; ) for(ll j = i; j--; ) if([&]() -> bool {
        // 軸 A[i] - A[j] 回りの 180 度未満の領域に全ての点があるか？
        if(A.size() <= 3) return 1;
        vector<Point> B = A;
        const Point U = B[i];
        B.erase(B.begin() + i);
        for(Point& x : B) x -= U;
        const Point V = B[j];
        B.erase(B.begin() + j);
        // 軸成分 (V) を取り除いたときの cross == cross の軸成分 なので、軸成分を取り除かなくても良い
        // 偏角ソート
        B = [&]{
            // 2 分割してループを切り開いてソート、同じ向きは unique
            // 2 分割の方法 : 基準のベクトルを持ってきて cross が正か / 0 の場合は V と cross を取って 0 と比較
            vector<Point> upper, lower;
            const Point W = B[0];
            for(Point x : B){
                ll c = dot(V, cross(W, x));
                if(c == 0) (cross(V, x) < Point{} ? lower : upper).push_back(x);
                else (c < 0 ? lower : upper).push_back(x);
            }
            sort(lower.begin(), lower.end(), [&](Point a, Point b){ return dot(cross(a, b), V) > 0; });
            uniq(lower, [&](Point a, Point b){ return dot(cross(a, b), V) == 0; });
            sort(upper.begin(), upper.end(), [&](Point a, Point b){ return dot(cross(a, b), V) > 0; });
            uniq(upper, [&](Point a, Point b){ return dot(cross(a, b), V) == 0; });
            upper.insert(upper.end(), lower.begin(), lower.end());
            return upper;
        }();
        // 180 度より大きく開いている部分が 1 箇所あれば OK
        assert(B.size());
        if(B.size() == 1) return 1;
        B.push_back(B.front());
        return adjacent_find(B.begin(), B.end(), [&](Point a, Point b){ return dot(cross(a, b), V) < 0; }) != B.end();
    }()){
        edge.emplace_back(i, j);
        g[i].push_back(j);
        g[j].push_back(i);
    }

    // 辺 2 本で面を作る
    vector<pair<Point, ll>> planes;
    auto add_plane = [&](Point a, Point b, Point c){
        // 面 abc の方程式を求める
        // dot(cross(b - a, c - a), x - a) == 0
        // dot(cross(b - a, c - a), x) == dot(cross(b - a, c - a), a)
        Point abc = cross(b - a, c - a);
        // 正規化
        if(abc < Point{}) abc = -abc;
        const ll g = gcd(abc.x, abc.y, abc.z);
        abc /= g;
        planes.emplace_back(abc, dot(abc, a)); // {2 乗, 3 乗} オーダー
    };

    for(ll i = 0; i < A.size(); i++) if(g[i].size()){
        // ある頂点から出る 1 つの辺を軸として偏角ソート
        auto& adj = g[i];
        const ll j = adj[0];
        const Point V = A[j] - A[i];
        sort(adj.begin() + 1, adj.end(), [&](ll j, ll k){ return dot(cross(A[j] - A[i], A[k] - A[i]), V) > 0; });
        const ll n = adj.size();
        adj.push_back(j);
        // 隣り合う 2 本の辺で面を作る
        for(ll j = 0; j < n; j++) add_plane(A[i], A[adj[j]], A[adj[j + 1]]);
    }
    sort(planes.begin(), planes.end());
    uniq(planes);

    // 全ての点が 1 平面上にある場合
    if(planes.size() == 1){
        vector<Point> poly;
        ll u = edge[0].first, v = edge[0].second;
        while(poly.size() < edge.size()){
            poly.push_back(A[u]);
            tie(u, v) = pair{v, u ^ g[v][0] ^ g[v][1]};
        }
        return solve2D(poly, K);
    }

    // 面のどちら側にあるか (dot(V, x) ≤ C)
    for(auto& [V, C] : planes){
        for(Point x : A){
            ll d = dot(V, x);
            if(d == C) continue;
            if(d > C){
                V = -V;
                C = -C;
            }
            break;
        }
    }

    // z 座標で切って、半平面の積
    remove_if(planes, [](pair<Point, ll> pl){ auto [V, C] = pl; return V.x == 0 && V.y == 0; });
    // x, y のみで偏角ソート
    sort(planes.begin(), planes.end(), [](pair<Point, ll> a, pair<Point, ll> b){
        auto [V1, C1] = a;
        auto [V2, C2] = b;
        const bool u1 = pair{V1.x, V1.y} < pair{0LL, 0LL};
        const bool u2 = pair{V2.x, V2.y} < pair{0LL, 0LL};
        if(u1 != u2) return u1 < u2;
        return V1.x * V2.y - V1.y * V2.x > 0; // 4 乗
    });

    // 格子点を数える
    ll zmin = INT_MAX, zmax = INT_MIN;
    for(auto [x, y, z] : A){
        chmin(zmin, z);
        chmax(zmax, z);
    }
    auto count_lattice = [&]() -> Modint {
        Modint ans = 0;
        for(ll z = zmin; z <= zmax; z++){
            vector<array<ll, 3>> lines; // Ax + By <= C
            for(auto& [V, C] : planes){
                // 傾きを unique
                array<ll, 3> l = {V.x, V.y, C - V.z * z};
                if(lines.size()){
                    auto& [a, b, c] = lines.back();
                    if(a * l[1] == b * l[0]){
                        if(c * (abs(l[0]) + abs(l[1])) > l[2] * (abs(a) + abs(b))){ // 5 乗
                            lines.back() = l;
                        }
                        continue;
                    }
                }
                lines.push_back(l);
            }

            // いらない半平面を捨てる
            auto is_needed = [](array<ll, 3> l0, array<ll, 3> l, array<ll, 3> l1) -> bool {
                const __int128_t W = l0[0] * l1[1] - l0[1] * l1[0]; // W 倍して整数点に (4 乗オーダー)
                if(W <= 0) return 1; // 角 l0 - l1 が 180 度以上
                // l0, l1 の交点
                array<__int128_t, 2> p01 = {l0[2] * l1[1] - l0[1] * l1[2], l0[0] * l1[2] - l0[2] * l1[0]}; // 5 乗オーダー
                return l[0] * p01[0] + l[1] * p01[1] > l[2] * W; // 7 乗オーダー
            };
            vector<array<ll, 3>> poly;
            for(auto& l : lines){
                while(poly.size() >= 2 && !is_needed(poly.rbegin()[1], poly.back(), l)) poly.pop_back();
                poly.push_back(l);
            }
            {
                auto p = poly.begin(), q = prev(poly.end());
                while(q - p > 2){
                    if(!is_needed(*q, *p, *next(p))){
                        p++;
                        continue;
                    }
                    if(!is_needed(*prev(q), *q, *p)){
                        q--;
                        continue;
                    }
                    break;
                }
                poly = vector(p, next(q));
            }

            // y 座標で切る
            ll ymin = INT_MAX, ymax = INT_MIN;
            auto calc_y = [&](array<ll, 3> l0, array<ll, 3> l1){
                const ll p = l0[0] * l1[2] - l0[2] * l1[0], q = l0[0] * l1[1] - l0[1] * l1[0]; // 5 乗, 4 乗
                chmin(ymin, div_ceil(p, q));
                chmax(ymax, div_floor(p, q));
            };
            for(ll i = 0; i + 1 < poly.size(); i++) calc_y(poly[i], poly[i + 1]);
            calc_y(poly.back(), poly.front());
            assert(ymin <= ymax + 1);

            remove_if(poly, [](array<ll, 3> l){ return l[0] == 0; });
            auto mid = partition_point(poly.begin(), poly.end(), [](array<ll, 3> l){ return l[0] > 0; });
            vector<array<ll, 3>> lower(poly.rbegin(), reverse_iterator{mid}), upper(poly.begin(), mid);
            auto lo = lower.begin(), hi = upper.begin();
            for(ll y = ymin; y <= ymax; y++){
                // しゃくとりで x 座標の範囲を計算
                auto cross_y = [&](array<ll, 3> l0, array<ll, 3> l1){
                    const ll p = l0[0] * l1[2] - l0[2] * l1[0], q = l0[0] * l1[1] - l0[1] * l1[0]; // 5 乗, 4 乗
                    return div_ceil(p, q);
                };
                while(next(lo) != lower.end() && cross_y(*next(lo), *lo) <= y) lo++;
                while(next(hi) != upper.end() && cross_y(*hi, *next(hi)) <= y) hi++;
                auto calc_x_lo = [&](array<ll, 3> l){ return div_ceil(y * l[1] - l[2], -l[0]); };
                auto calc_x_hi = [&](array<ll, 3> l){ return div_floor(l[2] - y * l[1], l[0]); };
                ll xmin = calc_x_lo(*lo);
                ll xmax = calc_x_hi(*hi);
                assert(xmin <= xmax + 1);
                ans += xmax - xmin + 1;
            }
        }
        return ans;
    };

    // エルハート多項式を補間
    const Modint a = count_lattice() - 1;
    for(auto& [V, C] : planes) C *= 2;
    zmin *= 2;
    zmax *= 2;
    const Modint b = count_lattice() - 1;
    for(auto& [V, C] : planes){
        C /= 2;
        C *= 3;
    }
    zmin /= 2;
    zmax /= 2;
    zmin *= 3;
    zmax *= 3;
    const Modint c = count_lattice() - 1;
    vector<Modint> ehrhart = {1, (18 * a - 9 * b + 2 * c) / 6, (-15 * a + 12 * b - 3 * c) / 6, (3 * a - 3 * b + 1 * c) / 6};

    return ((ehrhart[3] * K + ehrhart[2]) * K + ehrhart[1]) * K + ehrhart[0];
}
int main(){
    ll N, K;
    cin >> N >> K;
    vector<Point> A(N);
    for(auto& [x, y, z] : A) cin >> x >> y >> z;
    cout << solve(move(A), K).val() << endl;
}

Submission Info

Compile Error

./Main.cpp: In function ‘Modint solve(std::vector<Point>, Modint)’:
./Main.cpp:195:21: warning: comparison of integer expressions of different signedness: ‘ll’ {aka ‘long long int’} and ‘std::vector<Point>::size_type’ {aka ‘long unsigned int’} [-Wsign-compare]
  195 |     for(ll i = 0; i < A.size(); i++) if(g[i].size()){
      |                   ~~^~~~~~~~~~
./Main.cpp: In lambda function:
./Main.cpp:306:33: warning: comparison of integer expressions of different signedness: ‘ll’ {aka ‘long long int’} and ‘std::vector<std::array<long long int, 3> >::size_type’ {aka ‘long unsigned int’} [-Wsign-compare]
  306 |             for(ll i = 0; i + 1 < poly.size(); i++) calc_y(poly[i], poly[i + 1]);
      |                           ~~~~~~^~~~~~~~~~~~~

Judge Result