mgrs-convert 0.1.0

//! Grid geometry: MGRS coverage limits, the boundary nudge, and the single
//! float-to-integer truncation that produces 100 km indices and cell digits.
//
// This module performs the bounded float arithmetic that turns validated
// coordinates into integer grid fields, so `clippy::arithmetic_side_effects`
// is not denied here; see the errors section of the behavioral specification.
// The digit integers themselves are i64 throughout (mechanism M1).

/// The boundary epsilon, ldexp(1, -(53 - 25)) = 2^-28 m: the nudge that pulls a
/// coordinate sitting exactly on an excluded upper boundary just inside it.
/// Mechanism M4: written as a power of two next to the test pinning its value.
pub(crate) fn boundary_eps() -> f64 {
    2f64.powi(-28)
}

/// Cell size in metres at each precision, 10^(5-p).
pub(crate) const CELL: [f64; 12] = [
    1e5, 1e4, 1e3, 1e2, 1e1, 1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6,
];

/// 10^(11-p) as i64, the divisor that truncates micrometre offsets to p digits.
pub(crate) const POW10_11MP: [i64; 12] = [
    100_000_000_000,
    10_000_000_000,
    1_000_000_000,
    100_000_000,
    10_000_000,
    1_000_000,
    100_000,
    10_000,
    1_000,
    100,
    10,
    1,
];

/// Micrometres per 100 km tile.
const TILE_UM: i64 = 100_000_000_000;

/// Strict MGRS coverage, captured from the oracle at stage 0. Upper edges are
/// included and eps-nudged inside. Lower edges are closed.
pub(crate) const UTM_MGRS_EASTING: (f64, f64) = (100_000.0, 900_000.0);
/// Northing relative to the equator.
pub(crate) const UTM_MGRS_NORTHING_ABS: (f64, f64) = (-9_000_000.0, 9_500_000.0);
pub(crate) const UPS_NORTH_MGRS: (f64, f64) = (1_300_000.0, 2_700_000.0);
pub(crate) const UPS_SOUTH_MGRS: (f64, f64) = (800_000.0, 3_200_000.0);

/// One axis reduced to grid integers: the 100 km tile index and the offset
/// within the tile in micrometres.
///
/// The float is truncated exactly once, at micrometre resolution, and every
/// coarser digit count is an integer division of that one result, so the
/// twelve precision strings nest by construction (invariant 4).
pub(crate) fn axis_fields(value: f64) -> (i64, i64) {
    let um = (value * 1e6).floor() as i64;
    let tile = um.div_euclid(TILE_UM);
    let offset = um.rem_euclid(TILE_UM);
    (tile, offset)
}

/// Digits of an in-tile micrometre offset at p digit pairs.
pub(crate) fn digits_at(offset_um: i64, p: u8) -> i64 {
    let div = POW10_11MP.get(usize::from(p)).copied().unwrap_or(1);
    offset_um.div_euclid(div)
}

/// Cell size in metres at p digit pairs.
pub(crate) fn cell_size(p: u8) -> f64 {
    CELL.get(usize::from(p)).copied().unwrap_or(1.0)
}

/// MGRS latitude band index: floor(lat / 8) clamped to [-10, 9]. Division by 8
/// is exact, so a latitude one ulp from a band edge cannot round across it,
/// with one exception: a negative subnormal latitude underflows the division
/// to -0.0, where the mathematical floor is -1. Latitude 0 and -0.0 are both
/// band index 0, band N (`-0.0 < 0.0` is false, so no sign bit is consulted).
pub(crate) fn band_index(lat: f64) -> i64 {
    let f = (lat / 8.0).floor();
    let ib = if lat < 0.0 && f == 0.0 { -1 } else { f as i64 };
    ib.clamp(-10, 9)
}

/// The excluded-upper-edge nudge for UTM, applied before truncation: easting
/// exactly 900,000, frame northing exactly at the folded grid top, and, in the
/// southern frame only, a northing exactly on the equator, each move inside by
/// [`boundary_eps`]. Lower edges are closed and are not nudged. The northing
/// compares in the frame against the folded grid top, matching [`utm_extent`].
pub(crate) fn utm_nudge(north: bool, e: f64, n_frame: f64) -> (f64, f64) {
    let (_, ehi) = UTM_MGRS_EASTING;
    let (_, nhi) = UTM_MGRS_NORTHING_ABS;
    let fold = if north { 0.0 } else { 10_000_000.0 };
    let e2 = if e == ehi { e - boundary_eps() } else { e };
    let n2 = if n_frame == nhi + fold || (!north && n_frame == fold) {
        n_frame - boundary_eps()
    } else {
        n_frame
    };
    (e2, n2)
}

/// UTM MGRS coverage check. Returns the nudged coordinate for a value inside
/// the strict coverage (upper edges included) and `None` for the transfer
/// ring.
///
/// The northing compares in the caller's frame against the folded limits,
/// which are exact integer sums. Subtracting the fold from the coordinate
/// instead cancels: a southern-frame northing just under 1,000,000 m carries
/// ulps finer than the difference near -9,000,000 m can represent, so the
/// fold rounded the deficit away, the check passed, and the same value then
/// truncated to row -91. The row derivation reads the frame value, so the
/// check must too.
pub(crate) fn utm_extent(north: bool, e: f64, n_frame: f64) -> Option<(f64, f64)> {
    let (elo, ehi) = UTM_MGRS_EASTING;
    if !(elo..=ehi).contains(&e) {
        return None;
    }
    let fold = if north { 0.0 } else { 10_000_000.0 };
    let (nlo, nhi) = UTM_MGRS_NORTHING_ABS;
    if !(nlo + fold..=nhi + fold).contains(&n_frame) {
        return None;
    }
    Some(utm_nudge(north, e, n_frame))
}

/// The excluded-upper-edge nudge for UPS, both axes.
pub(crate) fn ups_nudge(north: bool, e: f64, n: f64) -> (f64, f64) {
    let (_, hi) = if north {
        UPS_NORTH_MGRS
    } else {
        UPS_SOUTH_MGRS
    };
    let e2 = if e == hi { e - boundary_eps() } else { e };
    let n2 = if n == hi { n - boundary_eps() } else { n };
    (e2, n2)
}

/// UPS MGRS coverage check, per cap, captured from the oracle at stage 0.
pub(crate) fn ups_extent(north: bool, e: f64, n: f64) -> Option<(f64, f64)> {
    let (lo, hi) = if north {
        UPS_NORTH_MGRS
    } else {
        UPS_SOUTH_MGRS
    };
    if !(lo..=hi).contains(&e) || !(lo..=hi).contains(&n) {
        return None;
    }
    Some(ups_nudge(north, e, n))
}

/// Center of a cell along one axis: tile origin plus (digits + 1/2) cells.
pub(crate) fn cell_center_axis(tile: i64, digits: i64, p: u8) -> f64 {
    (tile as f64) * 100_000.0 + ((digits as f64) + 0.5) * cell_size(p)
}

/// Frame northing from an equator-relative northing.
pub(crate) fn frame_northing(n_abs: f64, north: bool) -> f64 {
    if north {
        n_abs
    } else {
        n_abs + 10_000_000.0
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn boundary_eps_pinned() {
        // Mechanism M4. The reference constant is 2^-28 m, about 3.7 nm; the
        // incumbent shipped 1024 m here.
        assert_eq!(boundary_eps(), 3.725_290_298_461_914e-9);
        assert_eq!(
            boundary_eps(),
            (2f64).powi(25 - f64::MANTISSA_DIGITS as i32)
        );
    }

    #[test]
    fn utm_extent_compares_in_the_frame() {
        // One f64 below the closed southern limit. The ~0.116 nm deficit is
        // finer than the spacing near -9,000,000 m, so folding the coordinate
        // rounds the deficit away. The first assertion pins that hazard.
        let just_under = f64::from_bits(1_000_000.0_f64.to_bits() - 1);
        assert_eq!(just_under - 10_000_000.0, -9_000_000.0);
        assert!(utm_extent(false, 500_000.0, just_under).is_none());
        assert!(utm_extent(false, 500_000.0, 1_000_000.0).is_some());
        assert!(utm_extent(true, 500_000.0, -9_000_000.0).is_some());
        let below_north = f64::from_bits((-9_000_000.0_f64).to_bits() + 1);
        assert!(utm_extent(true, 500_000.0, below_north).is_none());
    }

    #[test]
    fn axis_fields_truncates_toward_negative_infinity() {
        // Binary-exact inputs so the micrometre product is exact.
        assert_eq!(axis_fields(585_664.25), (5, 85_664_250_000));
        // -2^-10 m sits just below zero and lands in the tile above -100 km.
        assert_eq!(axis_fields(-0.000_976_562_5), (-1, 99_999_999_023));
        assert_eq!(axis_fields(0.0), (0, 0));
    }

    #[test]
    fn digits_nest_by_construction() {
        let (_, off) = axis_fields(585_664.121_374_411);
        for p in 1..=11u8 {
            let fine = digits_at(off, p);
            let coarse = digits_at(off, p - 1);
            assert_eq!(fine / 10, coarse, "nesting at p {p}");
        }
    }
}