aoc/year2019/

day10.rs

//! # Monitoring Station
//!
//! Integer only solution, avoiding floating point or trigonometric lookups such as [`atan2`].
//!
//! ## Part One
//!
//! We compare each pair of points, first subtracting the current asteroid to get a relative vector.
//! Since all coordinates are integers we can check for multiple points on the same line by
//! reducing the vector by its [greatest common divisor](https://en.wikipedia.org/wiki/Greatest_common_divisor).
//! For example, looking from the origin `(0, 0)`, the points `(3, 5)`, `(6, 10)` and `(21, 35)`
//! are all on the same line, with gcds of 1, 2 and 7 respectively.
//!
//! For each point we build a set of previously seen values. Since we can see at most one asteroid
//! in a given direction, if a vector is already in the set then we ignore it. The final size of
//! the set is the number of visible asteroids.
//!
//! To speed things up a little, we rely on the fact that if asteroid `a` can see `b`, then `b`
//! can see `a`. The solution is still `O(n²)` complexity but we reduce the work by half. This
//! works by implicitly processing asteroids in the same order as the input, from left to right and
//! top to bottom, so that nearest asteroids are always encountered first.
//!
//! ## Part Two
//!
//! The key insight is that we only need the *relative* angle between two vectors to sort them
//! in clockwise order. The [vector cross product](https://en.wikipedia.org/wiki/Cross_product)
//! of `a` and `b` can be defined as `|a||b|sinθ` where `θ` is the angle between the vectors.
//! `θ` will be negative if anti-clockwise, zero if on the same line or positive if clockwise.
//!
//! This works for angles up to 90°. To handle the complete 360° rotation, we first order points
//! by [quadrant](https://en.wikipedia.org/wiki/Quadrant_(plane_geometry)) then by relative angle.
//!
//! Finally, we also order points from nearest to furthest, so that the total ordering is:
//! 1. Quadrant
//! 2. Clockwise angle
//! 3. Distance
//!
//! This results in something like (where letters indicate angle and numbers indicate distance):
//!
//! `a1 a2 a3 b1 c1 c2 c3 c4 d1 d2`
//!
//! We want to process asteroids in the order:
//!
//! `a1 b1 c1 d1 a2 c2 d2 a3 c3 c4`
//!
//! We do this by first numbering the position within each group, then numbering the group and
//! sorting a second time in this order.
//!
//! [`atan2`]: f64::atan2
use std::cmp::Ordering;

use crate::util::grid::*;
use crate::util::math::*;
use crate::util::point::*;

type Input = (i32, i32);

pub fn parse(input: &str) -> Input {
    let grid = Grid::parse(input);
    let width = grid.width;
    let height = grid.height;

    // Convert asteroids to `Point` objects.
    let mut points = Vec::new();

    for y in 0..height {
        for x in 0..width {
            let point = Point::new(x, y);
            if grid[point] == b'#' {
                points.push(point);
            }
        }
    }

    // Find asteroid with the highest visibility.
    let mut seen = Grid::new(2 * width, height, 0);
    let mut cache = seen.same_size_with(ORIGIN);
    let mut visible = vec![0; points.len()];
    let mut max_visible = 0;
    let mut max_index = 0;

    for i in 0..(points.len() - 1) {
        for j in (i + 1)..points.len() {
            let delta = points[j] - points[i];
            // Because points was populated in order from left to right and top to bottom, delta.y
            // will be non-negative. Although delta.x might be negative, shifting it by
            // width gives us a positive reference point that can then cache the mapping
            // to a representative index.
            let key = Point::new(width + delta.x, delta.y);

            if cache[key] == ORIGIN {
                // Key insight is that points on the same line are integer multiples of each other.
                let gcd = delta.x.abs().gcd(delta.y);
                let adjusted = Point::new(delta.x / gcd + width, delta.y / gcd);

                cache[key] = adjusted;
            }
            let index = cache[key];

            if seen[index] <= i {
                seen[index] = i + 1;
                visible[i] += 1;
                visible[j] += 1;
            }
        }

        if visible[i] > max_visible {
            max_visible = visible[i];
            max_index = i;
        }
    }

    // Remove our new base of operations, then sort remaining asteroids in clockwise order to
    // group by angle.
    let station = points.swap_remove(max_index);
    points.iter_mut().for_each(|p| *p -= station);
    points.sort_unstable_by(|&a, &b| clockwise(a, b));

    // Sort asteroids a second time, first by order within the group, then the group's order.
    let mut groups = Vec::with_capacity(points.len());
    let mut first = 0;
    let mut second = 0;

    groups.push((first, second, 0));

    for i in 1..points.len() {
        if angle(points[i], points[i - 1]) == Ordering::Greater {
            first = 0;
            second += 1;
        } else {
            first += 1;
        }
        groups.push((first, second, i));
    }

    // The 200th asteroid is at index 199. Quickselect is faster than a full sort.
    let index = groups.select_nth_unstable(199).1.2;
    let target = station + points[index];
    let result = 100 * target.x + target.y;

    (max_visible, result)
}

pub fn part1(input: &Input) -> i32 {
    input.0
}

pub fn part2(input: &Input) -> i32 {
    input.1
}

/// The [`then`] method chains [`Ordering`] results.
///
/// [`then`]: Ordering::then
fn clockwise(point: Point, other: Point) -> Ordering {
    quadrant(point)
        .cmp(&quadrant(other))
        .then(angle(point, other))
        .then(distance(point).cmp(&distance(other)))
}

/// Divide points into one of four quadrants. For points exactly on an axis, for example (1, 0)
/// or (-5, 0) we can choose either adjacent quadrant as long as we're consistent.
fn quadrant(point: Point) -> i32 {
    match (point.x >= 0, point.y >= 0) {
        (true, false) => 0,
        (true, true) => 1,
        (false, true) => 2,
        (false, false) => 3,
    }
}

/// Positive if clockwise, zero if on the same line, negative if anti-clockwise.
fn angle(point: Point, other: Point) -> Ordering {
    (other.x * point.y).cmp(&(other.y * point.x))
}

/// Euclidean distance squared. No need to take square root since we're only interested
/// in the *relative* distance.
fn distance(point: Point) -> i32 {
    point.x * point.x + point.y * point.y
}