delaunay 0.8.3

//! Spherical coordinate and metric backend.
//!
//! This module provides fallible coordinate and metric types for points on
//! `S^D` realized in `R^(D+1)`.

#![forbid(unsafe_code)]

use la_stack::{LaError, VectorAngle};
use thiserror::Error;

/// Typed point on `S^D` realized in `R^(D+1)`.
///
/// Rust's stable const generics do not allow a public `[f64; D + 1]` field, so
/// this topology backend stores the ambient coordinates in a length-checked
/// vector. Construction proves that the vector length is exactly `D + 1`,
/// every coordinate is finite, and the point has been normalized to the
/// requested radius.
/// Radius scaling must preserve each computed unit coordinate to within
/// `4 * f64::EPSILON` absolute error when divided back by the radius. This
/// rejects subnormal radii that cannot represent the direction accurately.
///
/// # Examples
///
/// ```rust
/// use delaunay::prelude::topology::spaces::SphericalPoint;
///
/// let point = SphericalPoint::<2>::try_new([3.0, 4.0, 0.0])?;
/// assert_eq!(point.ambient_dimension(), 3);
/// assert!((point.coords()[0] - 0.6).abs() < 1e-12);
/// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
/// ```
#[derive(Clone, Debug, PartialEq)]
pub struct SphericalPoint<const D: usize> {
    coords: Vec<f64>,
    radius: f64,
}

impl<const D: usize> SphericalPoint<D> {
    /// Creates a unit-radius spherical point from raw ambient coordinates.
    ///
    /// # Errors
    ///
    /// Returns [`SphericalPointError`] when `D + 1` is unrepresentable, the
    /// coordinate count is not `D + 1`, any coordinate is non-finite, or the
    /// input vector has zero norm.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let point = SphericalPoint::<2>::try_new([3.0, 4.0, 0.0])?;
    ///
    /// assert_eq!(point.ambient_dimension(), 3);
    /// assert!((point.coords()[0] - 0.6).abs() < 1e-12);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    pub fn try_new<const N: usize>(coords: [f64; N]) -> Result<Self, SphericalPointError> {
        Self::try_new_with_radius(coords, 1.0)
    }

    /// Creates a spherical point on radius `radius`.
    ///
    /// # Errors
    ///
    /// Returns [`SphericalPointError`] when `radius` is not finite and positive,
    /// `D + 1` is unrepresentable, the coordinate count is not `D + 1`, any
    /// coordinate is non-finite, or the input vector has zero norm. Returns
    /// [`SphericalPointError::UnrepresentableNormalization`] if radius scaling
    /// cannot preserve the computed unit direction within the type's tolerance.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let point = SphericalPoint::<2>::try_new_with_radius([0.0, 3.0, 4.0], 2.0)?;
    ///
    /// assert_eq!(point.radius(), 2.0);
    /// assert!((point.squared_norm() - 4.0).abs() < 1e-12);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    pub fn try_new_with_radius<const N: usize>(
        coords: [f64; N],
        radius: f64,
    ) -> Result<Self, SphericalPointError> {
        Self::try_from_slice_with_radius(&coords, radius)
    }

    /// Creates a unit-radius spherical point from a raw coordinate slice.
    ///
    /// # Errors
    ///
    /// Returns [`SphericalPointError::UnrepresentableDimension`] when `D + 1`
    /// cannot be represented as a `usize`, or
    /// [`SphericalPointError::InvalidAmbientCoordinateCount`] unless the slice
    /// has `D + 1` entries. Returns [`SphericalPointError::NonFiniteCoordinate`] if
    /// any entry is non-finite, or [`SphericalPointError::ZeroNorm`] if all entries
    /// are zero.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let raw = vec![0.0, 0.0, 2.0];
    /// let point = SphericalPoint::<2>::try_from_slice(&raw)?;
    ///
    /// assert_eq!(point.coords(), &[0.0, 0.0, 1.0]);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    pub fn try_from_slice(coords: &[f64]) -> Result<Self, SphericalPointError> {
        Self::try_from_slice_with_radius(coords, 1.0)
    }

    /// Creates a spherical point from a raw coordinate slice and radius.
    ///
    /// # Errors
    ///
    /// Returns [`SphericalPointError::InvalidRadius`] if `radius` is not finite
    /// and positive, [`SphericalPointError::InvalidAmbientCoordinateCount`]
    /// unless the slice has `D + 1` entries,
    /// [`SphericalPointError::NonFiniteCoordinate`] if any entry is non-finite,
    /// or [`SphericalPointError::ZeroNorm`] if all entries are zero.
    /// Returns [`SphericalPointError::UnrepresentableNormalization`] if scaling
    /// a finite nonzero direction to `radius` exceeds the unit-coordinate
    /// tolerance documented on [`SphericalPoint`].
    /// Returns [`SphericalPointError::UnrepresentableDimension`] when `D + 1`
    /// cannot be represented as a `usize`.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let raw = [1.0, 0.0, 0.0];
    /// let point = SphericalPoint::<2>::try_from_slice_with_radius(&raw, 3.0)?;
    ///
    /// assert_eq!(point.coords(), &[3.0, 0.0, 0.0]);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    pub fn try_from_slice_with_radius(
        coords: &[f64],
        radius: f64,
    ) -> Result<Self, SphericalPointError> {
        validate_radius(radius)?;
        let expected = Self::ambient_dimension_for_intrinsic()
            .ok_or(SphericalPointError::UnrepresentableDimension { dimension: D })?;
        if coords.len() != expected {
            return Err(SphericalPointError::InvalidAmbientCoordinateCount {
                dimension: D,
                expected,
                actual: coords.len(),
            });
        }

        let mut normalized = coords.to_vec();
        normalize_coordinates(&mut normalized, radius)?;
        Ok(Self {
            coords: normalized,
            radius,
        })
    }

    /// Returns the intrinsic sphere dimension `D`.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let point = SphericalPoint::<3>::try_new([1.0, 0.0, 0.0, 0.0])?;
    ///
    /// assert_eq!(point.intrinsic_dimension(), 3);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    #[must_use]
    pub const fn intrinsic_dimension(&self) -> usize {
        D
    }

    /// Returns the ambient coordinate dimension `D + 1`.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let point = SphericalPoint::<2>::try_new([1.0, 0.0, 0.0])?;
    ///
    /// assert_eq!(point.ambient_dimension(), 3);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    #[must_use]
    pub const fn ambient_dimension(&self) -> usize {
        self.coords.len()
    }

    /// Returns the ambient coordinate dimension for this intrinsic `D`.
    ///
    /// Returns `None` when `D + 1` is not representable as a `usize`.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// assert_eq!(SphericalPoint::<4>::ambient_dimension_for_intrinsic(), Some(5));
    /// assert_eq!(SphericalPoint::<{ usize::MAX }>::ambient_dimension_for_intrinsic(), None);
    /// ```
    #[must_use]
    pub const fn ambient_dimension_for_intrinsic() -> Option<usize> {
        D.checked_add(1)
    }

    /// Returns the radius this point was normalized to.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let point = SphericalPoint::<2>::try_new_with_radius([1.0, 0.0, 0.0], 2.5)?;
    ///
    /// assert_eq!(point.radius(), 2.5);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    #[must_use]
    pub const fn radius(&self) -> f64 {
        self.radius
    }

    /// Returns the normalized ambient coordinates.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let point = SphericalPoint::<2>::try_new([0.0, 0.0, 8.0])?;
    ///
    /// assert_eq!(point.coords(), &[0.0, 0.0, 1.0]);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    #[must_use]
    pub fn coords(&self) -> &[f64] {
        &self.coords
    }

    /// Returns the squared Euclidean norm of the stored ambient coordinates.
    ///
    /// This rounded scalar observation is in squared coordinate units. It can
    /// overflow to infinity for a large valid radius or underflow to zero for a
    /// small valid radius. It does not certify normalization or finite distances;
    /// use [`Self::radius`] for the radius established during construction.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalPoint;
    ///
    /// let point = SphericalPoint::<2>::try_new_with_radius([0.0, 1.0, 0.0], 4.0)?;
    ///
    /// assert!((point.squared_norm() - 16.0).abs() < 1e-12);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    #[must_use]
    pub fn squared_norm(&self) -> f64 {
        squared_norm_slice(&self.coords)
    }
}

/// Metric and canonicalization backend for `S^D`.
///
/// `D` is the intrinsic dimension. The backend accepts and emits
/// [`SphericalPoint<D>`] values whose ambient coordinate length is `D + 1`.
///
/// # Examples
///
/// ```rust
/// use delaunay::prelude::topology::spaces::{SphericalMetric, SphericalPoint};
/// use std::f64::consts::FRAC_PI_2;
///
/// let metric = SphericalMetric::<2>::unit();
/// let x = SphericalPoint::<2>::try_new([1.0, 0.0, 0.0])?;
/// let y = SphericalPoint::<2>::try_new([0.0, 1.0, 0.0])?;
/// assert!((metric.try_distance(&x, &y)? - FRAC_PI_2).abs() < 1e-12);
/// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
/// ```
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct SphericalMetric<const D: usize> {
    radius: f64,
}

impl<const D: usize> SphericalMetric<D> {
    /// Creates a unit-radius spherical metric.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalMetric;
    ///
    /// let metric = SphericalMetric::<2>::unit();
    ///
    /// assert_eq!(metric.radius(), 1.0);
    /// assert_eq!(metric.ambient_dimension(), Some(3));
    /// ```
    #[must_use]
    pub const fn unit() -> Self {
        Self { radius: 1.0 }
    }

    /// Creates a spherical metric with explicit radius.
    ///
    /// # Errors
    ///
    /// Returns [`SphericalPointError::InvalidRadius`] when `radius` is not
    /// finite and positive.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalMetric;
    ///
    /// let metric = SphericalMetric::<3>::try_new(2.0)?;
    ///
    /// assert_eq!(metric.radius(), 2.0);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    pub fn try_new(radius: f64) -> Result<Self, SphericalPointError> {
        validate_radius(radius)?;
        Ok(Self { radius })
    }

    /// Returns the sphere radius.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalMetric;
    ///
    /// assert_eq!(SphericalMetric::<2>::try_new(1.5)?.radius(), 1.5);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    #[must_use]
    pub const fn radius(self) -> f64 {
        self.radius
    }

    /// Returns the ambient coordinate dimension `D + 1`.
    ///
    /// Returns `None` when `D + 1` is not representable as a `usize`. A metric
    /// stores only its radius; unlike a [`SphericalPoint`], it does not prove
    /// that ambient coordinates of this dimension can be stored.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalMetric;
    ///
    /// assert_eq!(SphericalMetric::<3>::unit().ambient_dimension(), Some(4));
    /// assert_eq!(SphericalMetric::<{ usize::MAX }>::unit().ambient_dimension(), None);
    /// ```
    #[must_use]
    pub const fn ambient_dimension(self) -> Option<usize> {
        SphericalPoint::<D>::ambient_dimension_for_intrinsic()
    }

    /// Canonicalizes raw ambient coordinates onto this metric's radius.
    ///
    /// # Errors
    ///
    /// Returns [`SphericalPointError::InvalidAmbientCoordinateCount`] unless
    /// there are `D + 1` entries, [`SphericalPointError::NonFiniteCoordinate`] if
    /// any entry is non-finite, or [`SphericalPointError::ZeroNorm`] if all entries
    /// are zero. Returns [`SphericalPointError::UnrepresentableNormalization`]
    /// if this metric's radius cannot preserve the finite nonzero direction
    /// within the unit-coordinate tolerance documented on [`SphericalPoint`].
    /// Returns [`SphericalPointError::UnrepresentableDimension`] when `D + 1`
    /// cannot be represented as a `usize`.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalMetric;
    ///
    /// let metric = SphericalMetric::<2>::try_new(2.0)?;
    /// let point = metric.canonicalize([0.0, 0.0, 5.0])?;
    ///
    /// assert_eq!(point.coords(), &[0.0, 0.0, 2.0]);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    pub fn canonicalize<const N: usize>(
        self,
        coords: [f64; N],
    ) -> Result<SphericalPoint<D>, SphericalPointError> {
        SphericalPoint::try_new_with_radius(coords, self.radius)
    }

    /// Canonicalizes a raw ambient coordinate slice onto this metric's radius.
    ///
    /// # Errors
    ///
    /// Returns the unrepresentable-dimension, coordinate-count,
    /// non-finite-coordinate, zero-norm, or
    /// [`SphericalPointError::UnrepresentableNormalization`] errors described
    /// by [`Self::canonicalize`]. Even a finite nonzero input can fail when
    /// this metric's radius cannot preserve its direction within the
    /// unit-coordinate tolerance documented on [`SphericalPoint`].
    ///
    /// # Examples
    ///
    /// ```rust
    /// use delaunay::prelude::topology::spaces::SphericalMetric;
    ///
    /// let metric = SphericalMetric::<2>::unit();
    /// let raw = [3.0, 4.0, 0.0];
    /// let point = metric.canonicalize_slice(&raw)?;
    ///
    /// assert!((point.coords()[1] - 0.8).abs() < 1e-12);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    pub fn canonicalize_slice(
        self,
        coords: &[f64],
    ) -> Result<SphericalPoint<D>, SphericalPointError> {
        SphericalPoint::try_from_slice_with_radius(coords, self.radius)
    }

    /// Returns the geodesic distance between two points on this sphere.
    ///
    /// The returned arc length is in the same units as this metric's radius and
    /// lies in `[0, pi * radius]` up to floating-point roundoff.
    /// This is a rounded numerical distance, not an exact geometric predicate.
    /// Delegates the unsigned angle to la-stack's [`VectorAngle::angle`] on
    /// borrowed ambient coordinates, without allocation or a dimension cap.
    /// The exterior-product calculation takes `O((D + 1)^2)` work and constant
    /// auxiliary storage; it is intended for small ambient dimensions.
    /// Upstream power-of-two scaling and compensated exterior products avoid
    /// the endpoint sensitivity of `acos(dot)` and preserve representable
    /// subnormal angles. See the implementation provenance and background in
    /// [REFERENCES.md](https://github.com/acgetchell/delaunay/blob/main/REFERENCES.md#stable-vector-angles).
    ///
    /// # Errors
    ///
    /// Returns [`SphericalPointError::MismatchedRadius`] when either point was
    /// canonicalized for a different radius than this metric. Returns
    /// [`SphericalPointError::NonFiniteDistance`] when the arc length is not
    /// representable as a finite `f64`. An upstream numerical failure is
    /// preserved as [`SphericalPointError::AngleComputationFailed`] with its
    /// typed source.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use std::f64::consts::FRAC_PI_2;
    ///
    /// use delaunay::prelude::topology::spaces::{SphericalMetric, SphericalPoint};
    ///
    /// let metric = SphericalMetric::<2>::unit();
    /// let x = SphericalPoint::<2>::try_new([1.0, 0.0, 0.0])?;
    /// let y = SphericalPoint::<2>::try_new([0.0, 1.0, 0.0])?;
    ///
    /// assert!((metric.try_distance(&x, &y)? - FRAC_PI_2).abs() < 1e-12);
    /// # Ok::<(), delaunay::prelude::topology::spaces::SphericalPointError>(())
    /// ```
    pub fn try_distance(
        self,
        a: &SphericalPoint<D>,
        b: &SphericalPoint<D>,
    ) -> Result<f64, SphericalPointError> {
        for point in [a, b] {
            if point.radius().to_bits() != self.radius.to_bits() {
                return Err(SphericalPointError::MismatchedRadius {
                    expected: self.radius,
                    actual: point.radius(),
                });
            }
        }
        let angle = a
            .coords()
            .angle(b.coords())
            .map_err(|source| SphericalPointError::AngleComputationFailed { source })?;
        let distance = self.radius * angle;
        if distance.is_finite() {
            Ok(distance)
        } else {
            Err(SphericalPointError::NonFiniteDistance { distance })
        }
    }
}

/// Errors emitted while constructing, canonicalizing, or measuring spherical points.
#[derive(Clone, Debug, Error, PartialEq)]
#[non_exhaustive]
pub enum SphericalPointError {
    /// The intrinsic dimension's ambient arity cannot be represented by `usize`.
    #[error("spherical S^{dimension} ambient coordinate count is unrepresentable")]
    UnrepresentableDimension {
        /// Intrinsic sphere dimension whose ambient coordinate count overflows.
        dimension: usize,
    },

    /// Raw coordinate arity did not match `D + 1`.
    #[error("spherical S^{dimension} points require {expected} ambient coordinates, got {actual}")]
    InvalidAmbientCoordinateCount {
        /// Intrinsic sphere dimension.
        dimension: usize,
        /// Required ambient coordinate count.
        expected: usize,
        /// Actual coordinate count.
        actual: usize,
    },

    /// Radius was not finite and positive.
    #[error("spherical radius must be finite and positive, got {radius:?}")]
    InvalidRadius {
        /// Invalid radius.
        radius: f64,
    },

    /// A coordinate was non-finite.
    #[error("non-finite spherical coordinate {value:?} at ambient axis {axis}")]
    NonFiniteCoordinate {
        /// Ambient coordinate axis.
        axis: usize,
        /// Non-finite coordinate value.
        value: f64,
    },

    /// The raw vector had zero norm and cannot be projected onto a sphere.
    #[error("cannot normalize zero-length vector onto a sphere")]
    ZeroNorm,

    /// Scaling the coordinate vector produced a non-finite norm.
    #[error("cannot normalize vector with non-finite scaled norm {norm:?}")]
    NonFiniteNorm {
        /// Non-finite norm.
        norm: f64,
    },

    /// Radius scaling lost more than `4 * f64::EPSILON` in a unit coordinate.
    #[error(
        "spherical radius {radius:?} cannot represent unit coordinate {unit_coordinate:?} at axis {axis}: round trip produced {roundtrip_coordinate:?}"
    )]
    UnrepresentableNormalization {
        /// Ambient axis whose direction cannot be represented accurately.
        axis: usize,
        /// Requested finite positive radius.
        radius: f64,
        /// Computed coordinate before scaling to the requested radius.
        unit_coordinate: f64,
        /// Stored coordinate divided back by the requested radius.
        roundtrip_coordinate: f64,
    },

    /// Already-normalized points had different radii.
    #[error("spherical point radius {actual:?} does not match expected radius {expected:?}")]
    MismatchedRadius {
        /// Expected radius.
        expected: f64,
        /// Actual radius.
        actual: f64,
    },

    /// The upstream unsigned-angle calculation failed.
    #[error("spherical vector-angle calculation failed: {source}")]
    AngleComputationFailed {
        /// Typed numerical failure, including upstream location and provenance.
        #[source]
        source: LaError,
    },

    /// The requested geodesic distance is not representable as a finite `f64`.
    #[error("spherical geodesic distance is not finite: {distance:?}")]
    NonFiniteDistance {
        /// Non-finite computed distance.
        distance: f64,
    },
}

/// Converts validated spherical coordinates into an ambient fixed-size array.
///
/// Spherical points store their `D + 1` ambient coordinates in a vector because
/// stable const generics cannot expose `[f64; D + 1]` directly. This helper is
/// the checked bridge back to fixed-size arrays for ambient predicates and
/// convex-hull construction.
pub(crate) fn ambient_array_from_slice<const A: usize>(
    coords: &[f64],
) -> Result<[f64; A], SphericalPointError> {
    if coords.len() != A {
        return Err(SphericalPointError::InvalidAmbientCoordinateCount {
            dimension: A.saturating_sub(1),
            expected: A,
            actual: coords.len(),
        });
    }
    let mut out = [0.0; A];
    out.copy_from_slice(coords);
    Ok(out)
}

/// Validates that a radius can define a spherical metric.
///
/// Public point and metric constructors rely on this helper to keep every
/// stored sphere radius finite and strictly positive.
fn validate_radius(radius: f64) -> Result<(), SphericalPointError> {
    if radius.is_finite() && radius > 0.0 {
        return Ok(());
    }
    Err(SphericalPointError::InvalidRadius { radius })
}

/// Normalizes a finite nonzero coordinate slice onto radius `radius`.
///
/// This is the shared canonicalization path for the spherical coordinate and
/// metric backend. It classifies user-visible normalization failures while
/// scaling by the largest coordinate magnitude so finite vectors can be
/// projected without avoidable overflow or underflow.
/// Projection must accept a direction whose unscaled norm is unrepresentable.
/// Keep this radius-scaled, runtime-slice
/// reduction here: spherical storage has `D + 1` entries, which stable const
/// generics cannot pass directly to `la_stack::Vector<D>`.
fn normalize_coordinates(coords: &mut [f64], radius: f64) -> Result<(), SphericalPointError> {
    let mut max_abs = 0.0_f64;
    for (axis, coord) in coords.iter().copied().enumerate() {
        if !coord.is_finite() {
            return Err(SphericalPointError::NonFiniteCoordinate { axis, value: coord });
        }
        max_abs = max_abs.max(coord.abs());
    }
    if max_abs == 0.0 {
        return Err(SphericalPointError::ZeroNorm);
    }

    let sum_scaled_squares = coords.iter().fold(0.0, |acc, &coord| {
        let scaled = coord / max_abs;
        scaled.mul_add(scaled, acc)
    });
    let scale = sum_scaled_squares.sqrt();
    if !scale.is_finite() {
        return Err(SphericalPointError::NonFiniteNorm { norm: scale });
    }

    // Establish the output representation before mutating the slice: callers
    // of the in-place unit-sphere helper retain their input on failure. The
    // absolute unit-coordinate tolerance admits negligible component loss but
    // rejects severe subnormal distortion, including an all-zero projection.
    for (axis, &coord) in coords.iter().enumerate() {
        let unit_coordinate = (coord / max_abs) / scale;
        let roundtrip_coordinate = (unit_coordinate * radius) / radius;
        if (roundtrip_coordinate - unit_coordinate).abs() > 4.0 * f64::EPSILON {
            return Err(SphericalPointError::UnrepresentableNormalization {
                axis,
                radius,
                unit_coordinate,
                roundtrip_coordinate,
            });
        }
    }
    for coord in coords {
        *coord = ((*coord / max_abs) / scale) * radius;
    }
    Ok(())
}

/// Computes squared Euclidean norm for stored ambient coordinates.
///
/// Preserves [`SphericalPoint::squared_norm`]'s scalar observation contract over
/// `D + 1` ambient entries, including overflow to infinity and underflow to zero
/// at extreme radii. Normalization and distance certification use their own
/// scaled computations instead of this observation.
fn squared_norm_slice(coords: &[f64]) -> f64 {
    coords
        .iter()
        .fold(0.0, |acc, &coord| coord.mul_add(coord, acc))
}

/// Projects finite, nonzero coordinate vectors onto the unit sphere.
///
/// Returns `true` when the slice was normalized in place. Returns `false`
/// without modifying `coords` when any coordinate is non-finite or the vector
/// has zero length. The implementation scales by the largest absolute
/// coordinate before summing squares so very large and very small finite vectors
/// can be normalized without overflow or underflow.
///
/// This helper is crate-private glue for [`crate::topology::traits::GlobalTopology::Spherical`]
/// behavior-model paths that still operate on fixed-size coordinate arrays.
/// New spherical Delaunay construction should use [`SphericalPoint`] and
/// [`SphericalMetric`] directly so the intrinsic dimension `D` and ambient
/// dimension `D + 1` stay explicit.
pub(crate) fn normalize_unit_sphere_coordinates(coords: &mut [f64]) -> bool {
    normalize_coordinates(coords, 1.0).is_ok()
}

#[cfg(test)]
mod tests {
    use std::{
        assert_matches,
        f64::consts::{FRAC_1_SQRT_2, FRAC_PI_2},
    };

    use approx::assert_relative_eq;
    use la_stack::{NonFiniteLocation, NonFiniteOrigin, VectorOperand};

    use super::*;

    fn squared_norm(coords: &[f64]) -> f64 {
        coords
            .iter()
            .fold(0.0, |acc, &coord| coord.mul_add(coord, acc))
    }

    fn assert_unit_norm(coords: &[f64]) {
        assert_relative_eq!(squared_norm(coords), 1.0, epsilon = 1e-12);
    }

    /// Check arity publication and rejection across supported small dimensions.
    fn check_spherical_ambient_arity<const D: usize>() {
        let expected = D + 1;
        let mut coords = vec![0.0; expected];
        coords[0] = 1.0;
        let point = SphericalPoint::<D>::try_from_slice(&coords).unwrap();
        assert_eq!(point.ambient_dimension(), expected);
        assert_eq!(
            SphericalPoint::<D>::ambient_dimension_for_intrinsic(),
            Some(expected)
        );
        assert_eq!(
            SphericalMetric::<D>::unit().ambient_dimension(),
            Some(expected)
        );
        assert_eq!(
            SphericalPoint::<D>::try_from_slice(&coords[..D]),
            Err(SphericalPointError::InvalidAmbientCoordinateCount {
                dimension: D,
                expected,
                actual: D,
            })
        );
        assert_unit_norm(point.coords());
    }

    macro_rules! spherical_ambient_arity_tests {
        ($($dimension:literal),+ $(,)?) => {
            pastey::paste! {
                $(
                    #[test]
                    fn [<spherical_ambient_arity_ $dimension d>]() {
                        check_spherical_ambient_arity::<$dimension>();
                    }
                )+
            }
        };
    }

    spherical_ambient_arity_tests!(2, 3, 4, 5);

    #[test]
    fn spherical_coordinates_reject_unrepresentable_dimension() {
        let error = SphericalPointError::UnrepresentableDimension {
            dimension: usize::MAX,
        };
        assert_eq!(
            SphericalPoint::<{ usize::MAX }>::ambient_dimension_for_intrinsic(),
            None
        );
        assert_eq!(
            SphericalPoint::<{ usize::MAX }>::try_from_slice(&[]),
            Err(error.clone())
        );
        assert_eq!(
            SphericalPoint::<{ usize::MAX }>::try_new([]),
            Err(error.clone())
        );
        let metric = SphericalMetric::<{ usize::MAX }>::unit();
        assert_eq!(metric.ambient_dimension(), None);
        assert_eq!(metric.canonicalize_slice(&[]), Err(error));
    }

    #[test]
    fn spherical_point_handles_near_zero_vector() {
        let point = SphericalPoint::<2>::try_new([f64::MIN_POSITIVE, f64::MIN_POSITIVE, 0.0])
            .expect("finite nonzero ambient vector should normalize");
        let coords = point.coords();
        assert_relative_eq!(coords[0], FRAC_1_SQRT_2);
        assert_relative_eq!(coords[1], FRAC_1_SQRT_2);
        assert_relative_eq!(coords[2], 0.0);
        assert_unit_norm(coords);
    }

    #[test]
    fn spherical_point_handles_large_vector_without_overflow() {
        let point = SphericalPoint::<2>::try_new([f64::MAX, -f64::MAX, 0.0])
            .expect("finite nonzero ambient vector should normalize");
        let coords = point.coords();
        assert_relative_eq!(coords[0], FRAC_1_SQRT_2);
        assert_relative_eq!(coords[1], -FRAC_1_SQRT_2);
        assert_relative_eq!(coords[2], 0.0);
        assert_unit_norm(coords);
    }

    #[test]
    fn spherical_distance_uses_scaled_coordinates_for_large_radius() {
        let radius = f64::MAX / 4.0;
        let metric = SphericalMetric::<2>::try_new(radius)
            .expect("large finite radius should define a metric");
        let x = SphericalPoint::<2>::try_new_with_radius([1.0, 0.0, 0.0], radius)
            .expect("axis point should normalize onto large radius");
        let diagonal = SphericalPoint::<2>::try_new_with_radius([1.0, 1.0, 0.0], radius)
            .expect("diagonal point should normalize onto large radius");

        let distance = metric
            .try_distance(&x, &diagonal)
            .expect("large-radius dot product should avoid overflow");

        assert_relative_eq!(distance / radius, std::f64::consts::FRAC_PI_4);
    }

    #[test]
    fn spherical_distance_handles_subnormal_radius_without_reciprocal_overflow() {
        let radius = f64::MIN_POSITIVE / 2.0;
        let metric = SphericalMetric::<2>::try_new(radius)
            .expect("positive subnormal radius should define a metric");
        let x = SphericalPoint::<2>::try_new_with_radius([1.0, 0.0, 0.0], radius)
            .expect("axis point should normalize onto a tiny sphere");
        let y = SphericalPoint::<2>::try_new_with_radius([0.0, 1.0, 0.0], radius)
            .expect("orthogonal point should normalize onto a tiny sphere");

        let distance = metric
            .try_distance(&x, &y)
            .expect("tiny finite arc length should remain representable");

        assert_relative_eq!(distance / radius, FRAC_PI_2);
    }

    #[test]
    fn spherical_point_rejects_unrepresentable_subnormal_directions() {
        fn check<const D: usize>() {
            let mut coords = vec![1.0; D + 1];
            coords[0] = 2.0;
            for radius in [f64::from_bits(1), f64::from_bits(2)] {
                assert_matches!(
                    SphericalPoint::<D>::try_from_slice_with_radius(&coords, radius),
                    Err(SphericalPointError::UnrepresentableNormalization { axis: 0, .. })
                );
                let mut unchanged = coords.clone();
                assert!(normalize_coordinates(&mut unchanged, radius).is_err());
                assert_eq!(unchanged, coords);
            }
        }
        check::<2>();
        check::<3>();
        check::<4>();
        check::<5>();
        assert_matches!(
            SphericalPoint::<4>::try_new_with_radius([1.0; 5], f64::from_bits(1)),
            Err(SphericalPointError::UnrepresentableNormalization { .. })
        );
    }

    #[test]
    fn spherical_normalization_preserves_input_on_late_rejection() {
        let original = [1.0, 1.0, 2.0_f64.sqrt()];
        let mut coords = original;

        // The first two unit coordinates survive scaling to this radius, but
        // the last does not. Rejection must precede every in-place write.
        let error = normalize_coordinates(&mut coords, f64::from_bits(2)).unwrap_err();
        assert_matches!(
            error,
            SphericalPointError::UnrepresentableNormalization { axis: 2, .. }
        );
        assert_eq!(coords.map(f64::to_bits), original.map(f64::to_bits));
    }

    #[test]
    fn spherical_point_preserves_representable_subnormal_directions() {
        // Axis directions remain exactly representable even at the smallest radius.
        let radius = f64::from_bits(1);
        let point = SphericalPoint::<4>::try_new_with_radius([1.0, 0.0, 0.0, 0.0, 0.0], radius)
            .expect("the smallest positive axis coordinate is representable");
        assert_eq!(point.coords()[0].to_bits(), radius.to_bits());
        let metric = SphericalMetric::<4>::try_new(radius).unwrap();
        assert_eq!(metric.try_distance(&point, &point).unwrap().to_bits(), 0);

        // Four equal components also remain exact at a radius of two subnormal units.
        let diagonal = SphericalPoint::<3>::try_new_with_radius([1.0; 4], f64::from_bits(2))
            .expect("exactly representable subnormal diagonal should succeed");
        assert!(diagonal.coords().iter().all(|coord| coord.to_bits() == 1));

        // A large subnormal radius still resolves a general direction accurately.
        let radius = f64::MIN_POSITIVE / 2.0;
        let diagonal = SphericalPoint::<2>::try_new_with_radius([1.0, 1.0, 0.0], radius)
            .expect("subnormal scaling with sufficient precision should succeed");
        assert_relative_eq!(diagonal.coords()[0] / radius, FRAC_1_SQRT_2);

        // Losing an insignificant component does not invalidate the direction.
        let near_axis =
            SphericalPoint::<2>::try_new_with_radius([1.0, f64::MIN_POSITIVE, 0.0], radius)
                .expect("negligible component underflow stays within the direction tolerance");
        assert_eq!(near_axis.coords()[0].to_bits(), radius.to_bits());
    }

    #[test]
    fn spherical_distance_rejects_unrepresentable_arc_length() {
        let metric = SphericalMetric::<2>::try_new(f64::MAX)
            .expect("maximum finite radius should define a metric");
        let x = SphericalPoint::<2>::try_new_with_radius([1.0, 0.0, 0.0], f64::MAX)
            .expect("axis point should normalize onto maximum radius");
        let opposite = SphericalPoint::<2>::try_new_with_radius([-1.0, 0.0, 0.0], f64::MAX)
            .expect("opposite point should normalize onto maximum radius");

        assert_matches!(
            metric.try_distance(&x, &opposite),
            Err(SphericalPointError::NonFiniteDistance { distance }) if distance.is_infinite()
        );
    }

    #[test]
    fn spherical_distance_preserves_upstream_error_provenance() {
        let metric = SphericalMetric::<2>::unit();
        let valid = SphericalPoint::<2>::try_new([1.0, 0.0, 0.0]).unwrap();
        // Bypass the public constructor only to exercise defensive propagation;
        // ordinary SphericalPoint values already prove finite nonzero inputs.
        let corrupted = SphericalPoint::<2> {
            coords: vec![1.0, f64::NAN, 0.0],
            radius: 1.0,
        };
        assert_matches!(
            metric.try_distance(&valid, &corrupted),
            Err(SphericalPointError::AngleComputationFailed {
                source: LaError::NonFinite {
                    location: NonFiniteLocation::VectorOperandEntry {
                        operand: VectorOperand::Right,
                        index: 1,
                        ..
                    },
                    origin: NonFiniteOrigin::Input,
                    ..
                },
            })
        );
    }

    #[test]
    fn ambient_array_from_slice_checks_fixed_ambient_arity() {
        let coords = ambient_array_from_slice::<3>(&[1.0, 2.0, 3.0])
            .expect("matching ambient arity should copy into a fixed array");
        assert_relative_eq!(coords[0], 1.0);
        assert_relative_eq!(coords[1], 2.0);
        assert_relative_eq!(coords[2], 3.0);

        let err = ambient_array_from_slice::<3>(&[1.0, 2.0])
            .expect_err("wrong ambient arity should remain typed");
        assert_eq!(
            err,
            SphericalPointError::InvalidAmbientCoordinateCount {
                dimension: 2,
                expected: 3,
                actual: 2,
            }
        );
    }

    #[test]
    fn normalize_unit_sphere_coordinates_preserves_failed_inputs() {
        let mut zero = [0.0, 0.0, 0.0];
        assert!(!normalize_unit_sphere_coordinates(&mut zero));
        assert_relative_eq!(zero[0], 0.0);
        assert_relative_eq!(zero[1], 0.0);
        assert_relative_eq!(zero[2], 0.0);

        let mut coords = [1.0, f64::INFINITY, 0.0];
        assert!(!normalize_unit_sphere_coordinates(&mut coords));
        assert_relative_eq!(coords[0], 1.0);
        assert!(coords[1].is_infinite());
        assert_relative_eq!(coords[2], 0.0);

        let mut coords = [1.0, f64::NAN, -2.0];
        assert!(!normalize_unit_sphere_coordinates(&mut coords));
        assert_relative_eq!(coords[0], 1.0);
        assert!(coords[1].is_nan());
        assert_relative_eq!(coords[2], -2.0);
    }
}