rectgrid 0.4.1

use alloc::{rc::Rc, vec::Vec};
use core::{
    array::from_fn,
    cmp::Ordering,
    fmt::{self, Debug, Formatter},
    marker::PhantomData,
    ops::{Add, AddAssign, Deref, Div, Mul, Neg, Sub, SubAssign},
    primitive::{f64, u32, usize},
    result::Result,
};

use crate::Error;

/// An `f64` tagged with its unit system at the type level. The tag prevents mixing values
/// of different unit systems; see [`Px`], [`Unit`] and [`Parameter`].
#[repr(transparent)]
pub struct Value<Tag>(f64, PhantomData<Tag>);

impl<Tag> Clone for Value<Tag> {
    fn clone(&self) -> Self {
        *self
    }
}

impl<Tag> Copy for Value<Tag> {}

impl<Tag> Value<Tag> {
    /// Wraps a raw `f64`.
    pub fn new(v: f64) -> Self {
        Self(v, PhantomData)
    }

    /// Returns the raw `f64`.
    pub fn get(self) -> f64 {
        self.0
    }
}

impl<Tag> Debug for Value<Tag> {
    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
        f.debug_tuple("Value").field(&self.0).finish()
    }
}

impl<Tag> PartialEq for Value<Tag> {
    fn eq(&self, other: &Self) -> bool {
        self.0 == other.0
    }
}

impl<Tag> PartialOrd for Value<Tag> {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        self.0.partial_cmp(&other.0)
    }
}

impl<Tag> Default for Value<Tag> {
    fn default() -> Self {
        Self::new(0.0)
    }
}

impl<Tag> Neg for Value<Tag> {
    type Output = Value<Tag>;
    fn neg(self) -> Self::Output {
        Value::new(-self.0)
    }
}

impl<Tag> SubAssign for Value<Tag> {
    fn sub_assign(&mut self, rhs: Self) {
        self.0 -= rhs.0;
    }
}

impl<Tag> Mul<Value<Tag>> for f64 {
    type Output = Value<Tag>;
    fn mul(self, rhs: Value<Tag>) -> Self::Output {
        Value::new(self * rhs.0)
    }
}

impl<Tag> Div<f64> for Value<Tag> {
    type Output = Value<Tag>;
    fn div(self, rhs: f64) -> Self::Output {
        Value::new(self.0 / rhs)
    }
}

impl<Tag> From<f64> for Value<Tag> {
    fn from(v: f64) -> Self {
        Self::new(v)
    }
}

impl<Tag> Add for Value<Tag> {
    type Output = Value<Tag>;
    fn add(self, rhs: Self) -> Self::Output {
        Value::new(self.0 + rhs.0)
    }
}

impl<Tag> AddAssign for Value<Tag> {
    fn add_assign(&mut self, rhs: Self) {
        self.0 += rhs.0;
    }
}

impl<Tag> Sub for Value<Tag> {
    type Output = Value<Tag>;
    fn sub(self, rhs: Self) -> Self::Output {
        Value::new(self.0 - rhs.0)
    }
}

impl<Tag> Mul<f64> for Value<Tag> {
    type Output = Value<Tag>;
    fn mul(self, rhs: f64) -> Self::Output {
        Value::new(self.0 * rhs)
    }
}

impl<Tag> Div for Value<Tag> {
    type Output = f64;
    fn div(self, rhs: Self) -> f64 {
        self.0 / rhs.0
    }
}

/// Tag of [`Px`].
pub struct PxTag;
/// Tag of [`Unit`].
pub struct UnitTag;
/// Tag of [`Parameter`].
pub struct ParameterTag;

/// Orthogonal, axis-independent unit system extending infinitely in the positive direction from the origin.
pub type Px = Value<PxTag>;

/// Arbitrary orthogonal unit system, finite or infinite per axis, whose value is an ordinal from its own origin in the positive direction.
pub type Unit = Value<UnitTag>;

/// Unbounded local coordinate system for a single boundary box, where each side length is 1 and sign follows the unit coordinate.
pub type Parameter = Value<ParameterTag>;

/// A point in unit coordinates, one [`Unit`] per axis.
pub type Point<const D: usize> = [Unit; D];

/// Boundary box between `base` and `base + offset`, in unit coordinates.
///
/// Layout, one f64 per cell:
///
/// ```text
/// | base[0] .. base[D-1] | offset[0] .. offset[D-1] |
/// ```
#[derive(Clone, Copy)]
#[repr(C)]
pub struct BBox<const D: usize> {
    pub(crate) base:   Point<D>,
    pub(crate) offset: Point<D>,
}

impl<const D: usize> BBox<D> {
    /// offset is the side length from base toward the far corner. A negative axis is normalized by
    /// advancing base to the far corner and negating offset, so offset is always non-negative.
    ///
    /// ```
    /// use rectgrid::{BBox, Unit};
    /// let bx = BBox::new([Unit::new(0.0); 2], [Unit::new(1.0), Unit::new(3.0)]);
    /// assert!(bx.has_size());
    /// let flipped = BBox::new([Unit::new(0.0); 2], [Unit::new(-1.0), Unit::new(3.0)]);
    /// assert_eq!(flipped.base()[0].get(), -1.0);
    /// assert_eq!(flipped.offset()[0].get(), 1.0);
    /// ```
    pub fn new(base: Point<D>, offset: Point<D>) -> Self {
        let mut base = base;
        let mut offset = offset;
        for d in 0..D {
            if offset[d].get() < 0.0 {
                base[d] += offset[d];
                offset[d] = -offset[d];
            }
        }
        Self { base, offset }
    }

    /// Start point of the boundary box.
    pub fn base(&self) -> Point<D> {
        self.base
    }

    /// Vector from `base` to the far corner. Every axis is non-negative; see [`BBox::new`].
    pub fn offset(&self) -> Point<D> {
        self.offset
    }

    /// Floors base/offset to an integer grid. extend applies to base only, added before flooring.
    ///
    /// ```
    /// use rectgrid::{BBox, Unit};
    /// let mut bx = BBox::new(
    ///     [Unit::new(2.6), Unit::new(0.6)],
    ///     [Unit::new(1.9), Unit::new(3.0)],
    /// );
    /// bx.snap_floor(Some([Unit::new(-0.5), Unit::new(-0.5)]));
    /// assert_eq!(bx.base()[0].get(), 2.0);
    /// assert_eq!(bx.base()[1].get(), 0.0);
    /// assert_eq!(bx.offset()[0].get(), 1.0);
    /// assert_eq!(bx.offset()[1].get(), 3.0);
    /// ```
    pub fn snap_floor(&mut self, extend: Option<[Unit; D]>) -> &mut Self {
        for (d, u) in self.base.iter_mut().enumerate() {
            let v = if let Some(ext) = extend { *u + ext[d] } else { *u };
            *u = Unit::new(libm::floor(v.get()));
        }
        for u in self.offset.iter_mut() {
            *u = Unit::new(libm::floor(u.get()));
        }
        self
    }

    /// Whether every axis of offset is nonzero, i.e. the boundary box has area/volume.
    ///
    /// ```
    /// use rectgrid::{BBox, Unit};
    /// let area = BBox::new([Unit::new(0.0); 2], [Unit::new(1.0), Unit::new(3.0)]);
    /// assert!(area.has_size());
    /// let segment = BBox::new([Unit::new(0.0); 2], [Unit::new(0.0), Unit::new(3.0)]);
    /// assert!(!segment.has_size());
    /// let point = BBox::new([Unit::new(0.0); 2], [Unit::new(0.0); 2]);
    /// assert!(!point.has_size());
    /// ```
    pub fn has_size(&self) -> bool {
        self.offset.iter().all(|u| u.get() != 0.0)
    }
}

/// A pointer to a step closure `Fn(u32) -> Result<Px, Error>`; implemented for every such `Deref`.
/// The caller picks the pointer: `Rc`, `Arc<dyn Fn(..) + Send + Sync>`, `Box`, `&F`.
pub trait StepFn {
    /// Returns step `i`: the px length of the cell between unit `i` and unit `i + 1`.
    fn step(&self, i: u32) -> Result<Px, Error>;
}

impl<P> StepFn for P
where
    P: Deref,
    P::Target: Fn(u32) -> Result<Px, Error>,
{
    fn step(&self, i: u32) -> Result<Px, Error> {
        (**self)(i)
    }
}

/// Pointer type used when a grid does not name one.
pub type DefaultSteps = Rc<dyn Fn(u32) -> Result<Px, Error>>;

/// Definition of one axis of a grid: where its unit coordinates lie in px.
/// Built into an [`Accumulator`] by [`IncrementFunction::accumulate`].
/// `P` is the pointer to the step closure of `ForwardDifference` (see [`StepFn`]);
/// it defaults to [`DefaultSteps`].
pub enum IncrementFunction<P = DefaultSteps> {
    /// Fn(i) = points\[i+1\] - points\[i\]
    /// OutOfIndex means boundary; the argument is the difference's index (integer).
    /// The closure must, when out of range, saturate to the last valid index within range and return it as OutOfIndex.
    /// The domain ends at the first Err: units beyond it are unevaluable even if f returns Ok again later,
    /// and the reported index must be that of the first Err.
    ForwardDifference(P),
    /// Finite domain. Array enumerating unit values of the interval from the origin in the positive direction.
    VectorList(Vec<Px>),
    /// Infinite domain. NaN is invalid.
    Scale(f64),
}

impl<P: StepFn> IncrementFunction<P> {
    /// Builds forward/inverse accumulators from the definition.
    /// forward converts fractional units to px by linear interpolation.
    /// Returns InvalidDefinition if the definition cannot be evaluated, e.g. an empty VectorList.
    ///
    /// ```
    /// extern crate alloc;
    /// use rectgrid::{Accumulator, IncrementFunction, Px};
    ///
    /// let acc: Accumulator = IncrementFunction::Scale(10.0).accumulate().unwrap();
    /// assert_eq!(acc.forward(2.5).unwrap().get(), 25.0);
    ///
    /// let acc: Accumulator = IncrementFunction::VectorList(alloc::vec![Px::new(0.0), Px::new(10.0), Px::new(30.0)]).accumulate().unwrap();
    /// assert_eq!(acc.forward(1.0).unwrap().get(), 10.0);
    /// assert_eq!(acc.forward(1.5).unwrap().get(), 20.0);
    ///
    /// use rectgrid::Error;
    /// let empty: Result<Accumulator, _> = IncrementFunction::VectorList(alloc::vec![]).accumulate();
    /// assert!(matches!(empty, Err(Error::InvalidDefinition)));
    /// ```
    pub fn accumulate(self) -> Result<Accumulator<P>, Error> {
        match self {
            Self::Scale(s) => {
                if s.is_nan() {
                    return Err(Error::InvalidInput);
                }
                Ok(Accumulator::Scale(s))
            }
            Self::VectorList(pxs) => {
                if pxs.is_empty() {
                    return Err(Error::InvalidDefinition);
                }
                Ok(Accumulator::VectorList(pxs))
            }
            Self::ForwardDifference(f) => Ok(Accumulator::ForwardDifference(f)),
        }
    }
}

/// Sums the first `floor(x)` steps and interpolates linearly through the next one.
fn forward_difference_forward<P: StepFn>(f: &P, x: f64) -> Result<Px, Error> {
    let n = libm::floor(x) as u32;
    let frac = x - n as f64;
    let mut acc = Px::new(0.0);
    for k in 0..n {
        acc += f.step(k)?;
    }
    if frac != 0.0 {
        acc += f.step(n)? * frac;
    }
    Ok(acc)
}

/// Inverts px to unit for a ForwardDifference definition.
/// forward is piecewise linear (`S(k) + f(k) * (x - k)` on [k, k+1], S(k) = f(0) + ... + f(k-1)), so the
/// inverse scans the segments until S(k) <= target <= S(k+1) and solves the last one in closed form.
/// The sums are accumulated in the same order as forward, so `forward(k)` maps back to exactly `k`.
/// Cost is one call of f per segment up to the answer; an unbounded f whose sum never reaches target
/// is scanned up to u32::MAX segments.
/// Caller contract: f must be monotonically non-decreasing over Unit >= 0, matching IncrementFunction::ForwardDifference.
fn forward_difference_inverse<P: StepFn>(f: &P, target: Px) -> Result<Unit, Error> {
    let target = target.get();
    if !target.is_finite() {
        return Err(Error::InvalidInput);
    }
    if target == 0.0 {
        return Ok(Unit::new(0.0));
    }
    if target < 0.0 {
        let step = f.step(0)?.get();
        if step == 0.0 {
            return Err(Error::OutOfIndex(0));
        }
        return Ok(Unit::new(target / step));
    }

    let mut accumulated = 0.0;
    for k in 0..=u32::MAX {
        let step = f.step(k)?.get();
        let next = accumulated + step;
        if next >= target {
            // next == target is forward(k + 1) exactly
            let frac = if next == target {
                1.0
            } else if step == 0.0 {
                0.0
            } else {
                ((target - accumulated) / step).min(1.0)
            };
            return Ok(Unit::new(k as f64 + frac));
        }
        accumulated = next;
    }
    Err(Error::OutOfIndex(u32::MAX))
}

/// Forward/inverse coordinate conversion for a single axis, built from an IncrementFunction.
pub enum Accumulator<P = DefaultSteps> {
    /// Built from [`IncrementFunction::Scale`]. The inverse is solved directly as `target / s`.
    Scale(f64),
    /// Built from [`IncrementFunction::VectorList`]. The inverse is solved by `partition_point`
    /// and a linear interpolation within the located segment.
    VectorList(Vec<Px>),
    /// Built from [`IncrementFunction::ForwardDifference`]. The inverse scans the segments.
    ForwardDifference(P),
}

impl<P: StepFn> Accumulator<P> {
    /// unit coordinate -> px coordinate. `x` is measured from the origin of the axis; its
    /// fractional part is interpolated linearly within the cell. Returns `OutOfIndex` past a
    /// finite domain and `InvalidInput` for NaN.
    ///
    /// ```
    /// use rectgrid::{Accumulator, IncrementFunction};
    /// let acc: Accumulator = IncrementFunction::Scale(10.0).accumulate().unwrap();
    /// assert_eq!(acc.forward(2.5).unwrap().get(), 25.0);
    /// ```
    pub fn forward(&self, x: f64) -> Result<Px, Error> {
        if x.is_nan() {
            return Err(Error::InvalidInput);
        }
        match self {
            Self::Scale(s) => Ok(Px::new(s * x)),
            Self::VectorList(pxs) => {
                let last = (pxs.len() - 1) as u32;
                if x < 0.0 {
                    return Err(Error::OutOfIndex(0));
                }
                let n = libm::floor(x) as usize;
                let frac = x - n as f64;
                let lo = *pxs.get(n).ok_or(Error::OutOfIndex(last))?;
                if frac == 0.0 {
                    return Ok(lo);
                }
                let hi = *pxs.get(n + 1).ok_or(Error::OutOfIndex(last))?;
                Ok(lo + (hi - lo) * frac)
            }
            Self::ForwardDifference(f) => forward_difference_forward(f, x),
        }
    }

    /// px coordinate -> unit coordinate. For every variant, `inverse(forward(k))` is exactly the integer `k`.
    ///
    /// ```
    /// use rectgrid::{Accumulator, IncrementFunction, Px};
    /// let acc: Accumulator = IncrementFunction::Scale(10.0).accumulate().unwrap();
    /// assert_eq!(acc.inverse(Px::new(25.0)).unwrap().get(), 2.5);
    /// ```
    pub fn inverse(&self, target: Px) -> Result<Unit, Error> {
        if target.get().is_nan() {
            return Err(Error::InvalidInput);
        }
        match self {
            Self::Scale(s) => {
                // the quotient can miss an integer by an ulp (0.3 / 0.1 = 2.9999999999999996)
                let quotient = target.get() / s;
                let nearest = libm::round(quotient);
                Ok(Unit::new(if s * nearest == target.get() { nearest } else { quotient }))
            }
            Self::VectorList(pxs) => vector_list_inverse(pxs, target),
            Self::ForwardDifference(f) => forward_difference_inverse(f, target),
        }
    }
}

/// Inverse for VectorList: `partition_point` over the (non-empty, non-decreasing) array, then a linear solve in the located segment.
fn vector_list_inverse(pxs: &[Px], target: Px) -> Result<Unit, Error> {
    let last = (pxs.len() - 1) as u32;
    let t = target.get();
    if t < pxs[0].get() {
        return Err(Error::OutOfIndex(0));
    }
    if t > pxs[pxs.len() - 1].get() {
        return Err(Error::OutOfIndex(last));
    }
    let n = pxs.partition_point(|p| p.get() < t).saturating_sub(1);
    let lo = pxs[n].get();
    let hi = pxs.get(n + 1).map(|p| p.get()).unwrap_or(lo);
    if hi == lo {
        return Ok(Unit::new(n as f64));
    }
    let frac = (t - lo) / (hi - lo);
    Ok(Unit::new(n as f64 + frac))
}

/// Rectilinear grid with an independent accumulator per axis.
/// Px passed in from outside is global (origin not yet subtracted) and each method subtracts `origin` itself;
/// Px derived from a boundary box (`unit_to_px` and what is built on it) is local.
pub struct RectGrid<const D: usize, P = DefaultSteps> {
    /// Origin of the grid in global px. Every method subtracts it from the px passed in.
    pub origin:  [Px; D],
    accumulator: [Accumulator<P>; D],
}

impl<const D: usize, P: StepFn> RectGrid<D, P> {
    /// Builds a grid from its origin and one definition per axis.
    /// Returns the first error of [`IncrementFunction::accumulate`].
    pub fn new(origin: [Px; D], definitions: [IncrementFunction<P>; D]) -> Result<Self, Error> {
        let accumulator: Vec<_> =
            definitions.into_iter().map(|d| d.accumulate()).collect::<Result<_, _>>()?;
        let accumulator = match accumulator.try_into() {
            Ok(a) => a,
            Err(_) => unreachable!("definitions and accumulator share length D"),
        };
        Ok(Self { origin, accumulator })
    }

    /// Replaces the definition for axis d.
    ///
    /// ```
    /// use rectgrid::{RectGrid, IncrementFunction, Px, Unit};
    /// let mut grid = RectGrid::<1>::new([Px::new(0.0)], [IncrementFunction::Scale(10.0)]).unwrap();
    /// assert_eq!(grid.unit_to_px(0, &Unit::new(2.0)).unwrap().get(), 20.0);
    /// grid.set_definition(IncrementFunction::Scale(100.0), 0).unwrap();
    /// assert_eq!(grid.unit_to_px(0, &Unit::new(2.0)).unwrap().get(), 200.0);
    /// ```
    pub fn set_definition(
        &mut self,
        definition: IncrementFunction<P>,
        d: usize,
    ) -> Result<(), Error> {
        self.accumulator[d] = definition.accumulate()?;
        Ok(())
    }

    /// Inverts px to unit.
    /// Caller contract: each axis's accumulator must be monotonically non-decreasing over Unit >= 0
    /// (e.g. not a negative Scale or a decreasing difference); otherwise the result is not guaranteed.
    ///
    /// ```
    /// use rectgrid::{RectGrid, IncrementFunction, Px};
    /// let grid = RectGrid::<1>::new([Px::new(0.0)], [IncrementFunction::Scale(10.0)]).unwrap();
    /// let result = grid.point_to_unit([Px::new(25.0)]);
    /// assert!((result[0].unwrap().get() - 2.5).abs() < 1e-6);
    /// ```
    pub fn point_to_unit(&self, point: [Px; D]) -> [Result<Unit, Error>; D] {
        from_fn(|i| self.px_to_unit_axis(i, point[i] - self.origin[i]))
    }

    fn px_to_unit_axis(&self, i: usize, target: Px) -> Result<Unit, Error> {
        self.accumulator[i].inverse(target)
    }

    /// Converts a unit coordinate to px. unit must be an absolute value from the origin.
    ///
    /// ```
    /// use rectgrid::{RectGrid, IncrementFunction, Px, Unit};
    /// let grid = RectGrid::<1>::new([Px::new(0.0)], [IncrementFunction::Scale(200.0)]).unwrap();
    /// assert_eq!(grid.unit_to_px(0, &Unit::new(2.25)).unwrap().get(), 450.0);
    /// ```
    pub fn unit_to_px(&self, d: usize, unit: &Unit) -> Result<Px, Error> {
        self.accumulator[d].forward(unit.get())
    }

    /// Converts multiple unit coordinate points to px. Each axis is evaluated independently:
    /// a point gives one `Result` per axis (Err, e.g. OutOfIndex(last), for an unevaluable axis).
    ///
    /// ```
    /// extern crate alloc;
    /// use rectgrid::{RectGrid, IncrementFunction, Px, Unit};
    /// use rectgrid::Error;
    /// let grid = RectGrid::<2>::new(
    ///     [Px::new(0.0), Px::new(0.0)],
    ///     [
    ///         IncrementFunction::Scale(100.0),
    ///         IncrementFunction::VectorList(alloc::vec![Px::new(0.0), Px::new(10.0)]),
    ///     ],
    /// ).unwrap();
    /// let points = alloc::vec![[Unit::new(2.0), Unit::new(0.5)], [Unit::new(2.0), Unit::new(5.0)]];
    /// let px = grid.point_as_px(&points);
    /// assert_eq!(px[0][1].as_ref().unwrap().get(), 5.0);
    /// // the 2nd point's second axis exceeds the VectorList's domain (0..=1): only that axis is Err
    /// assert_eq!(px[1][0].as_ref().unwrap().get(), 200.0);
    /// assert!(matches!(px[1][1], Err(Error::OutOfIndex(1))));
    /// ```
    pub fn point_as_px(&self, points: &[Point<D>]) -> Vec<[Result<Px, Error>; D]> {
        points.iter().map(|pt| from_fn(|d| self.unit_to_px(d, &pt[d]))).collect()
    }

    /// unit_to_px for an `extend` edge: past the end of a finite domain the edge is clipped at the last valid unit.
    fn unit_to_px_clipped(&self, d: usize, unit: &Unit) -> Result<Px, Error> {
        match self.unit_to_px(d, unit) {
            Err(Error::OutOfIndex(last)) => self.accumulator[d].forward(last as f64),
            other => other,
        }
    }

    /// Returns None when the boundary box is not evaluable (some axis of box_as_px is Err): it is never hit.
    /// Otherwise (whether it hit, base_px without extend, offset_px without extend), the px reusable for the parameter.
    /// extend is added in unit space before conversion to px, so the width does not drift under a nonlinear accumulator;
    /// an edge past the end of a finite domain is clipped there (see unit_to_px_clipped).
    fn contains(
        &self,
        point: [Px; D],
        bx: &BBox<D>,
        extend: Option<([Unit; D], [Unit; D])>,
    ) -> Option<(bool, [Px; D], [Px; D])> {
        let mut base_px = [Px::new(0.0); D];
        let mut offset_px = [Px::new(0.0); D];
        let mut hit = true;
        for d in 0..D {
            base_px[d] = self.unit_to_px(d, &bx.base[d]).ok()?;
            offset_px[d] = self.unit_to_px(d, &(bx.base[d] + bx.offset[d])).ok()?;
            let (lo, hi) = match extend {
                Some((eb, eo)) => (
                    self.unit_to_px_clipped(d, &(bx.base[d] + eb[d])).ok()?,
                    self.unit_to_px_clipped(d, &(bx.base[d] + bx.offset[d] + eo[d])).ok()?,
                ),
                None => (base_px[d], offset_px[d]),
            };
            let local = point[d] - self.origin[d];
            hit &= local.get() >= lo.get() && local.get() <= hi.get();
        }
        Some((hit, base_px, offset_px))
    }

    fn parameter_from_px(point: [Px; D], base_px: [Px; D], offset_px: [Px; D]) -> [Parameter; D] {
        from_fn(|d| Self::parameter_axis(point[d], base_px[d], offset_px[d]))
    }

    fn parameter_axis(point: Px, base_px: Px, far_px: Px) -> Parameter {
        let width = far_px - base_px;
        if width.get() == 0.0 {
            Parameter::new(0.0)
        } else {
            Parameter::new((point - base_px) / width)
        }
    }

    /// Returns the highest index among the boundary boxes that point hits (a higher index has higher priority).
    /// A boundary box that is not evaluable (some axis of box_as_px is Err) is never hit; an extend past a finite domain end is clipped there.
    ///
    /// ```
    /// extern crate alloc;
    /// use rectgrid::{RectGrid, IncrementFunction, BBox, Px, Unit};
    /// let grid = RectGrid::<2>::new(
    ///     [Px::new(0.0), Px::new(0.0)],
    ///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
    /// ).unwrap();
    /// let boxes = alloc::vec![
    ///     BBox::new([Unit::new(0.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(1.0)]),
    ///     BBox::new([Unit::new(2.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(1.0)]),
    /// ];
    /// assert_eq!(grid.hit_test([Px::new(450.0), Px::new(10.0)], &boxes, None), Some(1));
    /// assert_eq!(grid.hit_test([Px::new(-100.0), Px::new(-100.0)], &boxes, None), None);
    /// ```
    pub fn hit_test(
        &self,
        point: [Px; D],
        boxes: &[BBox<D>],
        extend: Option<([Unit; D], [Unit; D])>,
    ) -> Option<usize> {
        boxes.iter().enumerate().rev().find_map(|(i, bx)| {
            matches!(self.contains(point, bx, extend), Some((true, _, _))).then_some(i)
        })
    }

    /// Like hit_test, with the get_parameter value of the hit boundary box. The parameter is unaffected by extend.
    ///
    /// ```
    /// extern crate alloc;
    /// use rectgrid::{RectGrid, IncrementFunction, BBox, Px, Unit};
    /// let grid = RectGrid::<2>::new(
    ///     [Px::new(0.0), Px::new(0.0)],
    ///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
    /// ).unwrap();
    /// let boxes = alloc::vec![BBox::new([Unit::new(0.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(1.0)])];
    /// let (i, parameter) = grid.hit_test_with_parameter([Px::new(100.0), Px::new(32.0)], &boxes, None).unwrap();
    /// assert_eq!(i, 0);
    /// assert!((parameter[0].get() - 0.5).abs() < 1e-9);
    /// assert!((parameter[1].get() - 0.5).abs() < 1e-9);
    /// ```
    pub fn hit_test_with_parameter(
        &self,
        point: [Px; D],
        boxes: &[BBox<D>],
        extend: Option<([Unit; D], [Unit; D])>,
    ) -> Option<(usize, [Parameter; D])> {
        let local: [Px; D] = from_fn(|d| point[d] - self.origin[d]);
        boxes.iter().enumerate().rev().find_map(|(i, bx)| {
            let (hit, base_px, offset_px) = self.contains(point, bx, extend)?;
            hit.then(|| (i, Self::parameter_from_px(local, base_px, offset_px)))
        })
    }

    /// Hit/no-hit per boundary box, in order. A boundary box that is not evaluable is no-hit (see hit_test).
    ///
    /// ```
    /// extern crate alloc;
    /// use rectgrid::{RectGrid, IncrementFunction, BBox, Px, Unit};
    /// let grid = RectGrid::<2>::new(
    ///     [Px::new(0.0), Px::new(0.0)],
    ///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
    /// ).unwrap();
    /// let boxes = alloc::vec![
    ///     BBox::new([Unit::new(0.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(1.0)]),
    ///     BBox::new([Unit::new(0.5), Unit::new(0.0)], [Unit::new(1.0), Unit::new(1.0)]),
    /// ];
    /// assert_eq!(grid.hit_tests([Px::new(150.0), Px::new(10.0)], &boxes, None), alloc::vec![true, true]);
    /// ```
    pub fn hit_tests(
        &self,
        point: [Px; D],
        boxes: &[BBox<D>],
        extend: Option<([Unit; D], [Unit; D])>,
    ) -> Vec<bool> {
        boxes
            .iter()
            .map(|bx| matches!(self.contains(point, bx, extend), Some((true, _, _))))
            .collect()
    }

    /// `ξ_d = (point_d − base_d) / offset_d`, evaluated in px.
    /// Each axis is evaluated independently and gives its own `Result`: an axis whose base or base + offset
    /// cannot be evaluated is Err (OutOfIndex(last) carries the last valid unit, enough to clamp with unit_to_px).
    /// Non-finite inputs are not checked and give a NaN or infinite parameter inside Ok.
    ///
    /// ```
    /// use rectgrid::{RectGrid, IncrementFunction, BBox, Px, Unit};
    /// let grid = RectGrid::<2>::new(
    ///     [Px::new(0.0), Px::new(0.0)],
    ///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
    /// ).unwrap();
    /// let bx = BBox::new([Unit::new(1.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(1.0)]);
    /// let parameter = grid.get_parameter([Px::new(300.0), Px::new(16.0)], bx);
    /// assert!((parameter[0].as_ref().unwrap().get() - 0.5).abs() < 1e-9);
    /// assert!((parameter[1].as_ref().unwrap().get() - 0.25).abs() < 1e-9);
    /// // outside the boundary box (base side), parameter goes negative
    /// let parameter = grid.get_parameter([Px::new(100.0), Px::new(0.0)], bx);
    /// assert!((parameter[0].as_ref().unwrap().get() - (-0.5)).abs() < 1e-9);
    /// ```
    pub fn get_parameter(&self, point: [Px; D], bx: BBox<D>) -> [Result<Parameter, Error>; D] {
        from_fn(|d| -> Result<Parameter, Error> {
            let base = self.unit_to_px(d, &bx.base[d])?;
            let far = self.unit_to_px(d, &(bx.base[d] + bx.offset[d]))?;
            Ok(Self::parameter_axis(point[d] - self.origin[d], base, far))
        })
    }

    /// Converts multiple boundary boxes to (base_px, offset_px) per axis, with offset_px = unit_to_px(base+offset) - unit_to_px(base)
    /// (correct under a nonlinear accumulator). Each axis is evaluated independently: a boundary box gives one `Result` per axis
    /// (Err, e.g. OutOfIndex(last), for an unevaluable axis).
    ///
    /// ```
    /// extern crate alloc;
    /// use rectgrid::{RectGrid, IncrementFunction, BBox, Px, Unit};
    /// let grid = RectGrid::<1>::new([Px::new(0.0)], [IncrementFunction::Scale(100.0)]).unwrap();
    /// let boxes = alloc::vec![BBox::new([Unit::new(1.0)], [Unit::new(2.0)])];
    /// let result = grid.box_as_px(&boxes);
    /// let (base_px, offset_px) = result[0][0].as_ref().unwrap();
    /// assert_eq!(base_px.get(), 100.0);
    /// assert_eq!(offset_px.get(), 200.0);
    /// ```
    pub fn box_as_px(&self, boxes: &[BBox<D>]) -> Vec<[Result<(Px, Px), Error>; D]> {
        boxes
            .iter()
            .map(|bx| {
                from_fn(|d| -> Result<(Px, Px), Error> {
                    let base = self.unit_to_px(d, &bx.base[d])?;
                    let far = self.unit_to_px(d, &(bx.base[d] + bx.offset[d]))?;
                    Ok((base, far - base))
                })
            })
            .collect()
    }

    /// Returns pointer's local coordinate with z subtracted. At drag start, z = base_px yields the drag offset;
    /// during drag, z = that offset yields the element's current reference position.
    ///
    /// ```
    /// use rectgrid::{RectGrid, IncrementFunction, Px};
    /// let grid = RectGrid::<2>::new(
    ///     [Px::new(10.0), Px::new(20.0)],
    ///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
    /// ).unwrap();
    /// let drag_offset = grid.offset([Px::new(230.0), Px::new(50.0)], [Px::new(200.0), Px::new(0.0)]);
    /// assert_eq!((drag_offset[0].get(), drag_offset[1].get()), (20.0, 30.0));
    /// ```
    pub fn offset(&self, pointer: [Px; D], z: [Px; D]) -> [Px; D] {
        from_fn(|d| (pointer[d] - self.origin[d]) - z[d])
    }
}

/// For a boundary box with area, determines per axis whether point is near an edge (within threshold).
/// Returns (each axis's parameter, corner result). A corner element is Some(true) near the base-side edge
/// ([0, threshold]), Some(false) near the offset-side edge ([1 - threshold, 1]) and None otherwise.
/// If any axis is Some it is a handle hit (all axes Some = corner, one = edge); all None returns None
/// (the caller falls back to e.g. a move drag).
/// extend (as in hit_test) is converted to a parameter-space margin per axis (extend / offset).
/// parameter is None when the boundary box has no size, an axis is unevaluable, or point is outside bx (extend included).
///
/// ```
/// use rectgrid::{RectGrid, IncrementFunction, BBox, Px, Unit, corner_test};
/// let grid = RectGrid::<2>::new(
///     [Px::new(0.0), Px::new(0.0)],
///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
/// ).unwrap();
/// let bx = BBox::new([Unit::new(2.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(3.0)]);
/// // clicking near bx's top-left corner (400, 0)px -> corner on the base side for both axes
/// let (_, corner) = corner_test(&grid, [Px::new(400.0), Px::new(0.0)], &bx, 0.1, None);
/// assert_eq!(corner, Some([Some(true), Some(true)]));
/// // clicking near the middle of bx's top edge (x center, y = base side) -> a y-only handle (edge drag)
/// let (_, corner) = corner_test(&grid, [Px::new(500.0), Px::new(0.0)], &bx, 0.1, None);
/// assert_eq!(corner, Some([None, Some(true)]));
/// // near the center of bx matches neither edge nor corner, but parameter is still obtainable
/// let (parameter, corner) = corner_test(&grid, [Px::new(500.0), Px::new(96.0)], &bx, 0.1, None);
/// assert!(parameter.is_some());
/// assert_eq!(corner, None);
/// // outside bx, parameter is also unobtainable
/// let (parameter, corner) = corner_test(&grid, [Px::new(400.0), Px::new(-10.0)], &bx, 0.1, None);
/// assert!(parameter.is_none());
/// assert_eq!(corner, None);
/// // just outside the base-side edge, but within extend -> still treated as a base-side handle
/// let extend = Some(([Unit::new(-0.1), Unit::new(-0.1)], [Unit::new(0.1), Unit::new(0.1)]));
/// let (parameter, corner) = corner_test(&grid, [Px::new(400.0), Px::new(-5.0)], &bx, 0.1, extend);
/// assert!(parameter.is_some());
/// assert_eq!(corner, Some([Some(true), Some(true)]));
/// ```
pub fn corner_test<const D: usize, P: StepFn>(
    grid: &RectGrid<D, P>,
    point: [Px; D],
    bx: &BBox<D>,
    threshold: f64,
    extend: Option<([Unit; D], [Unit; D])>,
) -> (Option<[Parameter; D]>, Option<[Option<bool>; D]>) {
    if !bx.has_size() {
        return (None, None);
    }
    let mut parameter = [Parameter::new(0.0); D];
    for (axis, result) in grid.get_parameter(point, *bx).into_iter().enumerate() {
        match result {
            Ok(value) => parameter[axis] = value,
            Err(_) => return (None, None),
        }
    }
    let margin: [(f64, f64); D] = from_fn(|d| match extend {
        Some((eb, eo)) => {
            let offset = bx.offset()[d].get();
            if offset == 0.0 { (0.0, 0.0) } else { (eb[d].get() / offset, eo[d].get() / offset) }
        }
        None => (0.0, 0.0),
    });
    let inside =
        parameter.iter().zip(margin).all(|(r, (lo, hi))| r.get() >= lo && r.get() <= 1.0 + hi);
    if !inside {
        return (None, None);
    }
    let corner: [Option<bool>; D] = from_fn(|d| {
        let r = parameter[d].get();
        if r <= threshold {
            Some(true)
        } else if r >= 1.0 - threshold {
            Some(false)
        } else {
            None
        }
    });
    let corner = if corner.iter().any(Option::is_some) { Some(corner) } else { None };
    (Some(parameter), corner)
}

/// Updates the boundary box's base/offset via a corner-handle drag. corner\[d\] = Some(base_side): true moves the
/// base-side edge, false the offset-side edge. The new offset is at least 1.0 unit (base/offset never cross).
///
/// ```
/// use rectgrid::{RectGrid, IncrementFunction, BBox, Px, Unit, drag_resize};
/// let grid = RectGrid::<2>::new(
///     [Px::new(0.0), Px::new(0.0)],
///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
/// ).unwrap();
/// let bx = BBox::new([Unit::new(2.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(3.0)]);
/// // dragging the top-left (base-side) corner to (150, 0)px -> the x axis's base shrinks to 0 unit
/// let resized = drag_resize(&grid, [Px::new(150.0), Px::new(0.0)], &bx, [Some(true), None]).unwrap();
/// assert_eq!(resized.base()[0].get(), 0.0);
/// assert_eq!(resized.offset()[0].get(), 3.0); // base 2.0 + offset 1.0 - new_base 0.0
/// assert_eq!(resized.offset()[1].get(), 3.0); // y axis unchanged
/// ```
pub fn drag_resize<const D: usize, P: StepFn>(
    grid: &RectGrid<D, P>,
    pointer: [Px; D],
    bx: &BBox<D>,
    corner: [Option<bool>; D],
) -> Result<BBox<D>, Error> {
    let unit = grid.point_to_unit(pointer);
    let mut resized = *bx;
    for d in 0..D {
        let Some(base_side) = corner[d] else { continue };
        let new_u = Unit::new(libm::floor(unit[d]?.get()));
        let base_u = bx.base[d];
        let offset_u = bx.offset[d];
        if base_side {
            let far_edge = base_u + offset_u;
            let clamped_u = Unit::new(new_u.get().min((far_edge - Unit::new(1.0)).get()));
            resized.base[d] = clamped_u;
            resized.offset[d] = Unit::new((far_edge - clamped_u).get().max(1.0));
        } else {
            let new_offset = (new_u - base_u).get().max(1.0);
            resized.offset[d] = Unit::new(new_offset);
        }
    }
    Ok(resized)
}

/// Computes base's px position during a move drag. A boundary box keeps base in Unit, so the px position is only
/// returned; it is committed to the boundary box by snap_bbox_to_unit (or snap_point_to_unit) at drag end.
///
/// ```
/// use rectgrid::{RectGrid, IncrementFunction, Px, drag_translate};
/// let grid = RectGrid::<2>::new(
///     [Px::new(10.0), Px::new(20.0)],
///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
/// ).unwrap();
/// let px = drag_translate(&grid, [Px::new(230.0), Px::new(50.0)], [Px::new(20.0), Px::new(30.0)]);
/// assert_eq!((px[0].get(), px[1].get()), (200.0, 0.0));
/// ```
pub fn drag_translate<const D: usize, P: StepFn>(
    grid: &RectGrid<D, P>,
    pointer: [Px; D],
    drag_offset: [Px; D],
) -> [Px; D] {
    grid.offset(pointer, drag_offset)
}

/// At drag end, snaps the move-drag result of a boundary box with area to the Unit grid and returns the updated boundary box.
/// drag_offset must be the value passed to drag_translate. extend is added before flooring.
///
/// ```
/// use rectgrid::{RectGrid, IncrementFunction, BBox, Px, Unit, snap_bbox_to_unit};
/// let grid = RectGrid::<2>::new(
///     [Px::new(0.0), Px::new(0.0)],
///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
/// ).unwrap();
/// let bx = BBox::new([Unit::new(0.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(1.0)]);
/// let snapped = snap_bbox_to_unit(
///     &grid, [Px::new(430.0), Px::new(70.0)], [Px::new(0.0), Px::new(0.0)], &bx,
///     Some([Unit::new(0.25), Unit::new(0.25)]),
/// ).unwrap();
/// assert_eq!(snapped.base()[0].get(), 2.0);
/// assert_eq!(snapped.base()[1].get(), 1.0);
/// assert_eq!(snapped.offset()[0].get(), 1.0); // offset is already floored, so it stays as-is
/// ```
pub fn snap_bbox_to_unit<const D: usize, P: StepFn>(
    grid: &RectGrid<D, P>,
    pointer: [Px; D],
    drag_offset: [Px; D],
    bx: &BBox<D>,
    extend: Option<[Unit; D]>,
) -> Result<BBox<D>, Error> {
    let mut snapped = *bx;
    for d in 0..D {
        let local = pointer[d] - grid.origin[d] - drag_offset[d];
        snapped.base[d] = grid.px_to_unit_axis(d, local)?;
    }
    snapped.snap_floor(extend);
    Ok(snapped)
}

/// At drag end, computes a boundary box snapped to the Unit grid from the move-drag result of a boundary box without area (a point).
/// drag_offset must be the value passed to drag_translate.
///
/// ```
/// use rectgrid::{RectGrid, IncrementFunction, Px, Unit, snap_point_to_unit};
/// let grid = RectGrid::<2>::new(
///     [Px::new(0.0), Px::new(0.0)],
///     [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
/// ).unwrap();
/// let snapped = snap_point_to_unit(
///     &grid, [Px::new(430.0), Px::new(70.0)], [Px::new(0.0), Px::new(0.0)], [Unit::new(0.25), Unit::new(0.25)],
/// ).unwrap();
/// assert_eq!(snapped.base()[0].get(), 2.0);
/// assert_eq!(snapped.base()[1].get(), 1.0);
/// assert!(!snapped.has_size());
/// ```
pub fn snap_point_to_unit<const D: usize, P: StepFn>(
    grid: &RectGrid<D, P>,
    pointer: [Px; D],
    drag_offset: [Px; D],
    snap: [Unit; D],
) -> Result<BBox<D>, Error> {
    let mut base: Point<D> = [Unit::new(0.0); D];
    for d in 0..D {
        let local = pointer[d] - grid.origin[d] - drag_offset[d];
        let u = grid.px_to_unit_axis(d, local)?;
        base[d] = Unit::new(libm::floor((u + snap[d]).get()));
    }
    Ok(BBox { base, offset: from_fn(|_| Unit::new(0.0)) })
}

#[cfg(test)]
mod tests {
    use alloc::rc::Rc;

    use super::*;

    #[test]
    fn point_to_unit_scale_roundtrip() {
        let grid = RectGrid::<2>::new(
            [Px::new(0.0), Px::new(0.0)],
            [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
        )
        .unwrap();
        let result = grid.point_to_unit([Px::new(450.0), Px::new(10.0)]);
        let x = result[0].as_ref().unwrap().get();
        let y = result[1].as_ref().unwrap().get();
        assert!((x - 2.25).abs() < 1e-6, "x = {}", x);
        assert!((y - 0.15625).abs() < 1e-6, "y = {}", y);
    }

    #[test]
    fn point_to_unit_vector_list_roundtrip() {
        let grid = RectGrid::<1>::new(
            [Px::new(0.0)],
            [IncrementFunction::VectorList(alloc::vec![
                Px::new(0.0),
                Px::new(10.0),
                Px::new(30.0),
                Px::new(60.0)
            ])],
        )
        .unwrap();
        let result = grid.point_to_unit([Px::new(45.0)]);
        let x = result[0].as_ref().unwrap().get();
        assert!((x - 2.5).abs() < 1e-6, "x = {}", x);
    }

    #[test]
    fn point_to_unit_vector_list_out_of_range() {
        let grid = RectGrid::<1>::new(
            [Px::new(0.0)],
            [IncrementFunction::VectorList(alloc::vec![
                Px::new(0.0),
                Px::new(10.0),
                Px::new(30.0),
                Px::new(60.0)
            ])],
        )
        .unwrap();
        let result = grid.point_to_unit([Px::new(100.0)]);
        assert!(matches!(result[0], Err(Error::OutOfIndex(3))));
    }

    #[test]
    fn point_to_unit_applies_origin_offset() {
        let grid = RectGrid::<2>::new(
            [Px::new(10.0), Px::new(20.0)],
            [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
        )
        .unwrap();
        let result = grid.point_to_unit([Px::new(230.0), Px::new(50.0)]);
        let x = result[0].as_ref().unwrap().get();
        let y = result[1].as_ref().unwrap().get();
        assert!((x - 1.1).abs() < 1e-6, "x = {}", x);
        assert!((y - 0.46875).abs() < 1e-6, "y = {}", y);
    }

    /// f(i) = (i+1)*10, so forward(2.5) = 10 + 20 + 0.5 * 30.
    #[test]
    fn accumulate_forward_difference() {
        let f =
            IncrementFunction::ForwardDifference(Rc::new(|i| Ok(Px::new((i + 1) as f64 * 10.0))));
        let acc = f.accumulate().unwrap();
        assert_eq!(acc.forward(0.0).unwrap().get(), 0.0);
        assert!((acc.forward(2.5).unwrap().get() - 45.0).abs() < 1e-9);
    }

    #[test]
    fn accumulator_scale_inverse_is_exact() {
        let acc = IncrementFunction::<DefaultSteps>::Scale(200.0).accumulate().unwrap();
        assert_eq!(acc.inverse(Px::new(450.0)).unwrap().get(), 2.25);
    }

    #[test]
    fn accumulator_vector_list_inverse_boundaries() {
        let acc = IncrementFunction::<DefaultSteps>::VectorList(alloc::vec![
            Px::new(0.0),
            Px::new(10.0),
            Px::new(30.0)
        ])
        .accumulate()
        .unwrap();
        assert_eq!(acc.inverse(Px::new(0.0)).unwrap().get(), 0.0);
        assert_eq!(acc.inverse(Px::new(30.0)).unwrap().get(), 2.0);
        assert!((acc.inverse(Px::new(20.0)).unwrap().get() - 1.5).abs() < 1e-9);
        assert!(matches!(acc.inverse(Px::new(-1.0)), Err(Error::OutOfIndex(0))));
        assert!(matches!(acc.inverse(Px::new(31.0)), Err(Error::OutOfIndex(2))));
    }

    #[test]
    fn bbox_new_normalizes_negative_offset() {
        let bx = BBox::new([Unit::new(0.0)], [Unit::new(-1.0)]);
        assert_eq!(bx.base()[0].get(), -1.0);
        assert_eq!(bx.offset()[0].get(), 1.0);

        let bx = BBox::new([Unit::new(0.0)], [Unit::new(0.0)]);
        assert_eq!(bx.base()[0].get(), 0.0);
        assert_eq!(bx.offset()[0].get(), 0.0);
    }

    #[test]
    fn snap_floor_no_extend() {
        let mut bx = BBox { base: [Unit::new(2.7)], offset: [Unit::new(1.9)] };
        bx.snap_floor(None);
        assert_eq!(bx.base[0].get(), 2.0);
        assert_eq!(bx.offset[0].get(), 1.0);
    }

    #[test]
    fn box_as_px_out_of_index() {
        let grid = RectGrid::<1>::new(
            [Px::new(0.0)],
            [IncrementFunction::VectorList(alloc::vec![Px::new(0.0), Px::new(10.0)])],
        )
        .unwrap();
        let boxes = alloc::vec![BBox { base: [Unit::new(0.0)], offset: [Unit::new(5.0)] }];
        assert!(matches!(grid.box_as_px(&boxes)[0][0], Err(Error::OutOfIndex(1))));
    }

    #[test]
    fn hit_test_with_parameter_no_hit() {
        let grid = RectGrid::<2>::new(
            [Px::new(0.0), Px::new(0.0)],
            [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
        )
        .unwrap();
        let boxes = alloc::vec![BBox {
            base:   [Unit::new(0.0), Unit::new(0.0)],
            offset: [Unit::new(1.0), Unit::new(1.0)],
        },];
        assert!(
            grid.hit_test_with_parameter([Px::new(500.0), Px::new(10.0)], &boxes, None).is_none()
        );
    }

    #[test]
    fn hit_tests_partial_and_no_hit() {
        let grid = RectGrid::<2>::new(
            [Px::new(0.0), Px::new(0.0)],
            [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
        )
        .unwrap();
        let boxes = alloc::vec![
            BBox {
                base:   [Unit::new(0.0), Unit::new(0.0)],
                offset: [Unit::new(1.0), Unit::new(1.0)],
            },
            BBox {
                base:   [Unit::new(2.0), Unit::new(0.0)],
                offset: [Unit::new(1.0), Unit::new(1.0)],
            },
        ];
        assert_eq!(
            grid.hit_tests([Px::new(50.0), Px::new(10.0)], &boxes, None),
            alloc::vec![true, false]
        );
        assert_eq!(
            grid.hit_tests([Px::new(-1.0), Px::new(-1.0)], &boxes, None),
            alloc::vec![false, false]
        );
    }

    /// Dragging the offset-side edge to 900px (unit 4.5, floored to 4) gives offset 4 - 2 with base unchanged.
    #[test]
    fn drag_resize_offset_side() {
        let grid = RectGrid::<2>::new(
            [Px::new(0.0), Px::new(0.0)],
            [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
        )
        .unwrap();
        let bx = BBox {
            base:   [Unit::new(2.0), Unit::new(0.0)],
            offset: [Unit::new(1.0), Unit::new(3.0)],
        };
        let resized =
            drag_resize(&grid, [Px::new(900.0), Px::new(0.0)], &bx, [Some(false), None]).unwrap();
        assert_eq!(resized.base[0].get(), 2.0);
        assert_eq!(resized.offset[0].get(), 2.0);
    }

    /// Dragging the base side past the far edge (unit 5.5) clamps base below it and offset to the minimum of 1.0.
    #[test]
    fn drag_resize_minimum_clamp() {
        let grid = RectGrid::<2>::new(
            [Px::new(0.0), Px::new(0.0)],
            [IncrementFunction::Scale(200.0), IncrementFunction::Scale(64.0)],
        )
        .unwrap();
        let bx = BBox {
            base:   [Unit::new(2.0), Unit::new(0.0)],
            offset: [Unit::new(2.0), Unit::new(3.0)],
        };
        let resized =
            drag_resize(&grid, [Px::new(1100.0), Px::new(0.0)], &bx, [Some(true), None]).unwrap();
        assert_eq!(resized.base[0].get(), 3.0);
        assert_eq!(resized.offset[0].get(), 1.0);
    }
    type Steps = Rc<dyn Fn(u32) -> Result<Px, Error>>;

    /// A non-uniform increment (1.0 ... 1.016) that has no closed form.
    fn wavy(k: u32) -> Result<Px, Error> {
        Ok(Px::new(1.0 + libm::fabs(libm::sin(k as f64 * 0.37)) * 0.01 + 1e-3 * (k % 7) as f64))
    }

    /// 1.0001^k: monotonically growing steps.
    fn growing(k: u32) -> Result<Px, Error> {
        Ok(Px::new(libm::pow(1.0001, k as f64)))
    }

    /// Hashed steps in [0.5, 1.5): irregular, deterministic.
    fn irregular(k: u32) -> Result<Px, Error> {
        let mut z = (k as u64).wrapping_add(0x9E3779B97F4A7C15);
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58476D1CE4E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D049BB133111EB);
        z ^= z >> 31;
        Ok(Px::new(0.5 + (z >> 11) as f64 / (1u64 << 53) as f64))
    }

    fn accumulator(steps: Steps) -> Accumulator {
        IncrementFunction::ForwardDifference(steps).accumulate().unwrap()
    }

    #[test]
    fn forward_difference_inverse_roundtrip_is_exact() {
        let datasets: [(&str, Steps, &[f64]); 6] = [
            (
                "wavy",
                Rc::new(wavy),
                &[0.25, 0.999, 1.0, 3.5, 99.75, 1000.125, 12345.0137, 50_000.5],
            ),
            ("constant 1e-9", Rc::new(|_| Ok(Px::new(1e-9))), &[0.25, 3.5, 1000.125, 20_000.75]),
            ("constant 1e9", Rc::new(|_| Ok(Px::new(1e9))), &[0.25, 3.5, 1000.125, 20_000.75]),
            ("growing 1.0001^k", Rc::new(growing), &[0.5, 17.25, 1000.5, 4000.75]),
            ("irregular", Rc::new(irregular), &[0.1, 1.9, 50.5, 1234.5678, 30_000.25]),
            (
                "first step zero",
                Rc::new(|k| Ok(Px::new(if k == 0 { 0.0 } else { 2.0 }))),
                &[1.5, 10.25, 999.5],
            ),
        ];
        for (name, steps, xs) in datasets {
            let acc = accumulator(steps);
            for &x in xs {
                let back = acc.inverse(acc.forward(x).unwrap()).unwrap().get();
                assert!((back - x).abs() <= 4.0 * f64::EPSILON * x, "{name}: x = {x}: {back}");
            }
        }
    }

    #[test]
    fn forward_difference_inverse_sweeps_many_x_exactly() {
        for steps in [Rc::new(wavy) as Steps, Rc::new(irregular)] {
            let acc = accumulator(steps);
            let mut worst = 0.0f64;
            for k in 0..2000 {
                let x = 12345.0 + k as f64 * 0.0137;
                let back = acc.inverse(acc.forward(x).unwrap()).unwrap().get();
                worst = worst.max((back - x).abs());
            }
            // ulp(12345) = 1.8e-12
            assert!(worst <= 1e-11, "worst roundtrip error {worst:e}");
        }
    }

    #[test]
    fn forward_difference_inverse_constant_steps_are_exact() {
        let acc = accumulator(Rc::new(|_| Ok(Px::new(10.0))));
        assert_eq!(acc.inverse(Px::new(25.0)).unwrap().get(), 2.5);
        assert_eq!(acc.inverse(Px::new(10.0)).unwrap().get(), 1.0);
        assert_eq!(acc.inverse(Px::new(0.0)).unwrap().get(), 0.0);
        assert_eq!(acc.inverse(Px::new(-5.0)).unwrap().get(), -0.5);
    }

    /// Steps 10, 0, 10: px 10 is reached at unit 1 and held until unit 2; the inverse returns the start of the plateau.
    #[test]
    fn forward_difference_inverse_zero_step_returns_start_of_plateau() {
        let acc = accumulator(Rc::new(|k| Ok(Px::new(if k == 1 { 0.0 } else { 10.0 }))));
        assert_eq!(acc.inverse(Px::new(10.0)).unwrap().get(), 1.0);
        assert_eq!(acc.inverse(Px::new(15.0)).unwrap().get(), 2.5);
    }

    /// Three steps of 10 (domain 0..=3); the closure reports OutOfIndex(3) past the end.
    #[test]
    fn forward_difference_inverse_finite_domain_boundaries() {
        let acc =
            accumulator(Rc::new(
                |k| {
                    if k < 3 { Ok(Px::new(10.0)) } else { Err(Error::OutOfIndex(3)) }
                },
            ));
        assert_eq!(acc.inverse(Px::new(25.0)).unwrap().get(), 2.5);
        assert_eq!(acc.inverse(Px::new(30.0)).unwrap().get(), 3.0);
        assert!(matches!(acc.inverse(Px::new(30.5)), Err(Error::OutOfIndex(3))));
    }

    #[test]
    fn nan_is_rejected_as_invalid_input_for_every_definition() {
        let steps: DefaultSteps = Rc::new(|_| Ok(Px::new(1.0)));
        let defs = [
            IncrementFunction::Scale(2.0),
            IncrementFunction::VectorList(alloc::vec![Px::new(0.0), Px::new(1.0)]),
            IncrementFunction::ForwardDifference(steps),
        ];
        for def in defs {
            let acc = def.accumulate().unwrap();
            assert!(matches!(acc.forward(f64::NAN), Err(Error::InvalidInput)));
            assert!(matches!(acc.inverse(Px::new(f64::NAN)), Err(Error::InvalidInput)));
        }
        let nan_scale: Result<Accumulator, _> = IncrementFunction::Scale(f64::NAN).accumulate();
        assert!(matches!(nan_scale, Err(Error::InvalidInput)));
    }

    #[test]
    fn forward_difference_negative_side_round_trips_through_the_first_step() {
        let acc = accumulator(Rc::new(|i| Ok(Px::new((i + 1) as f64 * 10.0))));
        for x in [-0.5, -1.0, -3.25] {
            let px = acc.forward(x).unwrap();
            assert_eq!(px.get(), 10.0 * x);
            assert!((acc.inverse(px).unwrap().get() - x).abs() < 1e-12);
        }
        let flat = accumulator(Rc::new(|_| Ok(Px::new(0.0))));
        assert!(matches!(flat.inverse(Px::new(-1.0)), Err(Error::OutOfIndex(0))));
        assert_eq!(flat.inverse(Px::new(0.0)).unwrap().get(), 0.0);
    }

    #[test]
    fn value_supports_comparison_debug_and_sign_operations() {
        use alloc::format;
        let a = Px::new(1.5);
        let b = Px::new(2.5);
        assert!(a < b && b > a && a <= a && a != b);
        assert_eq!(a, Px::new(1.5));
        assert_eq!(Px::default(), Px::new(0.0));
        assert_eq!(-a, Px::new(-1.5));
        assert_eq!(2.0 * a, Px::new(3.0));
        assert_eq!(b / 2.0, Px::new(1.25));
        let mut c = b;
        c -= a;
        assert_eq!(c, Px::new(1.0));
        assert_eq!(format!("{:?}", a), "Value(1.5)");
        assert_eq!(Px::new(f64::NAN).partial_cmp(&a), None);
        assert!(Px::new(f64::NAN) != Px::new(f64::NAN));
    }

    #[test]
    fn scale_accepts_negative_and_infinite_values() {
        let acc = IncrementFunction::<DefaultSteps>::Scale(-2.0).accumulate().unwrap();
        assert_eq!(acc.forward(3.0).unwrap().get(), -6.0);
        assert_eq!(acc.inverse(Px::new(-6.0)).unwrap().get(), 3.0);
        assert_eq!(acc.forward(f64::INFINITY).unwrap().get(), f64::NEG_INFINITY);
    }

    /// A non-finite target is rejected (an unbounded closure would be scanned for u32::MAX segments).
    #[test]
    fn forward_difference_inverse_rejects_non_finite_targets() {
        let acc = accumulator(Rc::new(|_| Ok(Px::new(1.0))));
        for target in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
            assert!(matches!(acc.inverse(Px::new(target)), Err(Error::InvalidInput)));
        }
    }

    #[test]
    fn forward_difference_inverse_cost_is_one_call_per_segment_up_to_the_answer() {
        use core::cell::Cell;
        let calls = Rc::new(Cell::new(0u32));
        let counter = calls.clone();
        let acc = accumulator(Rc::new(move |_| {
            counter.set(counter.get() + 1);
            Ok(Px::new(2.0))
        }));
        acc.inverse(Px::new(2000.5)).unwrap(); // segment 1000
        assert_eq!(calls.get(), 1001);
    }

    /// The caller chooses the pointer: an `Arc<dyn Fn + Send + Sync>` makes the grid Send + Sync.
    #[test]
    fn arc_pointer_makes_the_grid_send_and_sync() {
        use alloc::sync::Arc;

        type SyncSteps = Arc<dyn Fn(u32) -> Result<Px, Error> + Send + Sync>;

        fn assert_send_sync<T: Send + Sync>() {}
        assert_send_sync::<RectGrid<2, SyncSteps>>();

        let steps: SyncSteps = Arc::new(|k| Ok(Px::new(10.0 + k as f64)));
        let grid: RectGrid<2, SyncSteps> = RectGrid::new(
            [Px::new(0.0), Px::new(0.0)],
            [IncrementFunction::ForwardDifference(steps), IncrementFunction::Scale(2.0)],
        )
        .unwrap();
        assert_eq!(grid.unit_to_px(0, &Unit::new(2.0)).unwrap().get(), 21.0);
        assert_eq!(grid.unit_to_px(1, &Unit::new(3.0)).unwrap().get(), 6.0);
    }

    fn vector_list_grid(pxs: &[f64]) -> RectGrid<1> {
        RectGrid::<1>::new(
            [Px::new(0.0)],
            [IncrementFunction::VectorList(pxs.iter().map(|&v| Px::new(v)).collect())],
        )
        .unwrap()
    }

    fn unit_box(base: f64, offset: f64) -> BBox<1> {
        BBox::new([Unit::new(base)], [Unit::new(offset)])
    }

    /// VectorList [0, 10] has the domain units 0..=1; a boundary box at units 5..6 is never hit.
    #[test]
    fn hit_test_box_beyond_finite_domain_never_hits() {
        let grid = vector_list_grid(&[0.0, 10.0]);
        let boxes = alloc::vec![unit_box(5.0, 1.0)];
        for px in [-1.0, 0.0, 5.0, 10.0, 50.0] {
            assert_eq!(grid.hit_test([Px::new(px)], &boxes, None), None, "px {px}");
            assert_eq!(grid.hit_tests([Px::new(px)], &boxes, None), alloc::vec![false], "px {px}");
            assert!(grid.hit_test_with_parameter([Px::new(px)], &boxes, None).is_none(), "px {px}");
        }
        let boxes = alloc::vec![unit_box(0.0, 1.0), unit_box(5.0, 1.0)];
        assert_eq!(grid.hit_test([Px::new(0.0)], &boxes, None), Some(0));
        assert_eq!(grid.hit_tests([Px::new(5.0)], &boxes, None), alloc::vec![true, false]);
    }

    /// Base inside the domain and base + offset beyond it: box_as_px reports Err on that axis.
    #[test]
    fn hit_test_box_straddling_the_domain_end_is_not_evaluable() {
        let grid = vector_list_grid(&[0.0, 10.0]);
        let boxes = alloc::vec![unit_box(0.5, 1.0)];
        assert!(matches!(grid.box_as_px(&boxes)[0][0], Err(Error::OutOfIndex(1))));
        assert_eq!(grid.hit_test([Px::new(7.0)], &boxes, None), None);
    }

    #[test]
    fn get_parameter_of_unevaluable_box_is_err_on_that_axis() {
        let grid = vector_list_grid(&[0.0, 10.0]);
        for bx in [unit_box(5.0, 1.0), unit_box(0.5, 1.0)] {
            let parameter = grid.get_parameter([Px::new(5.0)], bx);
            assert!(matches!(parameter[0], Err(Error::OutOfIndex(1))));
        }
        let parameter = grid.get_parameter([Px::new(5.0)], unit_box(0.0, 1.0));
        assert_eq!(parameter[0].as_ref().unwrap().get(), 0.5);
    }

    #[test]
    fn corner_test_reports_nothing_for_an_unevaluable_box() {
        let grid = vector_list_grid(&[0.0, 10.0]);
        let bx = unit_box(5.0, 1.0);
        for px in [0.0, 5.0, 10.0] {
            assert!(corner_test(&grid, [Px::new(px)], &bx, 0.1, None).0.is_none());
            assert!(corner_test(&grid, [Px::new(px)], &bx, 0.1, None).1.is_none());
        }
    }

    /// The boundary box is the last cell (px 10..30) of VectorList [0, 10, 30] and extend reaches past the domain end: it is clipped there. The parameter is relative to the boundary box, not to extend.
    #[test]
    fn hit_test_extend_is_clipped_at_the_domain_end() {
        let grid = vector_list_grid(&[0.0, 10.0, 30.0]);
        let boxes = alloc::vec![unit_box(1.0, 1.0)];
        let extend = Some(([Unit::new(0.0)], [Unit::new(5.0)]));
        assert_eq!(grid.hit_test([Px::new(20.0)], &boxes, extend), Some(0));
        assert_eq!(grid.hit_test([Px::new(30.0)], &boxes, extend), Some(0));
        assert_eq!(grid.hit_test([Px::new(30.5)], &boxes, extend), None);
        assert_eq!(grid.hit_test([Px::new(9.0)], &boxes, extend), None);
        let extend = Some(([Unit::new(-0.5)], [Unit::new(0.0)]));
        assert_eq!(grid.hit_test([Px::new(7.0)], &boxes, extend), Some(0));
        assert_eq!(grid.hit_test([Px::new(4.0)], &boxes, extend), None);
        let extend = Some(([Unit::new(0.0)], [Unit::new(5.0)]));
        let (_, parameter) = grid.hit_test_with_parameter([Px::new(20.0)], &boxes, extend).unwrap();
        assert_eq!(parameter[0].get(), 0.5);
    }

    #[test]
    fn vector_list_below_the_start_is_out_of_index_zero_and_extend_clips_there() {
        let grid = vector_list_grid(&[0.0, 10.0, 30.0]);
        assert!(matches!(grid.unit_to_px(0, &Unit::new(-0.5)), Err(Error::OutOfIndex(0))));
        assert!(matches!(grid.point_to_unit([Px::new(-1.0)])[0], Err(Error::OutOfIndex(0))));
        let boxes = alloc::vec![unit_box(0.0, 1.0)];
        let extend = Some(([Unit::new(-5.0)], [Unit::new(0.0)]));
        assert_eq!(grid.hit_test([Px::new(0.0)], &boxes, extend), Some(0));
        assert_eq!(grid.hit_test([Px::new(-1.0)], &boxes, extend), None);
    }

    /// Three steps of 10 (domain 0..=3); extend is clipped at the domain end and a boundary box starting beyond it is not evaluable.
    #[test]
    fn hit_test_extend_is_clipped_for_forward_difference_domains() {
        let grid = RectGrid::<1>::new(
            [Px::new(0.0)],
            [IncrementFunction::ForwardDifference(Rc::new(|k| {
                if k < 3 { Ok(Px::new(10.0)) } else { Err(Error::OutOfIndex(3)) }
            }))],
        )
        .unwrap();
        let boxes = alloc::vec![unit_box(2.0, 1.0)];
        let extend = Some(([Unit::new(0.0)], [Unit::new(2.0)]));
        assert_eq!(grid.hit_test([Px::new(25.0)], &boxes, extend), Some(0));
        assert_eq!(grid.hit_test([Px::new(30.0)], &boxes, extend), Some(0));
        assert_eq!(grid.hit_test([Px::new(31.0)], &boxes, extend), None);
        assert_eq!(grid.hit_test([Px::new(0.0)], &alloc::vec![unit_box(4.0, 1.0)], None), None);
    }

    #[test]
    fn hit_test_one_unevaluable_axis_is_enough_to_miss() {
        let grid = RectGrid::<2>::new(
            [Px::new(0.0), Px::new(0.0)],
            [
                IncrementFunction::Scale(100.0),
                IncrementFunction::VectorList(alloc::vec![Px::new(0.0), Px::new(10.0)]),
            ],
        )
        .unwrap();
        let beyond = alloc::vec![BBox::new(
            [Unit::new(0.0), Unit::new(5.0)],
            [Unit::new(1.0), Unit::new(1.0)]
        )];
        assert_eq!(grid.hit_test([Px::new(50.0), Px::new(0.0)], &beyond, None), None);
        let inside = alloc::vec![BBox::new(
            [Unit::new(0.0), Unit::new(0.0)],
            [Unit::new(1.0), Unit::new(1.0)]
        )];
        assert_eq!(grid.hit_test([Px::new(50.0), Px::new(5.0)], &inside, None), Some(0));
    }

    /// Only OutOfIndex is clipped at a domain end: other evaluation errors make the boundary box unevaluable.
    #[test]
    fn hit_test_other_evaluation_errors_also_mean_not_evaluable() {
        let grid = RectGrid::<1>::new(
            [Px::new(0.0)],
            [IncrementFunction::ForwardDifference(Rc::new(|_| Err(Error::InvalidDefinition)))],
        )
        .unwrap();
        let boxes = alloc::vec![unit_box(1.0, 1.0)];
        assert_eq!(grid.hit_test([Px::new(0.0)], &boxes, None), None);
        assert!(matches!(
            grid.get_parameter([Px::new(0.0)], boxes[0])[0],
            Err(Error::InvalidDefinition)
        ));
        let extend = Some(([Unit::new(0.0)], [Unit::new(1.0)]));
        assert_eq!(grid.hit_test([Px::new(0.0)], &boxes, extend), None);
    }

    fn assert_grid_lines_exact(name: &str, acc: &Accumulator, count: u32) {
        for k in 0..count {
            let px = acc.forward(k as f64).unwrap();
            let back = acc.inverse(px).unwrap().get();
            assert_eq!(back, k as f64, "{name}: k = {k}: {back:?}");
        }
    }

    #[test]
    fn scale_inverse_returns_grid_lines_exactly() {
        for scale in [0.1, 1.0 / 3.0, 1.1, 0.7, 1.5, 64.0, 200.0, 3.0e-7, 4.0e6] {
            assert_grid_lines_exact(
                &alloc::format!("Scale({scale})"),
                &IncrementFunction::Scale(scale).accumulate().unwrap(),
                100_000,
            );
        }
    }

    #[test]
    fn forward_difference_inverse_returns_grid_lines_exactly() {
        for (name, steps) in [
            ("wavy", Rc::new(wavy) as Steps),
            ("growing", Rc::new(growing)),
            ("irregular", Rc::new(irregular)),
            ("constant 0.1", Rc::new(|_| Ok(Px::new(0.1)))),
            ("constant 1/3", Rc::new(|_| Ok(Px::new(1.0 / 3.0)))),
        ] {
            assert_grid_lines_exact(name, &accumulator(steps), 3000);
        }
    }

    #[test]
    fn vector_list_inverse_returns_grid_lines_exactly() {
        let tenths: alloc::vec::Vec<Px> = (0..2000).map(|k| Px::new(k as f64 * 0.1)).collect();
        assert_grid_lines_exact(
            "k * 0.1",
            &IncrementFunction::VectorList(tenths).accumulate().unwrap(),
            2000,
        );
        let mut position = 0.0;
        let irregular_positions: alloc::vec::Vec<Px> = (0..2000)
            .map(|k| {
                if k > 0 {
                    position += 0.1 + (k % 7) as f64 * 0.013;
                }
                Px::new(position)
            })
            .collect();
        assert_grid_lines_exact(
            "irregular positions",
            &IncrementFunction::VectorList(irregular_positions).accumulate().unwrap(),
            2000,
        );
    }

    /// 0.3 is not forward(3) = 0.30000000000000004, so it is not snapped to 3.
    #[test]
    fn scale_inverse_off_the_grid_line_is_still_the_plain_quotient() {
        let acc = IncrementFunction::<DefaultSteps>::Scale(0.1).accumulate().unwrap();
        assert_eq!(acc.inverse(Px::new(0.3)).unwrap().get(), 0.3 / 0.1);
        assert_eq!(acc.inverse(Px::new(0.25)).unwrap().get(), 2.5);
    }

    #[test]
    fn drag_without_movement_keeps_the_box_on_a_fractional_scale() {
        let grid = RectGrid::<1>::new([Px::new(0.0)], [IncrementFunction::Scale(0.1)]).unwrap();
        for k in 0..1000 {
            let bx = BBox::new([Unit::new(k as f64)], [Unit::new(2.0)]);
            let base_px = grid.unit_to_px(0, &Unit::new(k as f64)).unwrap();
            let snapped = snap_bbox_to_unit(&grid, [base_px], [Px::new(0.0)], &bx, None).unwrap();
            assert_eq!(snapped.base()[0].get(), k as f64, "k = {k}");
            let point =
                snap_point_to_unit(&grid, [base_px], [Px::new(0.0)], [Unit::new(0.0)]).unwrap();
            assert_eq!(point.base()[0].get(), k as f64, "k = {k}");
        }
    }

    #[test]
    fn drag_resize_to_a_grid_line_lands_on_that_line_on_a_fractional_scale() {
        let grid = RectGrid::<1>::new([Px::new(0.0)], [IncrementFunction::Scale(0.1)]).unwrap();
        let bx = BBox::new([Unit::new(2.0)], [Unit::new(3.0)]);
        for k in 3..1000 {
            let pointer = grid.unit_to_px(0, &Unit::new(k as f64)).unwrap();
            let resized = drag_resize(&grid, [pointer], &bx, [Some(false)]).unwrap();
            assert_eq!(resized.offset()[0].get(), (k - 2) as f64, "k = {k}");
        }
    }

    #[test]
    fn a_failed_axis_can_be_clamped_with_the_last_index_it_reports() {
        let grid = vector_list_grid(&[0.0, 10.0, 30.0]);
        let Err(Error::OutOfIndex(last)) = grid.unit_to_px(0, &Unit::new(7.0)) else {
            panic!("expected OutOfIndex");
        };
        assert_eq!(last, 2);
        assert_eq!(grid.unit_to_px(0, &Unit::new(last as f64)).unwrap().get(), 30.0);
    }

    /// Axis 0 would be inside the boundary box but axis 1 lies beyond its finite domain: nothing is reported.
    #[test]
    fn corner_test_reports_nothing_when_any_axis_cannot_be_evaluated() {
        let grid = RectGrid::<2>::new(
            [Px::new(0.0), Px::new(0.0)],
            [
                IncrementFunction::Scale(100.0),
                IncrementFunction::VectorList(alloc::vec![Px::new(0.0), Px::new(10.0)]),
            ],
        )
        .unwrap();
        let beyond = BBox::new([Unit::new(0.0), Unit::new(5.0)], [Unit::new(1.0), Unit::new(1.0)]);
        let (parameter, corner) =
            corner_test(&grid, [Px::new(50.0), Px::new(5.0)], &beyond, 0.1, None);
        assert!(parameter.is_none() && corner.is_none());
        let inside = BBox::new([Unit::new(0.0), Unit::new(0.0)], [Unit::new(1.0), Unit::new(1.0)]);
        let (parameter, _) = corner_test(&grid, [Px::new(50.0), Px::new(5.0)], &inside, 0.1, None);
        assert!(parameter.is_some());
    }

    /// The boundary box and point lists are taken as slices: arrays and sub-slices need no Vec.
    #[test]
    fn list_arguments_accept_arrays_and_sub_slices() {
        let grid = vector_list_grid(&[0.0, 10.0, 30.0]);
        let all = alloc::vec![unit_box(0.0, 1.0), unit_box(1.0, 1.0)];
        assert_eq!(grid.hit_tests([Px::new(5.0)], &[unit_box(0.0, 1.0)], None), alloc::vec![true]);
        assert_eq!(grid.hit_test([Px::new(20.0)], &all[1..], None), Some(0));
        assert_eq!(grid.box_as_px(&all[..1]).len(), 1);
        assert_eq!(grid.point_as_px(&[[Unit::new(1.0)]])[0][0].as_ref().unwrap().get(), 10.0);
    }
}