Official
F - Integer Convex Hull Editorial by KoD
点の座標が整数であるので、任意の \(3\) 点からなる三角形の面積は \(2\) 倍すると整数になります。凸包の面積は、ある点を固定し、凸包をその点を含む三角形に分解することによって計算することができます。「ある固定された点を必ず含むような凸多角形」を、三角形を一つずつ貼り合わせていくかのように構成することはできないでしょうか？
既知のアルゴリズムを適用できないか考えてみましょう。凸包を構成するアルゴリズムとしては、Monotone Chain と呼ばれるものが有名です。概略は以下のようになります :
- 点を \(x\) 座標の小さい順にソート する ( \(x\) 座標が同じものは \(y\) 座標の小さい順に並べる) 。ソート後において先頭にある点と末尾にある点は必ず凸包に含まれる。
- 点を順番に見ていき、前述の二点を結ぶように、上側凸包と下側凸包を構築する。
上の方法と同様に、まずは与えられた点をソートします。次に、凸多角形に含まれる左下の点 ( \(x\) 座標が最も小さいものの中で、\(y\) 座標が最も小さいもの) を固定し、\(O(N)\) 通り試すことにします。
固定された点を \(P_s\) とおきます。 \(P_s\) を一点にもつ三角形を貼り合わせていくようにして、上側凸包と下側凸包を構成しましょう。以下の動的計画法を考えます :
- \(\mathrm{upper}_{i, j, k}\)：頂点集合 \(S\) であって、その凸包が以下の条件を満たすものの総数。
- \(P_s\) から時計回りに点を見ていったとき、最後の \(2\) 点がそれぞれ \(P_i, P_j\) である。
- 全ての点が線分 \(P_sP_j\) の上側にある。
- 面積を \(2\) 倍した値を \(2\) で割ったあまりが \(k\) である。
- \(\mathrm{lower}_{i, j, k}\)：頂点集合 \(S\) であって、その凸包が以下の条件を満たすものの総数。
- \(P_s\) から反時計回りに点を見ていったとき、最後の \(2\) 点がそれぞれ \(P_i, P_j\) である。
- 全ての点が線分 \(P_s P_j\) の下側にある。
- 面積を \(2\) 倍した値を \(2\) で割ったあまりが \(k\) である。
各 \(j, k\) について \(\displaystyle \left( \sum_{i} \mathrm{upper}_{i, j, k} \right) \times \left( \sum_{i} \mathrm{lower}_{i, j, k} \right)\) を足し合わせたものが答えとなります。
前計算として、各三角形についてその面積と内部に含まれる点の個数を調べておく必要があることに注意してください。全体の計算量は \(O(N^4)\) となります。
実装例 (C++) :
#include <iostream>
#include <algorithm>
#include <atcoder/modint>
struct Point {
int x, y;
explicit Point(const int x = 0, const int y = 0): x(x), y(y) {}
bool operator < (const Point& other) const {
return x < other.x or (x == other.x and y < other.y);
}
Point operator - (const Point& other) const {
return Point(x - other.x, y - other.y);
}
};
int cross(const Point& p, const Point& q) {
return p.x * q.y - p.y * q.x;
}
using modint = atcoder::modint1000000007;
constexpr int MAX = 80;
int N;
Point P[MAX];
int inside[MAX][MAX][MAX], parity[MAX][MAX][MAX];
modint upper[MAX][MAX][2], lower[MAX][MAX][2];
modint pow2[MAX + 1];
// 三角形の面積を 2 倍した値
int area(const int i, const int j, const int k) {
return std::abs(cross(P[j] - P[i], P[k] - P[i]));
}
int main() {
std::cin >> N;
for (int i = 0; i < N; ++i) {
std::cin >> P[i].x >> P[i].y;
}
std::sort(P, P + N);
// 各三角形について、その内部に含まれる点の個数を計算する
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
for (int k = 0; k < N; ++k) {
if (i == j or j == k or k == i) {
continue;
}
parity[i][j][k] = area(i, j, k) % 2;
for (int l = 0; l < N; ++l) {
if (l != i and l != j and l != k and area(l, i, j) + area(l, j, k) + area(l, k, i) == area(i, j, k)) {
inside[i][j][k] += 1;
}
}
}
}
}
// 2 の冪乗を前計算しておく
pow2[0] = 1;
for (int i = 0; i < N; ++i) {
pow2[i + 1] = pow2[i] * 2;
}
modint ans = 0;
for (int must = N - 1; must >= 0; --must) {
// P[must] を左下の点とする凸多角形
for (int i = must; i < N; ++i) {
for (int j = must; j < N; ++j) {
for (int k = 0; k < 2; ++k) {
upper[i][j][k] = lower[i][j][k] = 0;
}
}
}
for (int i = must + 1; i < N; ++i) {
upper[must][i][0] = lower[must][i][0] = 1;
}
for (int i = must; i < N; ++i) {
for (int j = i + 1; j < N; ++j) {
for (int k = 0; k < 2; ++k) {
for (int l = j + 1; l < N; ++l) {
// 線分 P_iP_j と線分 P_jP_l の関係を調べる
if (cross(P[l] - P[j], P[j] - P[i]) > 0) {
// 新たに P_must, P_j, P_l を頂点にもつ三角形を貼り合わせる
upper[j][l][k ^ parity[must][j][l]] += upper[i][j][k] * pow2[inside[must][j][l]];
}
else {
lower[j][l][k ^ parity[must][j][l]] += lower[i][j][k] * pow2[inside[must][j][l]];
}
}
}
}
}
for (int j = must + 1; j < N; ++j) {
for (int k = 0; k < 2; ++k) {
modint up = 0, lo = 0;
for (int i = must; i < j; ++i) {
up += upper[i][j][k];
lo += lower[i][j][k];
}
ans += up * lo;
}
}
}
// 線分も凸多角形として数えてしまっているので、除いてから出力する
std::cout << (ans - N * (N - 1) / 2).val() << '\n';
return 0;
}
別解
凸包を構成するもう一つのアルゴリズムとして、Graham Scan と呼ばれるものがあります。これは、左下の点を中心に偏角ソートを行い、点を順番に見て凸包を構成するというものです。これを応用した動的計画法を行うことで、凸包を上下に分割することなく数え上げることができます。
実装例 (Rust) :
use proconio::{derive_readable, input};
use std::cmp::Ordering;
use std::ops::Sub;
#[derive_readable]
#[derive(Copy, Clone, Eq, PartialEq, Ord, PartialOrd)]
struct Point {
x: i32,
y: i32,
}
impl Sub for Point {
type Output = Self;
fn sub(self, other: Self) -> Self::Output {
Self {
x: self.x - other.x,
y: self.y - other.y,
}
}
}
impl Point {
fn cross(self, other: Self) -> i32 {
self.x * other.y - self.y * other.x
}
}
fn main() {
input! {
n: usize,
mut pts: [Point; n],
}
pts.sort();
let mut ans: u128 = 0;
for leftmost in 0..n {
let mut pts: Vec<_> = pts.iter().skip(leftmost).map(|&p| p - pts[leftmost]).collect();
pts.sort_by(|&a, &b| {
if a.cross(b) > 0 {
Ordering::Less
} else {
Ordering::Greater
}
});
let mut count = vec![vec![0; n]; n];
for i in 0..pts.len() {
for j in i + 1..pts.len() {
for k in i + 1..j {
if (pts[i] - pts[k]).cross(pts[j] - pts[k]) > 0 {
count[i][j] += 1;
}
}
}
}
let mut dp = vec![vec![[0; 2]; n]; n];
for j in 1..pts.len() {
dp[0][j][0] += 1;
}
for i in 0..pts.len() {
for j in i + 1..pts.len() {
for k in 0..2 {
for l in j + 1..pts.len() {
if (pts[j] - pts[i]).cross(pts[l] - pts[j]) > 0 {
dp[j][l][k ^ (pts[j].cross(pts[l]).rem_euclid(2) as usize)] +=
dp[i][j][k] * (1 << count[j][l]);
}
}
}
}
}
for i in 1..pts.len() {
for j in i + 1..pts.len() {
ans += dp[i][j][0];
}
}
}
println!("{}", ans % 1000000007);
}
類題
posted:
last update: