G - Conquest 解説 by en_translator

比較的簡潔な実装について

公式解説の「包絡線の性質を利用する解法」は、直線を傾きでソートしてスタックに入れつつ順に処理することで、簡潔に実装することができます。それぞれの線分は端点の \(x\) 座標が \(0\)\(1\)\(y\) 座標が(縦に \(10^6\) 倍すれば)\(0\) 以上 \(10^6\) 以下の整数なので、傾きは絶対値が \(10^6\) 以下の整数となることから正確にソート可能で、交点の \(x\) 座標は倍精度浮動小数点で計算しても大小関係は損なわれません。

後半パートに関しては、\((V_{i,1},V_{i,2},V_{i,3})\)\(O(1)\) 種類(正確には \(7\) 種類)しかないことから、平面 \[ z = V _ {i,1} x + V _ {i,2} y + V _ {i,3} (1-x-y) \] および \(x = 0, y = 0, x + y = 1\) の合計たかだか \(10\) 平面のうち \(3\) 平面の交点の \((x, y)\) 座標を全列挙し、そのうち \(\Delta_2\) に含まれる点での下包絡面の \(z\) 座標の最大値を求めれば答えが求められます。

実装例 (Rust)

use itertools::Itertools;
use ordered_float::OrderedFloat as F;
use proconio::input;

fn main() {
    let f = |&x: &f64| F(x);
    input!(t: usize);
    for _ in 0..t {
        input!(n: usize, x: [f64; 3 * n]);
        let mut v = vec![[0.; 3]; n];
        for j in 0..3 {
            let lines = || {
                (0..n).map(|i| {
                    let (pk1, pk2) = (x[i * 3 + j], x[i * 3 + (j + 1) % 3]);
                    (pk1 - pk2, pk2, i)
                })
            };
            let (min, s) = min_lines(lines());
            (0..n).for_each(|i| v[i][j] = min);
            for i in s {
                v[i][j] = min_lines(lines().filter(|x| x.2 != i)).0;
            }
        }
        v.sort_by_key(|x| x.map(F));
        v.dedup();
        let planes = (v.iter().map(|&[p, q, r]| ([r - p, r - q, 1.], r)))
            .chain([([1., 0., 0.], 0.), ([0., 1., 0.], 0.), ([1., 1., 0.], 1.)])
            .collect_vec();
        let ans = (0..planes.len())
            .flat_map(|i| (0..i).flat_map(move |j| (0..j).map(move |k| [i, j, k])))
            .filter_map(|ijk| {
                let [x, y, _] = intersection(ijk.map(|i| planes[i]))?;
                (x >= -1e-6 && y >= -1e-6 && x + y <= 1. + 1e-6).then(|| {
                    (0..planes.len() - 3)
                        .map(|i| {
                            let ([a, b, _], d) = planes[i];
                            d - a * x - b * y
                        })
                        .min_by_key(f)
                        .unwrap()
                })
            })
            .max_by_key(f)
            .unwrap();
        println!("{:.30}", ans / 1e6);
    }
}

fn min_lines(lines: impl Iterator<Item = (f64, f64, usize)>) -> (f64, Vec<usize>) {
    let mut lines = lines.collect_vec();
    lines.sort_by_key(|x| (F(x.0), F(-x.1)));
    lines.dedup_by_key(|x| x.0);

    let mut stack = vec![];
    'line: for abi @ (a, b, _) in lines {
        while let Some(&((a0, b0, _), x0)) = stack.last() {
            let x = (b0 - b) / (a - a0);
            if 1. <= x {
                continue 'line;
            } else if x <= x0 {
                stack.pop();
            } else {
                stack.push((abi, x));
                continue 'line;
            }
        }
        stack.push((abi, 0.));
    }

    if let Some(&((0.0.., b, i), _)) = stack.first() {
        (b, vec![i])
    } else if let Some(&((a @ ..0., b, i), _)) = stack.last() {
        (a + b, vec![i])
    } else {
        (stack.iter().tuple_windows())
            .find_map(|(&((aa, _, i), _), &((a, b, j), x))| {
                (aa.signum() * a.signum() <= 0.).then(|| (a * x + b, vec![i, j]))
            })
            .unwrap()
    }
}

fn intersection(planes: [([f64; 3], f64); 3]) -> Option<[f64; 3]> {
    let cross =
        |[[a, b, c], [x, y, z]]: [[f64; _]; _]| [b * z - c * y, c * x - a * z, a * y - b * x];
    let cross = |ij: [usize; _]| cross(ij.map(|i| planes[i % 3].0));
    let dot = |[[a, b, c], [x, y, z]]: [[f64; _]; _]| a * x + b * y + c * z;
    let det = dot([planes[0].0, cross([1, 2])]);
    (det.abs() > 1e-8).then(|| {
        let mut ret = [0.; 3];
        for i in 0..3 {
            let cross = cross([i + 1, i + 2]);
            for j in 0..3 {
                ret[j] += planes[i].1 * cross[j];
            }
        }
        ret.map(|x| x / det)
    })
}

投稿日時:
最終更新: