数学っぽい問題好き！！！

問題概要

$N$ 個のカフェオレがある。 $i$ 番目のカフェオレは、コーヒーが $c_{i}$ mg、ミルクが $m_{i}$ mg だけ入っている (ともに整数)。

これらを混ぜて作れるカフェオレのうち、コーヒー量 $x$ とミルク量が $y$ がともに正の整数となるものが何個あるか、1000000007 で割ったあまりを求めよ。

制約

$1 \le N \le 10^{5}$
$1 \le c_{i}, m_{i} \le 10^{9}$

考えたこと

$N$ 個のベクトル $(c_{i}, m_{i}$ の線形結合で表せる多角形の内部および周上にある格子点の個数を数え上げる問題！！！

線形結合は、 $N$ 個のベクトルを偏角が小さい順にソート (最初 double 型で m / c の大小比較でやったら誤差死した) して、その順にベクトルを累積和していくことで求められる！！

詳しくはふるやんさんの公式解説より！！

www.creativ.xyz

ピックの定理

問題は、凸多角形の内部および周上に格子点が何個あるかを求める問題へと帰着された。

ここでピックの定理が使える！ピックの定理とは

多角形の面積を $S$
多角形の内部の格子点の個数を $I$
多角形の周上の格子点の個数を $B$

としたとき、

$S = I + \frac{B}{2} - 1$

が成立するというものだ。今回は、

$S$ は簡単に求められる
$B$ も各線分の格子点の個数を gcd など駆使して求められる

となっているので、逆にこれらを用いて $I$ が求められる。 $I$ が求められたら、 $I+B-1$ が答えとなる (原点 $(0, 0)$ を除くので 1 を引く)。

なお、floor sum による別解もあるので下の方で。

コード

計算量は $O(N \log N + N \log V)$ となる ( $V$ は登場する整数値の最大値)。

#include <bits/stdc++.h>
using namespace std;
using pll = pair<long long, long long>;

// modint
template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) val += MOD;
    }
    constexpr int getmod() const { return MOD; }
    constexpr Fp operator - () const noexcept {
        return val ? MOD - val : 0;
    }
    constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
    constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
    constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
    constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
    constexpr Fp& operator += (const Fp& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator == (const Fp& r) const noexcept {
        return this->val == r.val;
    }
    constexpr bool operator != (const Fp& r) const noexcept {
        return this->val != r.val;
    }
    friend constexpr istream& operator >> (istream& is, Fp<MOD>& x) noexcept {
        is >> x.val;
        x.val %= MOD;
        if (x.val < 0) x.val += MOD;
        return is;
    }
    friend constexpr ostream& operator << (ostream& os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD>& r, long long n) noexcept {
        if (n == 0) return 1;
        if (n < 0) return modpow(modinv(r), -n);
        auto t = modpow(r, n / 2);
        t = t * t;
        if (n & 1) t = t * r;
        return t;
    }
    friend constexpr Fp<MOD> modinv(const Fp<MOD>& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        return Fp<MOD>(u);
    }
};

const int MOD = 1000000007;
using mint = Fp<MOD>;

mint calc_area(vector<pll> vp) {
    int N = vp.size();
    mint res = 0;
    for (int i = 0; i < N; ++i) {
        mint px = vp[i].first, py = vp[(i+1)%N].first;
        mint qx = vp[i].second, qy = vp[(i+1)%N].second;
        res += px * qy - py * qx;
    }
    return res;
}

long long GCD(long long x, long long y) {
    if (y == 0) return x;
    else return GCD(y, x % y);
}

int main() {
    int N;
    cin >> N;
    vector<long long> X(N), Y(N);
    for (int i = 0; i < N; ++i) cin >> X[i] >> Y[i];
    vector<int> ids(N);
    iota(ids.begin(), ids.end(), 0);
    sort(ids.begin(), ids.end(), [&](int i, int j) {
        return Y[i]*X[j] < Y[j]*X[i];
    });
    vector<pll> vp;
    vp.emplace_back(0, 0);
    for (auto i : ids) {
        auto p = vp.back();
        vp.emplace_back(p.first + X[i], p.second + Y[i]);
    }
    auto p = vp.back();
    int M = vp.size();
    for (int i = 1; i < M-1; ++i) {
        vp.emplace_back(p.first - vp[i].first, p.second - vp[i].second);
    }

    mint S = calc_area(vp);
    mint B = 0;
    for (int i = 0; i < vp.size(); ++i) {
        long long dx = vp[(i+1)%vp.size()].first - vp[i].first;
        long long dy = vp[(i+1)%vp.size()].second - vp[i].second;
        dx = abs(dx), dy = abs(dy);
        long long g = GCD(dx, dy);
        B += g;
    }
    mint res = (S + B) / 2;
    cout << res << endl;
}

追記：floor sum による別解

floor sum でも通してみた。

#include <bits/stdc++.h>
using namespace std;

// sum_{i=0}^{n-1} floor((a * i + b) / m)
// O(log(n + m + a + b))
long long floor_sum(long long n, long long a, long long b, long long m) {
    if (n == 0) return 0;
    long long res = 0;
    if (a >= m) {
        res += n * (n - 1) * (a / m) / 2;
        a %= m;
    }
    if (b >= m) {
        res += n * (b / m);
        b %= m;
    }
    if (a == 0) return res;
    long long ymax = (a * n + b) / m, xmax = ymax * m - b;
    if (ymax == 0) return res;
    res += (n - (xmax + a - 1) / a) * ymax;
    res += floor_sum(ymax, m, (a - xmax % a) % a, a);
    return res;
}

// #lp under (and on) the segment (x1, y1)-(x2, y2)
// not including y = 0, x = x2
long long num_lattice_points(long long x1, long long y1, long long x2, long long y2) {
    long long dx = x2 - x1;
    return floor_sum(dx, y2 - y1, dx * y1, dx);
}

long long GCD(long long x, long long y) {
    return y ? GCD(y, x % y) : x;
}

// modint
template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) val += MOD;
    }
    constexpr int getmod() const { return MOD; }
    constexpr Fp operator - () const noexcept {
        return val ? MOD - val : 0;
    }
    constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
    constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
    constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
    constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
    constexpr Fp& operator += (const Fp& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator == (const Fp& r) const noexcept {
        return this->val == r.val;
    }
    constexpr bool operator != (const Fp& r) const noexcept {
        return this->val != r.val;
    }
    friend constexpr istream& operator >> (istream& is, Fp<MOD>& x) noexcept {
        is >> x.val;
        x.val %= MOD;
        if (x.val < 0) x.val += MOD;
        return is;
    }
    friend constexpr ostream& operator << (ostream& os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD>& r, long long n) noexcept {
        if (n == 0) return 1;
        if (n < 0) return modpow(modinv(r), -n);
        auto t = modpow(r, n / 2);
        t = t * t;
        if (n & 1) t = t * r;
        return t;
    }
    friend constexpr Fp<MOD> modinv(const Fp<MOD>& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        return Fp<MOD>(u);
    }
};

// Binomial coefficient
template<class T> struct BiCoef {
    vector<T> fact_, inv_, finv_;
    constexpr BiCoef() {}
    constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) {
        init(n);
    }
    constexpr void init(int n) noexcept {
        fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1);
        int MOD = fact_[0].getmod();
        for(int i = 2; i < n; i++){
            fact_[i] = fact_[i-1] * i;
            inv_[i] = -inv_[MOD%i] * (MOD/i);
            finv_[i] = finv_[i-1] * inv_[i];
        }
    }
    constexpr T com(int n, int k) const noexcept {
        if (n < k || n < 0 || k < 0) return 0;
        return fact_[n] * finv_[k] * finv_[n-k];
    }
    constexpr T fact(int n) const noexcept {
        if (n < 0) return 0;
        return fact_[n];
    }
    constexpr T inv(int n) const noexcept {
        if (n < 0) return 0;
        return inv_[n];
    }
    constexpr T finv(int n) const noexcept {
        if (n < 0) return 0;
        return finv_[n];
    }
};

const int MOD = 1000000007;
using mint = Fp<MOD>;
BiCoef<mint> bc;

int main() {
    int N;
    cin >> N;
    vector<long long> X(N), Y(N);
    for (int i = 0; i < N; ++i) cin >> X[i] >> Y[i];
    vector<int> ids(N);
    iota(ids.begin(), ids.end(), 0);
    sort(ids.begin(), ids.end(), [&](int i, int j) {
        return Y[i]*X[j] < Y[j]*X[i];
    });
    mint res = 0;
    long long sy = 0;
    for (auto i : ids) {
        res -= mint(X[i]) * mint(sy);
        res -= num_lattice_points(0, 0, X[i], Y[i]);
        res += GCD(X[i], Y[i]);
        sy += Y[i];
    }
    reverse(ids.begin(), ids.end());
    sy = 0;
    for (auto i : ids) {
        res += mint(X[i]) * mint(sy);
        res += num_lattice_points(0, 0, X[i], Y[i]);
        sy += Y[i];
    }
    cout << res << endl;
}