最古の遺跡 - 2007年日本情報オリンピック本選問題3

問題へのリンク公式解説へのリンク AOJ AtCoder

難易度は5（JOI非公式難易度表 Ver. 1.9.9.5）。

キーワード

幾何・ベクトル

問題概要

二次元座標の格子点 $(x, y)$ に $N$ 本の柱が立っている。 4本の柱からなる正方形のうち、面積が最大のものの面積を出力せよ。

$1 \leq N \leq 3000$
$0 \leq x,y \leq 5000$

例
画像はhttps://www.ioi-jp.org/joi/2006/2007-ho-prob_and_sol/2007-ho.pdfより引用。アクセス日は2020年11月29日。

解法

まず愚直に考える。 4本の柱の組みを全探索して、正方形が作れるか、作れるならその面積を調べればよい。しかし時間計算量が $O(N ^ 4)$ となりTLEする。

2本の柱を最初に選ぼう。 2本の柱の位置が決まっていれば、その柱で正方形を作るために必要な残り2本の柱の位置は3通りに定まる。 これを利用して計算量を落とそう。

柱 $T_0, T_1$ を最初に選ぶ柱、柱 $T_2, T_3$ を残りの柱にする。 $T_2, T_3$ の位置は次図の3通りとなる。

f:id:hotarunx:20201129210439j:plain — T2, T3の位置

青🔵の場合

$T_2 = T_0 + \frac{1}{2} (T_1 - T_0) + R(-\frac{1}{2}\pi) (T_1 - T_0)$

$T_3 = T_1 + \frac{1}{2} (T_1 - T_0) + R(\frac{1}{2}\pi) (T_1 - T_0)$

赤💖の場合

$T_2 = T_0 + R(- \frac{1}{2}\pi) (T_1 - T_0)$

$T_3 = T_1 + R(- \frac{1}{2}\pi) (T_1 - T_0)$

緑🟩の場合

$T_2 = T_0 + R(\frac{1}{2}\pi) (T_1 - T_0)$

$T_3 = T_1 + R(\frac{1}{2}\pi) (T_1 - T_0)$

$T_2, T_3$ を計算して、その位置に柱があるかどうか高速に確認できたらよさそうだ。

柱の位置は？

柱の位置は必ず格子点。また $x,y \leq 5000$ である。 サイズ5000x5000の二次元配列gridを用意して、そこに柱があるかどうか管理する。

計算量

$H=max(y), W=max(x)$ とする。 計算量は $O(N ^ 2 + HW)$ である。

実装

言語はGCC C++17。

コード

// 最古の遺跡 - 2007年日本情報オリンピック本選問題3
// https://atcoder.jp/contests/joi2007ho/tasks/joi2007ho_c
// https://onlinejudge.u-aizu.ac.jp/challenges/sources/JOI/Final/0518?year=2007

// AtCoderとAOJはテストケースの形式が異なる
// ATCODERが定義されていればAtCoderのテストケースの形式に対応する

// オンラインジャッジがAtCoderならばATCODERを定義する
#define ATCODER

#include <array>
#include <complex>
#include <iostream>
#include <vector>
using namespace std;
constexpr int MAXX = 5000;
constexpr double EPS = 1e-3;
using GRID = array<array<bool, MAXX + 1>, MAXX + 1>;
constexpr bool between(int x, int lb, int ub) { return lb <= x && x <= ub; }

// grid[x][y]: グリッド(x,y)に柱がある？
GRID grid{};

// tに柱はある？
bool isTowerExist(const GRID &grid, complex<double> t);

// solve a testcase
void solve(int);

int main() {
    cin.tie(0);
    ios::sync_with_stdio(0);

    int n;

#ifdef ATCODER
// AtCoder テストケースは1つ
    cin >> n;
    solve(n);
#else
// AOJ テストケースは複数 n=0ならばテスト終了
    while (1) {
        cin >> n;
        if (n == 0) break;
        solve(n);
    }
#endif
}

// 1テストケースを解く
void solve(int n) {
    // 入力受け取り
    vector<complex<double>> towers(n);
    for (auto &&t : towers) {
        int x, y;
        cin >> x >> y;
        t = {(double)x, (double)y};
    }
    // gridを初期化
    for (auto &&row : grid) {
        row.fill(0);
    }
    // gridを更新
    for (auto &&t : towers) {
        grid[(int)t.real()][(int)t.imag()] = true;
    }

    double largest_area = 0;

    // 2本の柱t0, t1を全通り選ぶ
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            // 正方形を作るのに必要な柱の位置t2, t3と呼ぶ
            const complex<double> t0 = towers[i], t1 = towers[j];

            // t2, t3は3通り
            // t2, t3に柱があるか確認して、あれば最大面積を更新する

            complex<double> t2, t3;
            constexpr complex<double> c22 = complex<double>(2, 2);

            // その1
            t2 = t0 + (t1 - t0) / c22 + 1.0i * (t1 - t0) / c22;
            t3 = t0 + (t1 - t0) / c22 - 1.0i * (t1 - t0) / c22;
            if (isTowerExist(grid, t2) && isTowerExist(grid, t3)) {
                const double area = imag(conj(t2 - t0) * (t3 - t0));
                largest_area = max(area, largest_area);
            }

            // その2
            t2 = t0 - 1.0i * (t1 - t0);
            t3 = t1 - 1.0i * (t1 - t0);
            if (isTowerExist(grid, t2) && isTowerExist(grid, t3)) {
                const double area = imag(conj(t2 - t0) * (t3 - t0));
                largest_area = max(area, largest_area);
            }

            // その3
            t2 = t0 + 1.0i * (t1 - t0);
            t3 = t1 + 1.0i * (t1 - t0);
            if (isTowerExist(grid, t2) && isTowerExist(grid, t3)) {
                const double area = imag(conj(t2 - t0) * (t3 - t0));
                largest_area = max(area, largest_area);
            }
        }
    }

    cout << (long long)largest_area << "\n";
}

bool isTowerExist(const GRID &grid, complex<double> t) {
    // 座標が整数でないならNG roundと元の数の差がEPS以下なら整数ということにする
    if (abs(t.real() - round(t.real())) > EPS) return false;
    if (abs(t.imag() - round(t.imag())) > EPS) return false;
    t.real(round(t.real()));
    t.imag(round(t.imag()));

    // 座標がグリッド外ならNG
    if (!between(t.real(), 0, MAXX) || !between(t.imag(), 0, MAXX)) return false;

    return grid[(int)t.real()][(int)t.imag()];
}

実装で気をつけた点

二次元ベクトル $T_i$ はcomplex<double>で扱う。 complex<double>はメンバ変数や演算子が定義されてるから二次元座標を扱うときに使える。次の式も次のコードで表せる。

$T_2 = T_0 + R(\frac{1}{2}\pi) (T_1 - T_0)$ $T_3 = T_1 + R(\frac{1}{2}\pi) (T_1 - T_0)$

constexpr complex<double> c22 = complex<double>(2, 2);
t2 = t0 + 1.0i * (t1 - t0);
t3 = t1 + 1.0i * (t1 - t0);

complex<double>同士の加算減算ができるよ。
complex<double>とスカラの乗算・除算は用意されていないので代わりに実部と虚部が同じcomplex<double>を使う。
1.0iは実部が0、虚部が1.0のcomplex<double>リテラル。
ベクトルの90度回転は虚数単位を書けることで表現する。

doubleで表される座標が整数なのか調べたい。 round（最も近い整数に丸めた）した数と元の数の差がEPS未満ならば整数ということにする。

if (abs(t.real() - round(t.real())) > EPS) return false;
if (abs(t.imag() - round(t.imag())) > EPS) return false;
t.real(round(t.real()));
t.imag(round(t.imag()));