Hoi_koro memo

Editorialと異なる解法をしたら更新するblog

AtCoder codeFlyer2018 本戦 H 三角形と格子点

問題

H - 三角形と格子点

三角形ABCの3頂点を

  • A{\in[x_1,x_1+w)\times[y_1,y_1+h)}
  • B{\in[x_2,x_2+w)\times[y_2,y_2+h)}
  • C{\in[x_3,x_3+w)\times[y_3,y_3+h)}

となるようにとる.このようにして作られる{w^3h^3}個の三角形について三角形の内部に存在する格子点の個数の総和を求めよ.(mod {10^9+7})

制約

  • {0\leq X_i, Y_i \leq 10^{12}}
  • {1\leq w,h \leq 40000}

残りの制約はA,B,Cの位置関係について

x_a < x_c< x_b, y_a < y_b < y_c

が成り立つことを示している.

解法

ピックの定理から各三角形について

面積 = (内部の格子点数) + (周上の格子点数)/2-1

が成り立つ.これを{w^3h^3}個の三角形について和をとると

(面積の和) =(答) + (周上の格子点数の和)/2 - {w^3h^3}.

よって面積の和と周上の格子点数の和が分かれば解ける.

面積の和

A,B,Cの座標を

  • A{(x_1+w_1,y_1+h_1)}
  • B{(x_2+w_2,y_2+h_2)}
  • C{(x_3+w_3,y_3+h_3)}

とすると面積の和は

{\displaystyle{\sum_{w_1=0}^w\sum_{w_2=0}^w\sum_{w_3=0}^w\sum_{h_1=0}^h\sum_{h_2=0}^h\sum_{h_3=0}^h \frac{1}{2}}\left|\vec{AB}\times\vec{AC}\right| =\dots= \frac{1}{2}w^3h^3\left|(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\right|}.

面積の平均が\frac{1}{2}\left|(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\right|になることからもこれがわかる.

周上の格子点数の和

ある三角形ABCについて,

  • \vec{AB} = (x_2-x_1 + a, y_2-y_1+b)~~(-w+1\leq a \leq w-1, -h+1\leq b \leq h-1)となる三角形は(w-|a|)(h-|b|)hw個ある.
  • そのような三角形ABCの辺AB上(点Bを除く)にある格子点数は{\text gcd}(x_2-x_1 + a, y_2-y_1+b).

よって{w^3h^3}個の三角形について辺AB上(点Bを除く)にある格子点数の総和は以下のように求められる.

ただしS_{ik}=\{j~|~{\text gcd}(i,j) = k\}とする.

また,オイラー関数\phi(i)\sum_{i|n}\phi(i) = nを満たすことに注意.

\begin{align} &\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}\sum_{j = y_2-y_1 -h +1}^{y_2-y_1 +h -1} hw(w-|i-x_2+x_1|)(h-|j-y_2+y_1|){\text gcd}(i,j)\\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{k|i}\sum_{p\in S_{ik}}(h-|p-y_2+y_1|)k\\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{k|i}\sum_{p\in S_{ik}}(h-|p-y_2+y_1|)\sum_{l|k}\phi(l) \\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{k|i}\sum_{l|k}\phi(l)\sum_{p\in \cup_{l|m}S_{im}}(h-|p-y_2+y_1|) \\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{k|i}\sum_{l|k}\phi(l)\sum_{l|p}(h-|p-y_2+y_1|) \\ =&hw\sum_{i = x_2-x_1 -w +1}^{x_2-x_1 +w -1}(w-|i-x_2+x_1|)\sum_{l|i}\phi(l)\sum_{l|p}(h-|p-y_2+y_1|) \\ \end{align}

  • \displaystyle\sum_{l|p}(h-|p-y_2+y_1|) は等差数列の和なので, それぞれのlについてO(1)時間で求められる.
  • \phi(l)の計算のために,まず[ x_2-x_1 -w +1, x_2-x_1 +w )の範囲について,区間ふるいによる素数判定を行うと同時に素因数分解しておく.すると任意のi \in [ x_2-x_1 -w +1, x_2-x_1 +w )について iの約数pとオイラー関数の組(p, \phi (p))を約数の個数について線形時間で求められる.

同じ計算を辺BC, CAについても行えばよい.

code

#include <bits/stdc++.h>

using LL = long long;
constexpr LL LINF = 334ll << 53;
constexpr int INF = 15 << 26;
constexpr LL MOD = 1E9 + 7;

namespace Problem {
using namespace std;
template <long long Mod = MOD>
class Modint {
 public:
  long long v;

  Modint(const Modint &x) { v = x.v; }
  explicit Modint(int x) {
    if (x < 0 or Mod <= x) x %= Mod;
    if (x < 0) x += Mod;
    v = x;
  }
  explicit Modint(long long x) {
    if (x < 0 or Mod <= x) x %= Mod;
    if (x < 0) x += Mod;
    v = x;
  }
  explicit Modint(unsigned x) {
    if (Mod <= x) x %= Mod;
    v = x;
  }
  explicit Modint(unsigned long long x) {
    if (Mod <= x) x %= Mod;
    v = x;
  }
  template <typename T>
  explicit Modint(T x) {
    if (x < 0 or Mod <= x) x %= Mod;
    v = x;
  }
  Modint() : v(0) {}
  long long get() const { return v; }
  Modint mpow(Modint &n) const { return mpow(n.v); }
  Modint mpow(long long n) const {
    long long i = 1, p = v;
    Modint ret(1);
    while (i <= n) {
      if (i & n) ret *= p;
      i = (i << 1);
      p = (p * p) % Mod;
    }
    return ret;
  }
  Modint operator-() const { return (v ? Modint(Mod - v) : Modint(0)); }
  explicit operator int() const { return v; }
  explicit operator long long() const { return v; }
  Modint &operator+=(const Modint &a) {
    v = (v + a.v);
    if (v >= Mod) v -= Mod;
    return *this;
  }
  Modint &operator-=(const Modint &a) {
    v = (v - a.v);
    if (v < 0) v += Mod;
    return *this;
  }
  Modint &operator*=(const Modint &a) {
    v = (v * a.v) % Mod;
    return *this;
  }
  Modint &operator/=(const Modint &a) {
    v = (v * a.mpow(Mod - 2).v) % Mod;
    return *this;
  }
  template <class T>
  Modint &operator+=(const T &a) {
    v = (v + Modint(a).v);
    if (v >= Mod) v -= Mod;
    return *this;
  }
  template <class T>
  Modint &operator-=(const T &a) {
    v = (v - Modint(a).v);
    if (v < 0) v += Mod;
    return *this;
  }
  template <class T>
  Modint &operator*=(const T &a) {
    v = (v * Modint(a).v) % Mod;
    return *this;
  }
  template <class T>
  Modint &operator/=(const T &a) {
    v = (v * Modint(a).mpow(Mod - 2).v) % Mod;
    return *this;
  }
  friend Modint operator+(const Modint &a, const Modint &b) {
    return Modint(a) += b;
  }
  friend Modint operator-(const Modint &a, const Modint &b) {
    return Modint(a) -= b;
  }
  friend Modint operator*(const Modint &a, const Modint &b) {
    return Modint(a) *= b;
  }
  friend Modint operator/(const Modint &a, const Modint &b) {
    return Modint(a) /= b;
  }
  template <class T, class U>
  friend Modint operator+(const T &a, const U &b) {
    return Modint(a) += Modint(b);
  }
  template <class T, class U>
  friend Modint operator-(const T &a, const U &b) {
    return Modint(a) -= Modint(b);
  }
  template <class T, class U>
  friend Modint operator*(const T &a, const U &b) {
    return Modint(a) *= Modint(b);
  }
  template <class T, class U>
  friend Modint operator/(const T &a, const U &b) {
    return Modint(a) /= Modint(b);
  }
};
template <long long M>
ostream &operator<<(ostream &os, const Modint<M> m) {
  return os << m.v;
}
class SegmentSieve {
  long long lower, upper, range;
  vector<long long> primes;  // prime number in[2,sqrt(b)]
  vector<vector<pair<long long, int>>> prime_factor;
  vector<char> is_prime_seg;
  void sieve() {
    vector<char> isprime(range + 1, 1);
    vector<long long> table(upper - lower);
    iota(table.begin(), table.end(), lower);
    isprime[0] = isprime[1] = 0;
    for (long long i = 2; i <= range; i++) {
      if (isprime[i]) {
        primes.push_back(i);
        for (long long j = i * i; j <= range; j += i) isprime[j] = 0;
        for (long long j = lower - lower % i; j < upper; j += i) {
          if (j >= lower) {
            is_prime_seg[j - lower] = false;
            int ord = 0;
            while (table[j - lower] % i == 0) {
              table[j - lower] /= i;
              ord++;
            }
            prime_factor[j - lower].emplace_back(i, ord);
          }
        }
      }
    }
    for (int i = 0; i < (int)prime_factor.size(); ++i) {
      if (table[i] != 1) prime_factor[i].emplace_back(table[i], 1);
    }
  }

 public:
  SegmentSieve(long long a, long long b)
      : lower(a),
        upper(b),
        prime_factor(b - a),
        is_prime_seg(b - a, 1) {  //[a,b)
    range = max(b - a + 1, (long long)sqrt(b) + 1);
    sieve();
  };
  bool is_prime(long long x) {
    // assert(lower <= x && x < upper);
    return is_prime_seg[x - lower];
  }
  vector<pair<long long, long long>> divisor_euler(long long x) {
    // assert(lower <= x && x < upper);
    vector<pair<long long, long long>> ret = {{1ll, 1}};
    for (const auto p : prime_factor[x - lower]) {
      int sz = (int)ret.size();
      for (int i = 0; i < sz; ++i) {
        long long divisor = ret[i].first;
        long long euler_ = ret[i].second;
        for (int d = 1; d <= p.second; ++d) {
          divisor *= p.first;
          euler_ *= (d == 1 ? (p.first - 1) : p.first);
          ret.emplace_back(divisor, euler_);
        }
      }
    }
    return ret;
  }
};

class Solver {
 public:
  LL w, h;
  vector<LL> x, y;
  Solver() : x(3), y(3){};

  void solve() {
    for (int i = 0; i < 3; ++i) {
      cin >> x[i] >> y[i];
    }
    cin >> w >> h;
    Modint<> area(0), on_edge(0);
    area = Modint<>(h).mpow(3) * Modint<>(w).mpow(3) / 2 *
           (Modint<>(x[1] - x[0]) * Modint<>(y[2] - y[0]) -
            Modint<>(x[2] - x[0]) * Modint<>(y[1] - y[0]));
    on_edge += calc(x[0], y[0], x[1], y[1]);
    on_edge += calc(x[0], y[0], x[2], y[2]);
    on_edge += calc(x[1], y[1], x[2], y[2]);

    cout << area + Modint<>(h).mpow(3) * Modint<>(w).mpow(3) - on_edge / 2
         << endl;
  }
  Modint<> calc_sum(LL k, LL center, LL width) {
    Modint<> ret(0);
    //等差数列の和の部分

    long long r = center - center % k;
    long long l = (center - width + k - 1) / k * k;
    ret += Modint<>((r - l) / k + 1) *
           Modint<>(h - abs(center - r) + h - abs(center - l)) / 2;

    long long r2 = center + width - (center + width) % k;
    long long l2 = (center + 1 + k - 1) / k * k;
    ret += Modint<>((r2 - l2) / k + 1) *
           Modint<>(h - abs(center - r2) + h - abs(center - l2)) / 2;

    return ret;
  }
  Modint<> calc(LL a, LL b, LL c, LL d) {
    if (a > c) swap(a, c);
    if (b > d) swap(b, d);
    SegmentSieve segs(c - a - w + 1, c - a + w);
    Modint<> ret(0);
    /*for (LL i = c - a - w + 1; i < c - a + w; ++i) {
      Modint<> tmp(0);
      for (LL j = d - b - h + 1; j < d - b + h; ++j) {
        tmp += Modint<>(__gcd(i, j)) * (h - abs(j - d + b));
      }
      ret += Modint<>(h * w) * tmp * (w - abs(i - c + a));
    }*/
    for (LL i = c - a - w + 1; i < c - a + w; ++i) {
      Modint<> tmp(0);
      for (auto p : segs.divisor_euler(i)) {
        auto div = p.first;
        auto euler = p.second;
        tmp += calc_sum(div, d - b, h) * euler;
      }
      ret += Modint<>(h * w) * tmp * (w - abs(i - c + a));
    }
    return ret;
  }
};
}  // namespace Problem

int main() {
  std::cin.tie(0);
  std::ios_base::sync_with_stdio(false);

  Problem::Solver sol;
  sol.solve();
  return 0;
}

note

知らなかったこと

  • \sum_{i|n}\phi(i) = n
  • 長さ 𝑂(𝐻) の区間に含まれる整数の約数の個数の和は 𝑂(𝐻 log 𝐻) (解説より)

特にオイラー関数の性質を知らなかったので,自力でACは無理だった.