oxiproj-core 0.1.2

//! The `Scalar` trait — a numeric surface for generic projection math.
//!
//! Projection structs that implement [`crate::projection_generic::ProjectGeneric`]
//! write their forward/inverse formulae once, parameterised over `S: Scalar`.
//! At `S = f64` you get the ordinary numeric result; at `S = Dual1<2>` you
//! get the full 2×2 Jacobian for free via forward-mode AD; at `S = Dual2<2>`
//! you also get the diagonal Hessian.
//!
//! # Implementing `Scalar`
//!
//! The trait is sealed only by its many supertraits. Any type that satisfies
//! the arithmetic, ordering, and conversion bounds and provides the listed
//! transcendental operations is a valid scalar.
//!
//! Implementations for `f64`, [`Dual1<N>`](crate::autodiff::Dual1), and
//! [`Dual2<N>`](crate::autodiff::Dual2) are provided here.

use crate::autodiff::{Dual1, Dual2};
use core::ops::{Add, Div, Mul, Neg, Sub};

/// A numeric surface for generic projection computations.
///
/// Implementations must provide all standard arithmetic operators,
/// the common transcendental functions, and two conversion helpers.
/// Constants are provided as default methods backed by `from_f64`.
pub trait Scalar:
    Copy
    + Clone
    + core::fmt::Debug
    + core::fmt::Display
    + Add<Output = Self>
    + Sub<Output = Self>
    + Mul<Output = Self>
    + Div<Output = Self>
    + Neg<Output = Self>
    + PartialOrd
    + PartialEq
    + 'static
    + Send
    + Sync
{
    /// Additive identity.
    fn zero() -> Self;
    /// Multiplicative identity.
    fn one() -> Self;
    /// `one() + one()`
    fn two() -> Self {
        Self::one() + Self::one()
    }

    /// Lossless conversion from a compile-time constant.
    fn from_f64(v: f64) -> Self;
    /// Lossy projection to `f64`. For dual numbers, returns the primal value.
    fn to_f64(self) -> f64;

    // ---- Transcendentals ----

    /// Sine.
    fn sin(self) -> Self;
    /// Cosine.
    fn cos(self) -> Self;
    /// Tangent.
    fn tan(self) -> Self;
    /// Arcsine (inverse sine).
    fn asin(self) -> Self;
    /// Arccosine (inverse cosine).
    fn acos(self) -> Self;
    /// Arctangent (inverse tangent).
    fn atan(self) -> Self;
    /// Four-quadrant arctangent of `self / other`.
    fn atan2(self, other: Self) -> Self;
    /// Square root.
    fn sqrt(self) -> Self;
    /// Cube root.
    fn cbrt(self) -> Self;
    /// Exponential function, `e^self`.
    fn exp(self) -> Self;
    /// Natural logarithm.
    fn ln(self) -> Self;
    /// Absolute value.
    fn abs(self) -> Self;
    /// Raises `self` to an integer power.
    fn powi(self, n: i32) -> Self;
    /// Raises `self` to a floating-point power.
    fn powf(self, e: f64) -> Self;
    /// Hyperbolic sine.
    fn sinh(self) -> Self;
    /// Hyperbolic cosine.
    fn cosh(self) -> Self;
    /// Hyperbolic tangent.
    fn tanh(self) -> Self;
    /// Euclidean distance, `sqrt(self^2 + other^2)`, computed without
    /// intermediate overflow/underflow.
    fn hypot(self, other: Self) -> Self;

    // ---- Geodetic constants (default impls via `from_f64`) ----

    /// The constant pi (`π`).
    fn pi() -> Self {
        Self::from_f64(core::f64::consts::PI)
    }
    /// Two pi (`2π`, a full turn in radians).
    fn two_pi() -> Self {
        Self::from_f64(core::f64::consts::TAU)
    }
    /// Half pi (`π/2`, a quarter turn in radians).
    fn half_pi() -> Self {
        Self::from_f64(core::f64::consts::FRAC_PI_2)
    }
    /// Quarter pi (`π/4`, an eighth turn in radians).
    fn quarter_pi() -> Self {
        Self::from_f64(core::f64::consts::FRAC_PI_4)
    }
    /// Conversion factor from radians to degrees (`180/π`).
    fn rad_to_deg() -> Self {
        Self::from_f64(180.0 / core::f64::consts::PI)
    }
    /// Conversion factor from degrees to radians (`π/180`).
    fn deg_to_rad() -> Self {
        Self::from_f64(core::f64::consts::PI / 180.0)
    }
}

// ---- impl for f64 ----

impl Scalar for f64 {
    #[inline]
    fn zero() -> Self {
        0.0
    }
    #[inline]
    fn one() -> Self {
        1.0
    }
    #[inline]
    fn from_f64(v: f64) -> Self {
        v
    }
    #[inline]
    fn to_f64(self) -> f64 {
        self
    }
    #[inline]
    fn sin(self) -> Self {
        f64::sin(self)
    }
    #[inline]
    fn cos(self) -> Self {
        f64::cos(self)
    }
    #[inline]
    fn tan(self) -> Self {
        f64::tan(self)
    }
    #[inline]
    fn asin(self) -> Self {
        f64::asin(self)
    }
    #[inline]
    fn acos(self) -> Self {
        f64::acos(self)
    }
    #[inline]
    fn atan(self) -> Self {
        f64::atan(self)
    }
    #[inline]
    fn atan2(self, other: Self) -> Self {
        f64::atan2(self, other)
    }
    #[inline]
    fn sqrt(self) -> Self {
        f64::sqrt(self)
    }
    #[inline]
    fn cbrt(self) -> Self {
        f64::cbrt(self)
    }
    #[inline]
    fn exp(self) -> Self {
        f64::exp(self)
    }
    #[inline]
    fn ln(self) -> Self {
        f64::ln(self)
    }
    #[inline]
    fn abs(self) -> Self {
        f64::abs(self)
    }
    #[inline]
    fn powi(self, n: i32) -> Self {
        f64::powi(self, n)
    }
    #[inline]
    fn powf(self, e: f64) -> Self {
        f64::powf(self, e)
    }
    #[inline]
    fn sinh(self) -> Self {
        f64::sinh(self)
    }
    #[inline]
    fn cosh(self) -> Self {
        f64::cosh(self)
    }
    #[inline]
    fn tanh(self) -> Self {
        f64::tanh(self)
    }
    #[inline]
    fn hypot(self, other: Self) -> Self {
        f64::hypot(self, other)
    }
}

// ---- impl for Dual1<N> ----

impl<const N: usize> Scalar for Dual1<N> {
    #[inline]
    fn zero() -> Self {
        Self::constant(0.0)
    }
    #[inline]
    fn one() -> Self {
        Self::constant(1.0)
    }
    #[inline]
    fn from_f64(v: f64) -> Self {
        Self::constant(v)
    }
    #[inline]
    fn to_f64(self) -> f64 {
        self.v
    }
    #[inline]
    fn sin(self) -> Self {
        Dual1::sin(self)
    }
    #[inline]
    fn cos(self) -> Self {
        Dual1::cos(self)
    }
    #[inline]
    fn tan(self) -> Self {
        Dual1::tan(self)
    }
    #[inline]
    fn asin(self) -> Self {
        Dual1::asin(self)
    }
    #[inline]
    fn acos(self) -> Self {
        Dual1::acos(self)
    }
    #[inline]
    fn atan(self) -> Self {
        Dual1::atan(self)
    }
    #[inline]
    fn atan2(self, other: Self) -> Self {
        Dual1::atan2(self, other)
    }
    #[inline]
    fn sqrt(self) -> Self {
        Dual1::sqrt(self)
    }
    #[inline]
    fn cbrt(self) -> Self {
        Dual1::cbrt(self)
    }
    #[inline]
    fn exp(self) -> Self {
        Dual1::exp(self)
    }
    #[inline]
    fn ln(self) -> Self {
        Dual1::ln(self)
    }
    #[inline]
    fn abs(self) -> Self {
        Dual1::abs(self)
    }
    #[inline]
    fn powi(self, n: i32) -> Self {
        Dual1::powi(self, n)
    }
    #[inline]
    fn powf(self, e: f64) -> Self {
        Dual1::powf(self, e)
    }
    #[inline]
    fn sinh(self) -> Self {
        Dual1::sinh(self)
    }
    #[inline]
    fn cosh(self) -> Self {
        Dual1::cosh(self)
    }
    #[inline]
    fn tanh(self) -> Self {
        Dual1::tanh(self)
    }
    #[inline]
    fn hypot(self, other: Self) -> Self {
        Dual1::hypot(self, other)
    }
}

// ---- impl for Dual2<N> ----

impl<const N: usize> Scalar for Dual2<N> {
    #[inline]
    fn zero() -> Self {
        Self::constant(0.0)
    }
    #[inline]
    fn one() -> Self {
        Self::constant(1.0)
    }
    #[inline]
    fn from_f64(v: f64) -> Self {
        Self::constant(v)
    }
    #[inline]
    fn to_f64(self) -> f64 {
        self.v
    }
    #[inline]
    fn sin(self) -> Self {
        Dual2::sin(self)
    }
    #[inline]
    fn cos(self) -> Self {
        Dual2::cos(self)
    }
    #[inline]
    fn tan(self) -> Self {
        Dual2::tan(self)
    }
    #[inline]
    fn asin(self) -> Self {
        Dual2::asin(self)
    }
    #[inline]
    fn acos(self) -> Self {
        Dual2::acos(self)
    }
    #[inline]
    fn atan(self) -> Self {
        Dual2::atan(self)
    }
    #[inline]
    fn atan2(self, other: Self) -> Self {
        Dual2::atan2(self, other)
    }
    #[inline]
    fn sqrt(self) -> Self {
        Dual2::sqrt(self)
    }
    #[inline]
    fn cbrt(self) -> Self {
        Dual2::cbrt(self)
    }
    #[inline]
    fn exp(self) -> Self {
        Dual2::exp(self)
    }
    #[inline]
    fn ln(self) -> Self {
        Dual2::ln(self)
    }
    #[inline]
    fn abs(self) -> Self {
        Dual2::abs(self)
    }
    #[inline]
    fn powi(self, n: i32) -> Self {
        Dual2::powi(self, n)
    }
    #[inline]
    fn powf(self, e: f64) -> Self {
        Dual2::powf(self, e)
    }
    #[inline]
    fn sinh(self) -> Self {
        Dual2::sinh(self)
    }
    #[inline]
    fn cosh(self) -> Self {
        Dual2::cosh(self)
    }
    #[inline]
    fn tanh(self) -> Self {
        Dual2::tanh(self)
    }
    #[inline]
    fn hypot(self, other: Self) -> Self {
        Dual2::hypot(self, other)
    }
}

// ---- Tests ----

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

    fn check_scalar_mercator<S: Scalar>(lam_val: f64, phi_val: f64) -> (f64, f64) {
        let lam = S::from_f64(lam_val);
        let phi = S::from_f64(phi_val);
        let x = lam;
        let pi4 = S::quarter_pi();
        let half = S::from_f64(0.5);
        let y = (pi4 + phi * half).tan().ln();
        (x.to_f64(), y.to_f64())
    }

    #[test]
    fn scalar_f64_mercator_equator() {
        let (x, y) = check_scalar_mercator::<f64>(0.5, 0.0);
        assert!((x - 0.5).abs() < 1e-12);
        assert!(y.abs() < 1e-12);
    }

    #[test]
    fn scalar_dual1_mercator_matches_f64() {
        let phi = 0.3_f64;
        let lam = 1.2_f64;
        let (xf, yf) = check_scalar_mercator::<f64>(lam, phi);
        let (xd, yd) = check_scalar_mercator::<Dual1<2>>(lam, phi);
        assert!((xf - xd).abs() < 1e-12);
        assert!((yf - yd).abs() < 1e-12);
    }

    #[test]
    fn scalar_dual2_mercator_matches_f64() {
        let phi = 0.3_f64;
        let lam = 1.2_f64;
        let (xf, yf) = check_scalar_mercator::<f64>(lam, phi);
        let (xd, yd) = check_scalar_mercator::<Dual2<2>>(lam, phi);
        assert!((xf - xd).abs() < 1e-12);
        assert!((yf - yd).abs() < 1e-12);
    }

    #[test]
    fn scalar_constants() {
        use core::f64::consts::PI;
        assert!((f64::pi() - PI).abs() < 1e-15);
        assert!((f64::two_pi() - 2.0 * PI).abs() < 1e-15);
        assert!((f64::half_pi() - PI / 2.0).abs() < 1e-15);
        assert!((f64::quarter_pi() - PI / 4.0).abs() < 1e-15);
        assert!((f64::rad_to_deg() - 180.0 / PI).abs() < 1e-12);
        assert!((f64::deg_to_rad() - PI / 180.0).abs() < 1e-12);
    }

    #[test]
    fn scalar_zero_one_two() {
        assert!((f64::zero() - 0.0).abs() < 1e-15);
        assert!((f64::one() - 1.0).abs() < 1e-15);
        assert!((f64::two() - 2.0).abs() < 1e-15);
    }
}