ARC 082 E: ConvexScore

Keywords

convex hull DP

Time Complexity

\(O(N^3)\)

汎用的な処理が多数含まれているのでメモ。

#include <iostream>
#include <algorithm>
#include <vector>
#include <complex>
#include <cmath>
#include <cassert>
#include <tuple>

using namespace std;

constexpr long long mod = 998244353;

using P = complex<long long>;

long long dot(P a, P b) {
  return real(conj(a) * b);
}

long long cross(P a, P b) {
  return imag(conj(a) * b);
}

long long gcd(long long a, long long b) {
  a = abs(a);
  b = abs(b);
  if (b == 0) return a;
  return gcd(b, a % b);
}

// Get minimum element pointing in the same direction.
P tomato(P a) {
  long long g = gcd(a.real(), a.imag());
  return P(a.real() / g, a.imag() / g);
}

// Lexical order (argument, norm)
bool argcomp(P p1, P p2) {
  bool a = p1.imag() > 0 || (p1.imag() == 0 && p1.real() >= 0);
  bool b = p2.imag() > 0 || (p2.imag() == 0 && p2.real() >= 0);
  if (a != b) return a;
  long long c = p2.real() * p1.imag();
  long long d = p1.real() * p2.imag();
  if (c != d) return c < d;
  return norm(p1) < norm(p2);
}

namespace std {
  bool operator<(const P &a, const P &b) {
    return argcomp(a, b);
  }
}

int main() {
  int n;
  cin >> n;

  vector<P> ps(n);
  for (int i = 0; i < n; i++) {
    int x, y;
    cin >> x >> y;
    ps[i] = P(x, y);
  }

  vector<pair<int, int>> es;
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
      if (i == j) continue;
      es.emplace_back(i, j);
    }
  }
  // If two vectors are pointing same direction,
  // these are arranged such that both aren't taken.
  sort(es.begin(), es.end(), [&](pair<int, int> a, pair<int, int> b) {
    P d = ps[a.second] - ps[a.first];
    P e = ps[b.second] - ps[b.first];
    if (cross(d, e) == 0 && dot(d, e) >= 0) {
      return dot(ps[a.second] - ps[a.first], ps[b.first] - ps[a.first]) < 0;
    } else {
      return d < e;
    }
  });

  // Count the points inside (NOT on) the triangle(i,j,k)
  static int in[200][200][200];
  for (int j = 0; j < n; j++) {
    vector<pair<P, int>> qs;
    for (int i = 0; i < n; i++) {
      if (i == j) continue;
      qs.emplace_back(ps[i] - ps[j], i);
    }
    sort(qs.begin(), qs.end());
    for (int i = 0; i < n; i++) {
      if (i == j) continue;
      P s = tomato(ps[j] - ps[i]);
      int f = lower_bound(qs.begin(), qs.end(), make_pair(s, 0)) - qs.begin();
      for (int l = 0; l < n - 1; l++) {
        int k = qs[(f + l) % (n - 1)].second;
        in[i][j][k] -= l + 1;
        in[j][k][i] -= l + 1;
        in[k][i][j] -= l + 1;
      }
    }
  }
  static int on[200][200];
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
      if (i == j) continue;
      for (int k = 0; k < n; k++) {
        if (i == k || j == k) continue;
        on[i][j] += cross(ps[k] - ps[i], ps[j] - ps[i]) == 0 && dot(ps[k] - ps[i], ps[j] - ps[i]) > 0 && dot(ps[k] - ps[j], ps[i] - ps[j]) > 0;
        in[i][j][k] = max(0, in[i][j][k] + n);
      }
    }
  }

  vector<long long> two(n + 1);
  two[0] = 1;
  for (int i = 1; i <= n; i++) {
    two[i] = two[i - 1] * 2 % mod;
  }

  // Easy
  long long ans = 0;
  for (int i = 0; i < n; i++) {
    vector<long long> dp1(n);
    vector<long long> dp2(n);
    for (auto e : es) {
      int u = e.first;
      int v = e.second;

      if (u == i) {
        assert(dp1[v] == 0);
        dp1[v] = two[on[u][v]];
      } else if (v == i) {
        dp2[v] += dp2[u];
      } else {
        int c = on[i][v] + on[u][v] + in[i][u][v];
        dp2[v] += dp1[u] * two[c];
        dp2[v] += dp2[u] * two[c];
      }
      dp2[v] %= mod;
    }
    ans += dp2[i];
    ans %= mod;
  }

  cout << ans << endl;
}