oxiproj-core 0.1.2

//! Pure-Rust interval arithmetic with guaranteed (outward-rounded) enclosures.
//!
//! This module provides [`Interval`], a machine-interval type over `f64` that
//! computes a *rigorous enclosure* of the true mathematical result of every
//! operation. Given input intervals that each contain some (unknown) real
//! value, the output interval is guaranteed to contain the real result of the
//! same expression evaluated at those values. This underpins T7.4
//! (guaranteed-bounds transforms): running a projection over point-intervals
//! yields a certified box that provably contains the exact projected point.
//!
//! Unlike the SMT-solver route (which OxiZ 0.2.4 cannot yet serve — it has no
//! floating-point interval theory), this is a standalone, dependency-free
//! implementation that needs only stable `f64::next_up` / `f64::next_down`
//! (stable since Rust 1.85; MSRV here is 1.89) for outward rounding.
//!
//! # Rounding model (why the enclosures are valid)
//!
//! Every IEEE-754 basic operation (`+`, `-`, `*`, `/`, `sqrt`) is *correctly
//! rounded*: the computed `f64` differs from the exact real result by at most
//! half a ULP. Widening the low endpoint by one [`f64::next_down`] and the high
//! endpoint by one [`f64::next_up`] therefore brackets the exact result — the
//! standard "rounded-outward" construction. Basic ops here use a **1-ULP**
//! outward margin.
//!
//! Transcendental functions (`sin`, `cos`, `exp`, `ln`, `atan`, …) are *not*
//! guaranteed correctly rounded by `std`/`libm`, but high-quality libraries are
//! accurate to well under 1 ULP for the argument ranges used in geodesy. To
//! stay rigorous under that assumption we widen transcendental results by
//! [`TRANSCEND_ULP`] ULPs on each side. This is the conventional pragmatic
//! contract for interval arithmetic layered on a non-correctly-rounded libm and
//! is documented per function below.
//!
//! # `no_std`
//!
//! The type is `no_std`-compatible: `f64::next_up`/`next_down` live in `core`,
//! and every transcendental is routed through [`crate::math_compat`], which
//! dispatches to `std` or `libm` depending on the `no_std` feature.

use crate::math_compat as mc;
use core::fmt;
use core::ops::{Add, Div, Mul, Neg, Sub};

/// Outward-widening margin, in ULPs, applied to transcendental results.
///
/// Basic arithmetic (`+ - * /`, `sqrt`, `recip`, `powi`) is correctly rounded
/// and needs only a 1-ULP margin; the transcendental functions delegate to a
/// libm whose worst-case error is a few tenths of a ULP for our ranges, so a
/// 4-ULP margin gives a comfortable rigorous cushion.
pub const TRANSCEND_ULP: u32 = 4;

/// `f64` value nearest 2π. An interval whose width reaches `TWO_PI` certainly
/// spans at least one full sine/cosine period, so it is collapsed to `[-1, 1]`;
/// for narrower intervals at most one extremum of each sign can lie inside,
/// which the grid-point probe in [`contains_grid_point`] resolves exactly.
const TWO_PI: f64 = core::f64::consts::TAU;
const PI: f64 = core::f64::consts::PI;
const HALF_PI: f64 = core::f64::consts::FRAC_PI_2;

/// Widen `x` outward toward −∞ by `ulps` steps (NaN collapses to −∞).
#[inline]
fn wide_down(x: f64, ulps: u32) -> f64 {
    if x.is_nan() {
        return f64::NEG_INFINITY;
    }
    let mut r = x;
    let mut i = 0;
    while i < ulps {
        r = r.next_down();
        i += 1;
    }
    r
}

/// Widen `x` outward toward +∞ by `ulps` steps (NaN collapses to +∞).
#[inline]
fn wide_up(x: f64, ulps: u32) -> f64 {
    if x.is_nan() {
        return f64::INFINITY;
    }
    let mut r = x;
    let mut i = 0;
    while i < ulps {
        r = r.next_up();
        i += 1;
    }
    r
}

/// Rigorous hull of four candidate values: `(min, max)`, with any NaN
/// collapsing the hull to the whole real line (a safe over-approximation).
#[inline]
fn hull4(a: f64, b: f64, c: f64, d: f64) -> (f64, f64) {
    let mut lo = f64::INFINITY;
    let mut hi = f64::NEG_INFINITY;
    for v in [a, b, c, d] {
        if v.is_nan() {
            return (f64::NEG_INFINITY, f64::INFINITY);
        }
        if v < lo {
            lo = v;
        }
        if v > hi {
            hi = v;
        }
    }
    (lo, hi)
}

/// Conservative test: does `base + k*period` (for some integer `k`) fall inside
/// `[a, b]`? The caller guarantees `b - a < period` so at most one grid point
/// can lie inside; we probe the three integers nearest the midpoint and apply a
/// small outward guard. The test is intentionally biased to **over-report**
/// (return `true` when a grid point is merely *near* the boundary): for the
/// trig functions this only ever widens the result to the extremal `±1` /
/// unbounded value, which preserves the enclosure guarantee.
#[inline]
fn contains_grid_point(a: f64, b: f64, base: f64, period: f64) -> bool {
    let mid = 0.5 * (a + b);
    let kf = mc::round((mid - base) / period);
    let guard = 1e-9 * if period > 1.0 { period } else { 1.0 };
    for dk in [-1.0, 0.0, 1.0] {
        let x = base + (kf + dk) * period;
        if x >= a - guard && x <= b + guard {
            return true;
        }
    }
    false
}

/// A closed real interval `[lo, hi]` with rigorous outward-rounded arithmetic.
///
/// Invariant: `lo <= hi` and neither endpoint is NaN. Constructors normalise
/// any violating input to [`Interval::whole`] (the entire real line), which is
/// the safe over-approximation.
#[derive(Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Interval {
    lo: f64,
    hi: f64,
}

impl Interval {
    /// Construct `[lo, hi]`, normalising the arguments: endpoints out of order
    /// are swapped, and any NaN yields [`Interval::whole`].
    #[inline]
    pub fn new(lo: f64, hi: f64) -> Interval {
        if lo.is_nan() || hi.is_nan() {
            return Interval::whole();
        }
        if lo <= hi {
            Interval { lo, hi }
        } else {
            Interval { lo: hi, hi: lo }
        }
    }

    /// A degenerate (point) interval `[x, x]` holding the exact `f64` value `x`.
    ///
    /// No widening is applied: `x` is represented exactly. NaN yields
    /// [`Interval::whole`].
    #[inline]
    pub fn point(x: f64) -> Interval {
        if x.is_nan() {
            Interval::whole()
        } else {
            Interval { lo: x, hi: x }
        }
    }

    /// The whole real line `[-∞, +∞]` — the top element / no-information value.
    #[inline]
    pub const fn whole() -> Interval {
        Interval {
            lo: f64::NEG_INFINITY,
            hi: f64::INFINITY,
        }
    }

    /// Lower endpoint.
    #[inline]
    pub const fn lo(self) -> f64 {
        self.lo
    }

    /// Upper endpoint.
    #[inline]
    pub const fn hi(self) -> f64 {
        self.hi
    }

    /// Width `hi - lo`, rounded up so the returned value is never an
    /// under-estimate of the true width.
    #[inline]
    pub fn width(self) -> f64 {
        wide_up(self.hi - self.lo, 1)
    }

    /// Midpoint `(lo + hi) / 2`.
    #[inline]
    pub fn mid(self) -> f64 {
        0.5 * (self.lo + self.hi)
    }

    /// Magnitude: `max(|lo|, |hi|)`.
    #[inline]
    pub fn mag(self) -> f64 {
        let a = mc::abs(self.lo);
        let b = mc::abs(self.hi);
        if a > b {
            a
        } else {
            b
        }
    }

    /// Both endpoints finite?
    #[inline]
    pub fn is_finite(self) -> bool {
        self.lo.is_finite() && self.hi.is_finite()
    }

    /// Does the interval contain the real value `x`?
    #[inline]
    pub fn contains(self, x: f64) -> bool {
        !x.is_nan() && self.lo <= x && x <= self.hi
    }

    /// Does the interval contain all of `other`?
    #[inline]
    pub fn contains_interval(self, other: Interval) -> bool {
        self.lo <= other.lo && other.hi <= self.hi
    }

    /// Smallest interval containing both `self` and `other`.
    #[inline]
    pub fn hull(self, other: Interval) -> Interval {
        Interval::new(self.lo.min(other.lo), self.hi.max(other.hi))
    }

    // ---- reciprocal & power ----

    /// Reciprocal `1 / self` (enclosure).
    ///
    /// If the interval straddles zero the reciprocal is unbounded, so the result
    /// is [`Interval::whole`]. Otherwise `1/x` is monotone on the interval and
    /// the endpoints map to `[1/hi, 1/lo]`, widened outward by 1 ULP.
    #[inline]
    pub fn recip(self) -> Interval {
        if self.lo <= 0.0 && self.hi >= 0.0 {
            return Interval::whole();
        }
        let a = 1.0 / self.hi;
        let b = 1.0 / self.lo;
        let (lo, hi) = if a < b { (a, b) } else { (b, a) };
        Interval::new(wide_down(lo, 1), wide_up(hi, 1))
    }

    /// Integer power `self^n` (enclosure).
    ///
    /// Enclosure argument: for an *odd* exponent `x^p` is strictly increasing,
    /// so the range is `[lo^p, hi^p]`. For an *even* exponent `x^p` is a
    /// symmetric U shape: the maximum is `max(|lo|,|hi|)^p`, and the minimum is
    /// `0` when the interval straddles zero, else `min(|lo|,|hi|)^p`. Negative
    /// exponents compose the natural power with [`Interval::recip`].
    pub fn powi(self, n: i32) -> Interval {
        if n == 0 {
            return Interval::point(1.0);
        }
        let p = n.unsigned_abs();
        let nat = self.pow_nat(p);
        if n < 0 {
            nat.recip()
        } else {
            nat
        }
    }

    /// `self^p` for `p >= 1` (natural exponent). See [`Interval::powi`].
    fn pow_nat(self, p: u32) -> Interval {
        if p == 1 {
            return self;
        }
        let pi = p as i32;
        if p.is_multiple_of(2) {
            let a = mc::abs(self.lo);
            let b = mc::abs(self.hi);
            let (mn, mx) = if a < b { (a, b) } else { (b, a) };
            let lo = if self.lo <= 0.0 && self.hi >= 0.0 {
                0.0
            } else {
                mc::powi(mn, pi)
            };
            let hi = mc::powi(mx, pi);
            Interval::new(wide_down(lo, 1), wide_up(hi, 1))
        } else {
            let lo = mc::powi(self.lo, pi);
            let hi = mc::powi(self.hi, pi);
            Interval::new(wide_down(lo, 1), wide_up(hi, 1))
        }
    }

    // ---- elementary functions ----

    /// Absolute value (exact; no rounding needed).
    ///
    /// Enclosure argument: `|x|` is 0 at the origin and increases outward, so on
    /// a zero-straddling interval the range is `[0, max(|lo|,|hi|)]`; otherwise
    /// it is the reflected endpoints.
    #[inline]
    pub fn abs(self) -> Interval {
        if self.lo >= 0.0 {
            self
        } else if self.hi <= 0.0 {
            Interval::new(-self.hi, -self.lo)
        } else {
            let m = if -self.lo > self.hi {
                -self.lo
            } else {
                self.hi
            };
            Interval::new(0.0, m)
        }
    }

    /// Square root (enclosure over the non-negative part of the domain).
    ///
    /// Enclosure argument: `sqrt` is increasing and correctly rounded, so a
    /// 1-ULP outward margin is rigorous. Negative inputs are outside the domain;
    /// the lower endpoint is clamped to 0 (the boundary of the domain), and a
    /// fully-negative interval — no real root — yields [`Interval::whole`].
    #[inline]
    pub fn sqrt(self) -> Interval {
        if self.hi < 0.0 {
            return Interval::whole();
        }
        let lo_in = if self.lo < 0.0 { 0.0 } else { self.lo };
        let lo = wide_down(mc::sqrt(lo_in), 1);
        let hi = wide_up(mc::sqrt(self.hi), 1);
        Interval::new(lo.max(0.0), hi)
    }

    /// `hypot(self, other) = sqrt(self² + other²)` (enclosure).
    ///
    /// Enclosure argument: `hypot` is increasing in each of `|self|` and
    /// `|other|`, so the extremes are attained at the extremal magnitudes.
    #[inline]
    pub fn hypot(self, other: Interval) -> Interval {
        let a = self.abs();
        let b = other.abs();
        let lo = wide_down(mc::hypot(a.lo, b.lo), TRANSCEND_ULP);
        let hi = wide_up(mc::hypot(a.hi, b.hi), TRANSCEND_ULP);
        Interval::new(lo.max(0.0), hi)
    }

    /// Natural exponential (enclosure).
    ///
    /// Enclosure argument: `exp` is strictly increasing on all of ℝ, so the
    /// range is `[exp(lo), exp(hi)]`, widened by [`TRANSCEND_ULP`].
    #[inline]
    pub fn exp(self) -> Interval {
        let lo = wide_down(mc::exp(self.lo), TRANSCEND_ULP);
        let hi = wide_up(mc::exp(self.hi), TRANSCEND_ULP);
        Interval::new(lo.max(0.0), hi)
    }

    /// Natural logarithm (enclosure over the positive part of the domain).
    ///
    /// Enclosure argument: `ln` is strictly increasing on `(0, ∞)`. If the
    /// interval reaches down to or below 0 the lower bound is `-∞` (the limit as
    /// `x → 0⁺`); a fully non-positive interval has no real logarithm and yields
    /// [`Interval::whole`].
    #[inline]
    pub fn ln(self) -> Interval {
        if self.hi <= 0.0 {
            return Interval::whole();
        }
        let lo = if self.lo <= 0.0 {
            f64::NEG_INFINITY
        } else {
            wide_down(mc::ln(self.lo), TRANSCEND_ULP)
        };
        let hi = wide_up(mc::ln(self.hi), TRANSCEND_ULP);
        Interval::new(lo, hi)
    }

    /// Inverse hyperbolic sine `asinh(x) = ln(x + sqrt(x² + 1))` (enclosure).
    ///
    /// Enclosure follows by rigorous composition of `powi`, `sqrt`, `add`, and
    /// `ln`, each of which is itself an enclosure. `asinh` is defined on all of
    /// ℝ so no domain handling is required.
    #[inline]
    pub fn asinh(self) -> Interval {
        let s = (self.powi(2) + Interval::point(1.0)).sqrt();
        (self + s).ln()
    }

    /// Inverse hyperbolic tangent `atanh(x) = ½·ln((1+x)/(1-x))` (enclosure).
    ///
    /// Valid for `|x| < 1`; if the interval reaches `±1` the internal division
    /// straddles zero and the result widens toward the unbounded limit.
    #[inline]
    pub fn atanh(self) -> Interval {
        let num = Interval::point(1.0) + self;
        let den = Interval::point(1.0) - self;
        Interval::point(0.5) * (num / den).ln()
    }

    /// Sine (enclosure with full quadrant analysis).
    ///
    /// Enclosure argument: `sin` attains a local maximum `+1` only at
    /// `π/2 + 2kπ` and a local minimum `-1` only at `-π/2 + 2kπ`. If the
    /// interval spans no full period (`width < 2π`) and contains none of those
    /// critical points, `sin` is monotone across it and the extremes are the
    /// endpoint values. When a critical point *is* enclosed we widen that side
    /// to the exact extremum `±1`. Intervals of width `≥ 2π`, or with an
    /// infinite endpoint, collapse to `[-1, 1]`.
    pub fn sin(self) -> Interval {
        if !self.is_finite() || self.width() >= TWO_PI {
            return Interval::new(-1.0, 1.0);
        }
        let sa = mc::sin(self.lo);
        let sb = mc::sin(self.hi);
        let mut lo = sa.min(sb);
        let mut hi = sa.max(sb);
        if contains_grid_point(self.lo, self.hi, HALF_PI, TWO_PI) {
            hi = 1.0;
        }
        if contains_grid_point(self.lo, self.hi, -HALF_PI, TWO_PI) {
            lo = -1.0;
        }
        let lo = wide_down(lo, TRANSCEND_ULP).max(-1.0);
        let hi = wide_up(hi, TRANSCEND_ULP).min(1.0);
        Interval::new(lo, hi)
    }

    /// Cosine (enclosure with full quadrant analysis).
    ///
    /// Enclosure argument: `cos` peaks at `+1` on `2kπ` and troughs at `-1` on
    /// `π + 2kπ`. Same monotone-between-critical-points reasoning as
    /// [`Interval::sin`].
    pub fn cos(self) -> Interval {
        if !self.is_finite() || self.width() >= TWO_PI {
            return Interval::new(-1.0, 1.0);
        }
        let ca = mc::cos(self.lo);
        let cb = mc::cos(self.hi);
        let mut lo = ca.min(cb);
        let mut hi = ca.max(cb);
        if contains_grid_point(self.lo, self.hi, 0.0, TWO_PI) {
            hi = 1.0;
        }
        if contains_grid_point(self.lo, self.hi, PI, TWO_PI) {
            lo = -1.0;
        }
        let lo = wide_down(lo, TRANSCEND_ULP).max(-1.0);
        let hi = wide_up(hi, TRANSCEND_ULP).min(1.0);
        Interval::new(lo, hi)
    }

    /// Tangent (enclosure; unbounded across a pole).
    ///
    /// Enclosure argument: `tan` has poles at `π/2 + kπ` and is strictly
    /// increasing on every pole-free branch. If the interval could enclose a
    /// pole (its width reaches `π`, an endpoint is infinite, or a pole location
    /// lies inside) the tangent is unbounded and the result is
    /// [`Interval::whole`]. Otherwise the endpoints give `[tan(lo), tan(hi)]`.
    pub fn tan(self) -> Interval {
        if !self.is_finite() || self.width() >= PI {
            return Interval::whole();
        }
        if contains_grid_point(self.lo, self.hi, HALF_PI, PI) {
            return Interval::whole();
        }
        let lo = wide_down(mc::tan(self.lo), TRANSCEND_ULP);
        let hi = wide_up(mc::tan(self.hi), TRANSCEND_ULP);
        Interval::new(lo, hi)
    }

    /// Arctangent (enclosure).
    ///
    /// Enclosure argument: `atan` is strictly increasing and bounded in
    /// `(-π/2, π/2)`, so the range is simply the widened endpoint images.
    #[inline]
    pub fn atan(self) -> Interval {
        let lo = wide_down(mc::atan(self.lo), TRANSCEND_ULP);
        let hi = wide_up(mc::atan(self.hi), TRANSCEND_ULP);
        Interval::new(lo, hi)
    }

    /// Arcsine (enclosure over the clamped domain `[-1, 1]`).
    ///
    /// Enclosure argument: `asin` is strictly increasing on `[-1, 1]`. The input
    /// is intersected with the domain before evaluation; an interval disjoint
    /// from `[-1, 1]` has no real arcsine and yields the full range
    /// `[-π/2, π/2]`.
    pub fn asin(self) -> Interval {
        let a = self.lo.max(-1.0);
        let b = self.hi.min(1.0);
        if a > b {
            return Interval::new(-HALF_PI, HALF_PI);
        }
        let lo = wide_down(mc::asin(a), TRANSCEND_ULP).max(-HALF_PI);
        let hi = wide_up(mc::asin(b), TRANSCEND_ULP).min(HALF_PI);
        Interval::new(lo, hi)
    }

    /// Arccosine (enclosure over the clamped domain `[-1, 1]`).
    ///
    /// Enclosure argument: `acos` is strictly *decreasing* on `[-1, 1]`, so the
    /// endpoints map to `[acos(hi), acos(lo)]`. Domain handling matches
    /// [`Interval::asin`].
    pub fn acos(self) -> Interval {
        let a = self.lo.max(-1.0);
        let b = self.hi.min(1.0);
        if a > b {
            return Interval::new(0.0, PI);
        }
        let lo = wide_down(mc::acos(b), TRANSCEND_ULP).max(0.0);
        let hi = wide_up(mc::acos(a), TRANSCEND_ULP).min(PI);
        Interval::new(lo, hi)
    }

    /// Two-argument arctangent `atan2(self, other)` where `self` is the `y`
    /// range and `other` the `x` range (enclosure).
    ///
    /// Enclosure argument: `atan2` is smooth except at the origin and across the
    /// branch cut on the negative `x`-axis. We special-case those:
    /// - the box straddling the **origin** (both ranges contain 0) or the
    ///   **branch cut** (`x < 0` with `y` spanning 0) can reach every angle, so
    ///   the result is `[-π, π]`;
    /// - otherwise, splitting the box at the coordinate axes yields sub-boxes
    ///   confined to a single quadrant, where `∂atan2/∂x` and `∂atan2/∂y` keep
    ///   constant sign, so each sub-box's extremes lie at its corners. The hull
    ///   of all corner evaluations is a rigorous enclosure.
    pub fn atan2(self, other: Interval) -> Interval {
        let y = self;
        let x = other;
        let y_has_zero = y.lo <= 0.0 && y.hi >= 0.0;
        let x_has_zero = x.lo <= 0.0 && x.hi >= 0.0;
        // Origin inside, or negative-x branch cut crossed: every angle reachable.
        if y_has_zero && (x_has_zero || x.lo < 0.0) {
            return Interval::new(wide_down(-PI, 1), wide_up(PI, 1));
        }
        // Split at the axes so each sub-box stays in one quadrant.
        let mut xs = [x.lo, x.hi, 0.0];
        let mut ys = [y.lo, y.hi, 0.0];
        let nx = if x.lo < 0.0 && x.hi > 0.0 { 3 } else { 2 };
        let ny = if y.lo < 0.0 && y.hi > 0.0 { 3 } else { 2 };
        // Order the split coordinates ascending within the used prefix.
        sort_prefix(&mut xs, nx);
        sort_prefix(&mut ys, ny);
        let mut lo = f64::INFINITY;
        let mut hi = f64::NEG_INFINITY;
        let mut i = 0;
        while i + 1 < nx {
            let mut j = 0;
            while j + 1 < ny {
                for &yy in &[ys[j], ys[j + 1]] {
                    for &xx in &[xs[i], xs[i + 1]] {
                        let a = mc::atan2(yy, xx);
                        if a < lo {
                            lo = a;
                        }
                        if a > hi {
                            hi = a;
                        }
                    }
                }
                j += 1;
            }
            i += 1;
        }
        Interval::new(wide_down(lo, TRANSCEND_ULP), wide_up(hi, TRANSCEND_ULP))
    }
}

/// Ascending insertion-sort of the first `n` elements of `arr` (n ≤ 3).
#[inline]
fn sort_prefix(arr: &mut [f64; 3], n: usize) {
    let mut i = 1;
    while i < n {
        let mut j = i;
        while j > 0 && arr[j - 1] > arr[j] {
            arr.swap(j - 1, j);
            j -= 1;
        }
        i += 1;
    }
}

impl From<f64> for Interval {
    #[inline]
    fn from(x: f64) -> Interval {
        Interval::point(x)
    }
}

impl Add for Interval {
    type Output = Interval;
    /// Sum enclosure: `[lo₁+lo₂, hi₁+hi₂]`, each endpoint widened 1 ULP.
    #[inline]
    fn add(self, o: Interval) -> Interval {
        Interval::new(wide_down(self.lo + o.lo, 1), wide_up(self.hi + o.hi, 1))
    }
}

impl Sub for Interval {
    type Output = Interval;
    /// Difference enclosure: `[lo₁-hi₂, hi₁-lo₂]`, each endpoint widened 1 ULP.
    #[inline]
    fn sub(self, o: Interval) -> Interval {
        Interval::new(wide_down(self.lo - o.hi, 1), wide_up(self.hi - o.lo, 1))
    }
}

impl Mul for Interval {
    type Output = Interval;
    /// Product enclosure: hull of the four endpoint products, widened 1 ULP.
    #[inline]
    fn mul(self, o: Interval) -> Interval {
        let (lo, hi) = hull4(
            self.lo * o.lo,
            self.lo * o.hi,
            self.hi * o.lo,
            self.hi * o.hi,
        );
        Interval::new(wide_down(lo, 1), wide_up(hi, 1))
    }
}

impl Div for Interval {
    type Output = Interval;
    /// Quotient enclosure. If the divisor straddles zero the quotient is
    /// unbounded and the result is [`Interval::whole`]; otherwise it is the hull
    /// of the four endpoint quotients, widened 1 ULP.
    #[inline]
    fn div(self, o: Interval) -> Interval {
        if o.lo <= 0.0 && o.hi >= 0.0 {
            return Interval::whole();
        }
        let (lo, hi) = hull4(
            self.lo / o.lo,
            self.lo / o.hi,
            self.hi / o.lo,
            self.hi / o.hi,
        );
        Interval::new(wide_down(lo, 1), wide_up(hi, 1))
    }
}

impl Neg for Interval {
    type Output = Interval;
    /// Negation is exact: `[-hi, -lo]`.
    #[inline]
    fn neg(self) -> Interval {
        Interval {
            lo: -self.hi,
            hi: -self.lo,
        }
    }
}

impl fmt::Debug for Interval {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "[{}, {}]", self.lo, self.hi)
    }
}

impl fmt::Display for Interval {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "[{}, {}]", self.lo, self.hi)
    }
}

/// A 2-D coordinate whose components are guaranteed-enclosure [`Interval`]s.
///
/// Returned by the interval forward transforms (merc, tmerc). Each component
/// provably contains the exact projected value.
#[derive(Clone, Copy, PartialEq, Debug)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct IntervalCoord {
    /// Easting / `x` enclosure.
    pub x: Interval,
    /// Northing / `y` enclosure.
    pub y: Interval,
}

impl IntervalCoord {
    /// Construct from two component intervals.
    #[inline]
    pub const fn new(x: Interval, y: Interval) -> IntervalCoord {
        IntervalCoord { x, y }
    }

    /// Does this box contain the point `(x, y)`?
    #[inline]
    pub fn contains(self, x: f64, y: f64) -> bool {
        self.x.contains(x) && self.y.contains(y)
    }
}

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

    /// Reference `next_up`/`next_down`-based enclosure sanity: an exact value is
    /// always inside its own point interval.
    #[test]
    fn point_contains_itself() {
        for &v in &[0.0, -1.5, 3.25e12, -7e-8] {
            assert!(Interval::point(v).contains(v));
        }
    }

    #[test]
    fn add_sub_mul_div_enclose_exact() {
        // Random-ish deterministic sweep: (a op b) must be inside the interval op.
        let vals = [-3.7, -1.0, -0.25, 0.0, 0.5, 2.0, 9.9, 123.456];
        for &a in &vals {
            for &b in &vals {
                let ia = Interval::point(a);
                let ib = Interval::point(b);
                assert!((ia + ib).contains(a + b), "add {a} {b}");
                assert!((ia - ib).contains(a - b), "sub {a} {b}");
                assert!((ia * ib).contains(a * b), "mul {a} {b}");
                if b != 0.0 {
                    assert!((ia / ib).contains(a / b), "div {a} {b}");
                }
            }
        }
    }

    #[test]
    fn known_rational_enclosures_are_tight() {
        // Exactly representable operations: enclosure width stays at the 1-ULP
        // rounding margin (a few ULPs), and the exact result is inside.
        let two = Interval::point(2.0);
        let three = Interval::point(3.0);
        let sum = two + three; // exactly 5
        assert!(sum.contains(5.0));
        assert!(sum.width() < 1e-14);
        let prod = two * three; // exactly 6
        assert!(prod.contains(6.0));
        assert!(prod.width() < 1e-14);
        // [2,3]^2 = [4,9] tight (no zero crossing).
        let sq = Interval::new(2.0, 3.0).powi(2);
        assert!(sq.contains(4.0) && sq.contains(9.0));
        assert!(sq.lo() <= 4.0 && sq.hi() >= 9.0);
        assert!(sq.lo() > 3.99 && sq.hi() < 9.01);
        // Even power straddling zero has minimum exactly 0.
        let sqz = Interval::new(-2.0, 3.0).powi(2);
        assert!(sqz.contains(0.0));
        assert!(sqz.lo() <= 0.0 && sqz.hi() >= 9.0);
    }

    #[test]
    fn monotone_widening_under_input_widening() {
        // Widening an input interval must not shrink the output enclosure.
        let narrow = Interval::new(0.30, 0.31);
        let wide = Interval::new(0.20, 0.41);
        assert!(wide.contains_interval(narrow));
        for f in [
            |i: Interval| i.sin(),
            |i: Interval| i.cos(),
            |i: Interval| i.exp(),
            |i: Interval| i.atan(),
            |i: Interval| i * i,
            |i: Interval| i.sqrt(),
        ] {
            let fn_narrow = f(narrow);
            let fn_wide = f(wide);
            assert!(
                fn_wide.contains_interval(fn_narrow),
                "widening shrank output: narrow={fn_narrow:?} wide={fn_wide:?}"
            );
        }
    }

    #[test]
    fn transcendentals_enclose_f64() {
        // Point interval evaluated through each transcendental must contain the
        // corresponding std f64 result.
        let xs = [-2.3, -0.7, 0.0, 0.4, 1.1, 2.9, 5.5];
        for &x in &xs {
            assert!(Interval::point(x).sin().contains(x.sin()), "sin {x}");
            assert!(Interval::point(x).cos().contains(x.cos()), "cos {x}");
            assert!(Interval::point(x).atan().contains(x.atan()), "atan {x}");
            assert!(Interval::point(x).exp().contains(x.exp()), "exp {x}");
            if x > 0.0 {
                assert!(Interval::point(x).ln().contains(x.ln()), "ln {x}");
                assert!(Interval::point(x).sqrt().contains(x.sqrt()), "sqrt {x}");
            }
            if x.abs() <= 1.0 {
                assert!(Interval::point(x).asin().contains(x.asin()), "asin {x}");
                assert!(Interval::point(x).acos().contains(x.acos()), "acos {x}");
            }
        }
    }

    #[test]
    fn sin_cos_critical_points_reach_extrema() {
        use core::f64::consts::PI;
        // An interval that spans π/2 encloses the sine maximum +1.
        let s = Interval::new(PI / 2.0 - 0.1, PI / 2.0 + 0.1).sin();
        assert!(s.hi() >= 1.0 - 1e-12, "sin max not reached: {s:?}");
        // An interval spanning -π/2 encloses the sine minimum -1.
        let s2 = Interval::new(-PI / 2.0 - 0.1, -PI / 2.0 + 0.1).sin();
        assert!(s2.lo() <= -1.0 + 1e-12, "sin min not reached: {s2:?}");
        // cos maximum at 0, minimum at π.
        let c = Interval::new(-0.2, 0.2).cos();
        assert!(c.hi() >= 1.0 - 1e-12, "cos max not reached: {c:?}");
        let c2 = Interval::new(PI - 0.2, PI + 0.2).cos();
        assert!(c2.lo() <= -1.0 + 1e-12, "cos min not reached: {c2:?}");
        // Wide interval collapses to [-1, 1].
        let full = Interval::new(-10.0, 10.0).sin();
        assert!(full.lo() <= -1.0 + 1e-12 && full.hi() >= 1.0 - 1e-12);
    }

    #[test]
    fn sin_encloses_true_extrema_densely() {
        // Over a set of sub-intervals, the enclosure must contain the true sin
        // value at every interior sample.
        let a = 0.3;
        let b = 2.9; // spans π/2
        let iv = Interval::new(a, b).sin();
        let mut t = a;
        while t <= b {
            assert!(iv.contains(t.sin()), "sin({t}) escaped {iv:?}");
            t += 0.01;
        }
    }

    #[test]
    fn tan_pole_gives_whole_line() {
        use core::f64::consts::PI;
        // Interval straddling π/2 encloses a pole → unbounded.
        let t = Interval::new(PI / 2.0 - 0.05, PI / 2.0 + 0.05).tan();
        assert!(!t.is_finite());
        // Pole-free branch stays finite and encloses the f64 value.
        let t2 = Interval::new(0.1, 0.2).tan();
        assert!(t2.is_finite());
        assert!(t2.contains(0.15_f64.tan()));
    }

    #[test]
    fn atan2_matches_quadrants_and_encloses() {
        // First-quadrant box: enclosure contains all sampled corner angles.
        let y = Interval::new(1.0, 2.0);
        let x = Interval::new(3.0, 4.0);
        let a = y.atan2(x);
        for &yy in &[1.0_f64, 1.5, 2.0] {
            for &xx in &[3.0_f64, 3.5, 4.0] {
                assert!(a.contains(yy.atan2(xx)), "atan2 {yy} {xx} escaped {a:?}");
            }
        }
        // Origin-containing box spans the full angular range.
        let full = Interval::new(-1.0, 1.0).atan2(Interval::new(-1.0, 1.0));
        assert!(full.lo() <= -PI + 1e-9 && full.hi() >= PI - 1e-9);
    }

    #[test]
    fn recip_and_negative_powi() {
        let r = Interval::new(2.0, 4.0).recip();
        assert!(r.contains(0.25) && r.contains(0.5));
        let p = Interval::new(2.0, 2.0).powi(-2);
        assert!(p.contains(0.25));
        // Reciprocal across zero is unbounded.
        assert!(!Interval::new(-1.0, 1.0).recip().is_finite());
    }

    #[test]
    fn asinh_atanh_compose_correctly() {
        for &x in &[-0.9_f64, -0.3, 0.0, 0.2, 0.8] {
            let expect = (
                (x + (x * x + 1.0).sqrt()).ln(),
                0.5 * ((1.0 + x) / (1.0 - x)).ln(),
            );
            assert!(Interval::point(x).asinh().contains(expect.0), "asinh {x}");
            assert!(Interval::point(x).atanh().contains(expect.1), "atanh {x}");
        }
    }

    #[test]
    fn hull_and_intersect_helpers() {
        let a = Interval::new(0.0, 2.0);
        let b = Interval::new(1.0, 3.0);
        let h = a.hull(b);
        assert_eq!((h.lo(), h.hi()), (0.0, 3.0));
        assert!(a.contains_interval(Interval::new(0.5, 1.5)));
    }
}