xad-rs 0.8.1

//! `Jet2` - dedicated second-order forward-mode AD type.
//!
//! A `Jet2<T>` carries three values:
//!   - `value`: f(x)
//!   - `d1`:    f'(x)  (first derivative w.r.t. a seed direction)
//!   - `d2`:    f''(x) (second derivative w.r.t. the same seed direction)
//!
//! Unlike full Hessian machinery, `Jet2` tracks derivatives with respect to
//! a *single* scalar variable. This is exactly what you need to compute
//! diagonal elements of a Hessian (e.g. "own-gamma" in financial risk) in a
//! single forward pass, with no tape and no finite-difference error.
//!
//! Usage:
//! ```
//! use xad_rs::forward::jet2::Jet2;
//! // Compute f(x) = x^3 and its first/second derivatives at x = 2.
//! let x = Jet2::variable(2.0_f64);
//! let y = x * x * x;
//! assert_eq!(y.value(), 8.0);
//! assert_eq!(y.first_derivative(), 12.0); // 3x^2
//! assert_eq!(y.second_derivative(), 12.0); // 6x
//! ```
//!
//! This module also contains [`NamedJet2<T>`] — a named wrapper over
//! seeded `Jet2<T>`.
//!
//! Does NOT carry an `Arc<VarRegistry>` field in
//! release builds. The struct layout is `{ inner: Jet2<T>, seeded: Option<usize> }`
//! plus, under `#[cfg(debug_assertions)]` only, a `gen_id: u64` stamped by the
//! owning [`NamedForwardTape`] scope for the cross-registry debug guard.
//! Release builds carry zero atomic-refcount cost per operator.
//!
//! `Jet2<T>` tracks first AND second derivatives along ONE seed direction.
//! The named form gives that direction a name at construction time and
//! exposes [`first_derivative`](NamedJet2::first_derivative) /
//! [`second_derivative`](NamedJet2::second_derivative) accessors that
//! return the seeded values when the name matches and `T::zero()`
//! otherwise. The `seeded: Option<usize>` field is NOT the registry — it
//! is the positional index of the currently-seeded variable, a per-value
//! attribute that survives the minimal struct layout.
//!
//! Unlike `NamedJet1Vec` and `NamedJet1`, every binary op on
//! `NamedJet2` must propagate the `seeded: Option<usize>` field through
//! a private [`merge_seeded`] combinator. Operations between two
//! `NamedJet2<T>` values with different seeds panic in debug builds
//! (two-direction `Jet2` is a semantic violation).
//!
//! The only way to obtain a `NamedJet2<T>` is via the
//! `NamedForwardTape` constructor API (see module docs in
//! `src/forward_tape.rs`).

use crate::passive::Passive;
use std::fmt;
use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};

// ============================================================================
// Jet2<T> — positional second-order forward-mode dual number
// ============================================================================

/// Second-order forward-mode dual number.
#[derive(Clone, Copy)]
pub struct Jet2<T: Passive> {
    value: T,
    d1: T,
    d2: T,
}

impl<T: Passive> Jet2<T> {
    /// Create a `Jet2` with explicit value, first and second derivative.
    #[inline]
    pub fn new(value: T, d1: T, d2: T) -> Self {
        Jet2 { value, d1, d2 }
    }

    /// Create a constant (derivative-free) `Jet2`.
    #[inline]
    pub fn constant(value: T) -> Self {
        Jet2 {
            value,
            d1: T::zero(),
            d2: T::zero(),
        }
    }

    /// Create the active variable: value with `d1 = 1`, `d2 = 0`.
    /// Use this to seed the single input direction to differentiate against.
    #[inline]
    pub fn variable(value: T) -> Self {
        Jet2 {
            value,
            d1: T::one(),
            d2: T::zero(),
        }
    }

    /// Underlying value.
    #[inline]
    pub fn value(&self) -> T {
        self.value
    }

    /// First derivative (tangent) along the seeded direction.
    #[inline]
    pub fn first_derivative(&self) -> T {
        self.d1
    }

    /// Second derivative along the seeded direction.
    #[inline]
    pub fn second_derivative(&self) -> T {
        self.d2
    }

    /// Raise `self` to a scalar power `n`: `self^n`.
    ///
    /// Applies the chain rule for 2nd order: for `g(u) = u^n`,
    /// `g'(u) = n u^{n-1}`, `g''(u) = n (n-1) u^{n-2}`, then
    /// `result.d1 = g'(v) * self.d1` and
    /// `result.d2 = g''(v) * self.d1^2 + g'(v) * self.d2`.
    pub fn powf(self, n: T) -> Jet2<T> {
        let v = self.value;
        let two = T::from(2.0).unwrap();
        let vn = v.powf(n);
        let gp = n * v.powf(n - T::one());
        let gpp = n * (n - T::one()) * v.powf(n - two);
        Jet2 {
            value: vn,
            d1: gp * self.d1,
            d2: gpp * self.d1 * self.d1 + gp * self.d2,
        }
    }

    /// Exponential `exp(self)`.
    pub fn exp(self) -> Jet2<T> {
        let e = self.value.exp();
        Jet2 {
            value: e,
            d1: e * self.d1,
            d2: e * self.d1 * self.d1 + e * self.d2,
        }
    }

    /// Natural logarithm `ln(self)`.
    pub fn ln(self) -> Jet2<T> {
        // g(u) = ln(u), g'(u) = 1/u, g''(u) = -1/u^2
        let v = self.value;
        let inv = T::one() / v;
        let inv2 = inv * inv;
        Jet2 {
            value: v.ln(),
            d1: inv * self.d1,
            d2: -inv2 * self.d1 * self.d1 + inv * self.d2,
        }
    }

    /// Square root.
    pub fn sqrt(self) -> Jet2<T> {
        // g(u) = u^{1/2}, g'(u) = 1/(2*sqrt(u)), g''(u) = -1/(4*u^{3/2})
        let v = self.value;
        let s = v.sqrt();
        let half = T::from(0.5).unwrap();
        let quarter = T::from(0.25).unwrap();
        let gp = half / s;
        let gpp = -quarter / (s * v);
        Jet2 {
            value: s,
            d1: gp * self.d1,
            d2: gpp * self.d1 * self.d1 + gp * self.d2,
        }
    }

    /// Sine.
    pub fn sin(self) -> Jet2<T> {
        // g(u) = sin(u), g'(u) = cos(u), g''(u) = -sin(u)
        let v = self.value;
        let (s, c) = (v.sin(), v.cos());
        Jet2 {
            value: s,
            d1: c * self.d1,
            d2: -s * self.d1 * self.d1 + c * self.d2,
        }
    }

    /// Cosine.
    pub fn cos(self) -> Jet2<T> {
        // g(u) = cos(u), g'(u) = -sin(u), g''(u) = -cos(u)
        let v = self.value;
        let (s, c) = (v.sin(), v.cos());
        Jet2 {
            value: c,
            d1: -s * self.d1,
            d2: -c * self.d1 * self.d1 - s * self.d2,
        }
    }

    /// Error function `erf(self)`.
    ///
    /// Derivatives:
    ///   `erf'(x)  = (2/√π) · exp(-x²)`
    ///   `erf''(x) = -2x · erf'(x)`
    pub fn erf(self) -> Jet2<T> {
        let v = self.value;
        let r = crate::math::erf(v);
        let two = T::from(2.0).unwrap();
        // Convert the compile-time `f64` constant `FRAC_2_SQRT_PI` into the
        // generic scalar `T` exactly once — cheaper than computing a sqrt.
        let two_over_sqrt_pi =
            T::from(std::f64::consts::FRAC_2_SQRT_PI).unwrap();
        let gp = two_over_sqrt_pi * (-v * v).exp();
        let gpp = -two * v * gp;
        Jet2 {
            value: r,
            d1: gp * self.d1,
            d2: gpp * self.d1 * self.d1 + gp * self.d2,
        }
    }

    /// Standard normal CDF: `Φ(x) = 0.5 · (1 + erf(x / √2))`.
    ///
    /// - `g'(x) = φ(x) = (1/√(2π)) · exp(-x²/2)`
    /// - `g''(x) = -x · φ(x)`
    pub fn norm_cdf(self) -> Jet2<T> {
        let v = self.value;
        let r = crate::math::norm_cdf(v);
        let gp = crate::math::norm_pdf(v);
        let gpp = -v * gp;
        Jet2 {
            value: r,
            d1: gp * self.d1,
            d2: gpp * self.d1 * self.d1 + gp * self.d2,
        }
    }

    /// Inverse standard normal CDF: `Φ⁻¹(p)`.
    ///
    /// - `g'(p) = 1 / φ(Φ⁻¹(p))`
    /// - `g''(p) = Φ⁻¹(p) · g'(p)²`
    pub fn inv_norm_cdf(self) -> Jet2<T> {
        let v = self.value;
        let r = crate::math::inv_norm_cdf(v);
        let gp = T::one() / crate::math::norm_pdf(r);
        let gpp = r * gp * gp;
        Jet2 {
            value: r,
            d1: gp * self.d1,
            d2: gpp * self.d1 * self.d1 + gp * self.d2,
        }
    }
}

// ============================================================================
// Operator implementations for Jet2<T>
// ============================================================================

// --- Add ---
impl<T: Passive> Add for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn add(self, rhs: Self) -> Self {
        Jet2 {
            value: self.value + rhs.value,
            d1: self.d1 + rhs.d1,
            d2: self.d2 + rhs.d2,
        }
    }
}

impl<T: Passive> Add<T> for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn add(self, rhs: T) -> Self {
        Jet2 {
            value: self.value + rhs,
            d1: self.d1,
            d2: self.d2,
        }
    }
}

impl Add<Jet2<f64>> for f64 {
    type Output = Jet2<f64>;
    #[inline]
    fn add(self, rhs: Jet2<f64>) -> Jet2<f64> {
        Jet2 {
            value: self + rhs.value,
            d1: rhs.d1,
            d2: rhs.d2,
        }
    }
}

impl Add<Jet2<f32>> for f32 {
    type Output = Jet2<f32>;
    #[inline]
    fn add(self, rhs: Jet2<f32>) -> Jet2<f32> {
        Jet2 {
            value: self + rhs.value,
            d1: rhs.d1,
            d2: rhs.d2,
        }
    }
}

// --- Sub ---
impl<T: Passive> Sub for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn sub(self, rhs: Self) -> Self {
        Jet2 {
            value: self.value - rhs.value,
            d1: self.d1 - rhs.d1,
            d2: self.d2 - rhs.d2,
        }
    }
}

impl<T: Passive> Sub<T> for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn sub(self, rhs: T) -> Self {
        Jet2 {
            value: self.value - rhs,
            d1: self.d1,
            d2: self.d2,
        }
    }
}

impl Sub<Jet2<f64>> for f64 {
    type Output = Jet2<f64>;
    #[inline]
    fn sub(self, rhs: Jet2<f64>) -> Jet2<f64> {
        Jet2 {
            value: self - rhs.value,
            d1: -rhs.d1,
            d2: -rhs.d2,
        }
    }
}

impl Sub<Jet2<f32>> for f32 {
    type Output = Jet2<f32>;
    #[inline]
    fn sub(self, rhs: Jet2<f32>) -> Jet2<f32> {
        Jet2 {
            value: self - rhs.value,
            d1: -rhs.d1,
            d2: -rhs.d2,
        }
    }
}

// --- Mul ---
// (a*b)' = a'*b + a*b'
// (a*b)'' = a''*b + 2*a'*b' + a*b''
impl<T: Passive> Mul for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn mul(self, rhs: Self) -> Self {
        let two = T::from(2.0).unwrap();
        Jet2 {
            value: self.value * rhs.value,
            d1: self.d1 * rhs.value + self.value * rhs.d1,
            d2: self.d2 * rhs.value + two * self.d1 * rhs.d1 + self.value * rhs.d2,
        }
    }
}

impl<T: Passive> Mul<T> for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn mul(self, rhs: T) -> Self {
        Jet2 {
            value: self.value * rhs,
            d1: self.d1 * rhs,
            d2: self.d2 * rhs,
        }
    }
}

impl Mul<Jet2<f64>> for f64 {
    type Output = Jet2<f64>;
    #[inline]
    fn mul(self, rhs: Jet2<f64>) -> Jet2<f64> {
        Jet2 {
            value: self * rhs.value,
            d1: self * rhs.d1,
            d2: self * rhs.d2,
        }
    }
}

impl Mul<Jet2<f32>> for f32 {
    type Output = Jet2<f32>;
    #[inline]
    fn mul(self, rhs: Jet2<f32>) -> Jet2<f32> {
        Jet2 {
            value: self * rhs.value,
            d1: self * rhs.d1,
            d2: self * rhs.d2,
        }
    }
}

// --- Div ---
// a/b = a * (1/b)
// (1/b)'  = -b'/b^2
// (1/b)'' = 2 b'^2/b^3 - b''/b^2
impl<T: Passive> Div for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn div(self, rhs: Self) -> Self {
        let two = T::from(2.0).unwrap();
        let inv_b = T::one() / rhs.value;
        let inv_b2 = inv_b * inv_b;
        let inv_b3 = inv_b2 * inv_b;
        let recip = Jet2 {
            value: inv_b,
            d1: -rhs.d1 * inv_b2,
            d2: two * rhs.d1 * rhs.d1 * inv_b3 - rhs.d2 * inv_b2,
        };
        self * recip
    }
}

impl<T: Passive> Div<T> for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn div(self, rhs: T) -> Self {
        let inv = T::one() / rhs;
        Jet2 {
            value: self.value * inv,
            d1: self.d1 * inv,
            d2: self.d2 * inv,
        }
    }
}

impl Div<Jet2<f64>> for f64 {
    type Output = Jet2<f64>;
    #[inline]
    fn div(self, rhs: Jet2<f64>) -> Jet2<f64> {
        Jet2::constant(self) / rhs
    }
}

impl Div<Jet2<f32>> for f32 {
    type Output = Jet2<f32>;
    #[inline]
    fn div(self, rhs: Jet2<f32>) -> Jet2<f32> {
        Jet2::constant(self) / rhs
    }
}

// --- Neg ---
impl<T: Passive> Neg for Jet2<T> {
    type Output = Jet2<T>;
    #[inline]
    fn neg(self) -> Self {
        Jet2 {
            value: -self.value,
            d1: -self.d1,
            d2: -self.d2,
        }
    }
}

// --- Compound assignment ---
impl<T: Passive> AddAssign for Jet2<T> {
    #[inline]
    fn add_assign(&mut self, rhs: Self) {
        *self = *self + rhs;
    }
}

impl<T: Passive> AddAssign<T> for Jet2<T> {
    #[inline]
    fn add_assign(&mut self, rhs: T) {
        self.value = self.value + rhs;
    }
}

impl<T: Passive> SubAssign for Jet2<T> {
    #[inline]
    fn sub_assign(&mut self, rhs: Self) {
        *self = *self - rhs;
    }
}

impl<T: Passive> SubAssign<T> for Jet2<T> {
    #[inline]
    fn sub_assign(&mut self, rhs: T) {
        self.value = self.value - rhs;
    }
}

impl<T: Passive> MulAssign for Jet2<T> {
    #[inline]
    fn mul_assign(&mut self, rhs: Self) {
        *self = *self * rhs;
    }
}

impl<T: Passive> MulAssign<T> for Jet2<T> {
    #[inline]
    fn mul_assign(&mut self, rhs: T) {
        *self = *self * rhs;
    }
}

impl<T: Passive> DivAssign for Jet2<T> {
    #[inline]
    fn div_assign(&mut self, rhs: Self) {
        *self = *self / rhs;
    }
}

impl<T: Passive> DivAssign<T> for Jet2<T> {
    #[inline]
    fn div_assign(&mut self, rhs: T) {
        *self = *self / rhs;
    }
}

// --- Display / Debug / Default ---
impl<T: Passive> fmt::Display for Jet2<T> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{}", self.value)
    }
}

impl<T: Passive> fmt::Debug for Jet2<T> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "Jet2(v={}, d1={}, d2={})",
            self.value, self.d1, self.d2
        )
    }
}

impl<T: Passive> Default for Jet2<T> {
    fn default() -> Self {
        Jet2::constant(T::zero())
    }
}

impl<T: Passive> From<T> for Jet2<T> {
    fn from(value: T) -> Self {
        Jet2::constant(value)
    }
}

impl From<i32> for Jet2<f64> {
    fn from(value: i32) -> Self {
        Jet2::constant(value as f64)
    }
}

impl From<i32> for Jet2<f32> {
    fn from(value: i32) -> Self {
        Jet2::constant(value as f32)
    }
}

impl<T: Passive> PartialEq for Jet2<T> {
    fn eq(&self, other: &Self) -> bool {
        self.value == other.value
    }
}

impl<T: Passive> PartialOrd for Jet2<T> {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        self.value.partial_cmp(&other.value)
    }
}

// ============================================================================
// NamedJet2<T> — named wrapper over seeded Jet2<T>
// ============================================================================

/// Labeled wrapper around the positional [`Jet2<T>`] type.
///
/// # Example
///
/// The constructor API for `NamedJet2<T>` is owned by
/// `NamedForwardTape<T>` via the `declare_jet2` / `into_scope` /
/// `scope.jet2(handle)` pattern. See the module docs in
/// `src/forward_tape.rs` for the full rationale.
///
/// ```
/// use xad_rs::{NamedForwardScope, NamedForwardTape};
///
/// let mut ft = NamedForwardTape::<f64>::new();
/// let x_h = ft.declare_jet2("x", 3.0);
/// let scope: NamedForwardScope<f64> = ft.into_scope();
///
/// let x = scope.jet2(x_h);
/// let f = x * x; // f(x) = x^2, f'(x) = 2x = 6, f''(x) = 2
///
/// assert_eq!(f.value(), 9.0);
/// assert_eq!(f.first_derivative("x"), 6.0);
/// assert_eq!(f.second_derivative("x"), 2.0);
/// ```
#[derive(Clone)]
pub struct NamedJet2<T: Passive> {
    pub(crate) inner: Jet2<T>,
    pub(crate) seeded: Option<usize>,
    // NOTE: field name is `gen_id` — `gen` alone is a reserved keyword in
    // Rust 2024 edition.
    #[cfg(debug_assertions)]
    pub(crate) gen_id: u64,
}

impl<T: Passive> NamedJet2<T> {
    /// Internal constructor used by `NamedForwardTape` input / freeze
    /// paths. Reads the TLS active generation (debug builds only) to stamp
    /// the `gen_id` field. Not part of the public API.
    #[inline]
    pub(crate) fn __from_parts(inner: Jet2<T>, seeded: Option<usize>) -> Self {
        Self {
            inner,
            seeded,
            #[cfg(debug_assertions)]
            gen_id: crate::forward_tape::current_gen(),
        }
    }

    /// Value part.
    #[inline]
    pub fn value(&self) -> T {
        self.inner.value()
    }

    /// First derivative with respect to `name`.
    ///
    /// Returns the seeded first derivative if `name` matches the seeded
    /// variable, else `T::zero()`. Reads the active registry from the
    /// `NamedForwardTape` thread-local slot. Panics if `name` is not in
    /// the registry, or if called outside a frozen `NamedForwardTape`
    /// scope.
    pub fn first_derivative(&self, name: &str) -> T {
        let idx = crate::forward_tape::with_active_registry(|r| {
            let r = r.expect(
                "NamedJet2::first_derivative called outside a frozen NamedForwardTape scope",
            );
            r.index_of(name).unwrap_or_else(|| {
                panic!(
                    "NamedJet2::first_derivative: name {:?} not present in registry",
                    name
                )
            })
        });
        if self.seeded == Some(idx) {
            self.inner.first_derivative()
        } else {
            T::zero()
        }
    }

    /// Second derivative (diagonal Hessian entry) with respect to `name`.
    ///
    /// Returns the seeded second derivative if `name` matches the seeded
    /// variable, else `T::zero()`. Panics if `name` is not in the registry,
    /// or if called outside a frozen `NamedForwardTape` scope.
    pub fn second_derivative(&self, name: &str) -> T {
        let idx = crate::forward_tape::with_active_registry(|r| {
            let r = r.expect(
                "NamedJet2::second_derivative called outside a frozen NamedForwardTape scope",
            );
            r.index_of(name).unwrap_or_else(|| {
                panic!(
                    "NamedJet2::second_derivative: name {:?} not present in registry",
                    name
                )
            })
        });
        if self.seeded == Some(idx) {
            self.inner.second_derivative()
        } else {
            T::zero()
        }
    }

    /// Escape hatch: direct access to the inner positional `Jet2<T>`.
    #[inline]
    pub fn inner(&self) -> &Jet2<T> {
        &self.inner
    }

    // ============ Elementary math delegations ============
    // Each method forwards to the inherent `Jet2<T>` elementary (which
    // takes `self` by value — `Jet2<T>` is `Copy`), preserves the
    // `seeded` field because unary elementaries cannot change which
    // direction is active, and stamps the parent's generation explicitly.

    /// Natural exponential, preserving the seed and the parent's generation.
    #[inline]
    pub fn exp(&self) -> Self {
        Self {
            inner: self.inner.exp(),
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }

    /// Natural logarithm, preserving the seed and the parent's generation.
    #[inline]
    pub fn ln(&self) -> Self {
        Self {
            inner: self.inner.ln(),
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }

    /// Square root, preserving the seed and the parent's generation.
    #[inline]
    pub fn sqrt(&self) -> Self {
        Self {
            inner: self.inner.sqrt(),
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }

    /// Sine, preserving the seed and the parent's generation.
    #[inline]
    pub fn sin(&self) -> Self {
        Self {
            inner: self.inner.sin(),
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }

    /// Cosine, preserving the seed and the parent's generation.
    #[inline]
    pub fn cos(&self) -> Self {
        Self {
            inner: self.inner.cos(),
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }

    /// Standard normal CDF, preserving the seed and the parent's generation.
    #[inline]
    pub fn norm_cdf(&self) -> Self {
        Self {
            inner: self.inner.norm_cdf(),
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }

    /// Inverse standard normal CDF, preserving the seed and the parent's generation.
    #[inline]
    pub fn inv_norm_cdf(&self) -> Self {
        Self {
            inner: self.inner.inv_norm_cdf(),
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }
}

impl<T: Passive> fmt::Debug for NamedJet2<T> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.debug_struct("NamedJet2")
            .field("value", &self.inner.value())
            .field("first", &self.inner.first_derivative())
            .field("second", &self.inner.second_derivative())
            .field("seeded", &self.seeded)
            .finish()
    }
}

impl<T: Passive> fmt::Display for NamedJet2<T> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "NamedJet2({})", self.inner.value())
    }
}

/// Merge two seeded directions. Same seed or one constant is fine;
/// two different seeds is a `Jet2` semantic violation — debug-panic,
/// release-pick-LHS.
#[inline]
pub(crate) fn merge_seeded(a: Option<usize>, b: Option<usize>) -> Option<usize> {
    match (a, b) {
        (None, None) => None,
        (Some(x), None) | (None, Some(x)) => Some(x),
        (Some(x), Some(y)) if x == y => Some(x),
        (Some(_), Some(_)) => {
            #[cfg(debug_assertions)]
            panic!(
                "NamedJet2: operation between two differently-seeded variables; \
                 seeded Jet2 supports only one active direction"
            );
            #[cfg(not(debug_assertions))]
            a
        }
    }
}

// ============================ Operator impls for NamedJet2 ============================
// Hand-written local macro (no shared stamping macro) because of the
// `seeded` field merge logic.
// `Jet2<T>` is `Copy`, so inner-value ops are by-value throughout.
// Each wrapper-vs-wrapper impl performs a debug-only `check_gen` between
// the two operands' generations, then merges `seeded` via `merge_seeded`,
// then constructs the result preserving the LHS's generation stamp.
// Pattern: 6 variants per binary op × 4 ops = 24 impls, plus 2 Neg impls.

macro_rules! __named_d2_binop {
    ($trait_:ident, $method:ident, $op:tt) => {
        // owned + owned
        impl<T: Passive> ::core::ops::$trait_<NamedJet2<T>> for NamedJet2<T> {
            type Output = NamedJet2<T>;
            #[inline]
            fn $method(self, rhs: NamedJet2<T>) -> NamedJet2<T> {
                #[cfg(debug_assertions)]
                crate::forward_tape::check_gen(self.gen_id, rhs.gen_id);
                NamedJet2 {
                    inner: self.inner $op rhs.inner,
                    seeded: merge_seeded(self.seeded, rhs.seeded),
                    #[cfg(debug_assertions)]
                    gen_id: self.gen_id,
                }
            }
        }
        // ref + ref
        impl<T: Passive> ::core::ops::$trait_<&NamedJet2<T>> for &NamedJet2<T> {
            type Output = NamedJet2<T>;
            #[inline]
            fn $method(self, rhs: &NamedJet2<T>) -> NamedJet2<T> {
                #[cfg(debug_assertions)]
                crate::forward_tape::check_gen(self.gen_id, rhs.gen_id);
                NamedJet2 {
                    inner: self.inner $op rhs.inner,
                    seeded: merge_seeded(self.seeded, rhs.seeded),
                    #[cfg(debug_assertions)]
                    gen_id: self.gen_id,
                }
            }
        }
        // owned + ref
        impl<T: Passive> ::core::ops::$trait_<&NamedJet2<T>> for NamedJet2<T> {
            type Output = NamedJet2<T>;
            #[inline]
            fn $method(self, rhs: &NamedJet2<T>) -> NamedJet2<T> {
                #[cfg(debug_assertions)]
                crate::forward_tape::check_gen(self.gen_id, rhs.gen_id);
                NamedJet2 {
                    inner: self.inner $op rhs.inner,
                    seeded: merge_seeded(self.seeded, rhs.seeded),
                    #[cfg(debug_assertions)]
                    gen_id: self.gen_id,
                }
            }
        }
        // ref + owned
        impl<T: Passive> ::core::ops::$trait_<NamedJet2<T>> for &NamedJet2<T> {
            type Output = NamedJet2<T>;
            #[inline]
            fn $method(self, rhs: NamedJet2<T>) -> NamedJet2<T> {
                #[cfg(debug_assertions)]
                crate::forward_tape::check_gen(self.gen_id, rhs.gen_id);
                NamedJet2 {
                    inner: self.inner $op rhs.inner,
                    seeded: merge_seeded(self.seeded, rhs.seeded),
                    #[cfg(debug_assertions)]
                    gen_id: self.gen_id,
                }
            }
        }
        // X op T (scalar, no cross-registry check)
        impl<T: Passive> ::core::ops::$trait_<T> for NamedJet2<T> {
            type Output = NamedJet2<T>;
            #[inline]
            fn $method(self, rhs: T) -> NamedJet2<T> {
                NamedJet2 {
                    inner: self.inner $op rhs,
                    seeded: self.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: self.gen_id,
                }
            }
        }
        impl<T: Passive> ::core::ops::$trait_<T> for &NamedJet2<T> {
            type Output = NamedJet2<T>;
            #[inline]
            fn $method(self, rhs: T) -> NamedJet2<T> {
                NamedJet2 {
                    inner: self.inner $op rhs,
                    seeded: self.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: self.gen_id,
                }
            }
        }
    };
}

__named_d2_binop!(Add, add, +);
__named_d2_binop!(Sub, sub, -);
__named_d2_binop!(Mul, mul, *);
__named_d2_binop!(Div, div, /);

impl<T: Passive> ::core::ops::Neg for NamedJet2<T> {
    type Output = NamedJet2<T>;
    #[inline]
    fn neg(self) -> NamedJet2<T> {
        NamedJet2 {
            inner: -self.inner,
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }
}
impl<T: Passive> ::core::ops::Neg for &NamedJet2<T> {
    type Output = NamedJet2<T>;
    #[inline]
    fn neg(self) -> NamedJet2<T> {
        NamedJet2 {
            inner: -self.inner,
            seeded: self.seeded,
            #[cfg(debug_assertions)]
            gen_id: self.gen_id,
        }
    }
}

// ============ Scalar-on-LHS impls (f64 and f32) ============
// The inner `Jet2<T>` only provides owned `f64 op Jet2<f64>` / `f32 op
// Jet2<f32>` (no ref variants). The named ref variants dereference the
// `Copy` inner to call the owned op.

macro_rules! __named_d2_scalar_lhs {
    ($scalar:ty) => {
        impl ::core::ops::Add<NamedJet2<$scalar>> for $scalar {
            type Output = NamedJet2<$scalar>;
            #[inline]
            fn add(self, rhs: NamedJet2<$scalar>) -> NamedJet2<$scalar> {
                NamedJet2 {
                    inner: self + rhs.inner,
                    seeded: rhs.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: rhs.gen_id,
                }
            }
        }
        impl ::core::ops::Add<&NamedJet2<$scalar>> for $scalar {
            type Output = NamedJet2<$scalar>;
            #[inline]
            fn add(self, rhs: &NamedJet2<$scalar>) -> NamedJet2<$scalar> {
                NamedJet2 {
                    inner: self + rhs.inner,
                    seeded: rhs.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: rhs.gen_id,
                }
            }
        }
        impl ::core::ops::Sub<NamedJet2<$scalar>> for $scalar {
            type Output = NamedJet2<$scalar>;
            #[inline]
            fn sub(self, rhs: NamedJet2<$scalar>) -> NamedJet2<$scalar> {
                NamedJet2 {
                    inner: self - rhs.inner,
                    seeded: rhs.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: rhs.gen_id,
                }
            }
        }
        impl ::core::ops::Sub<&NamedJet2<$scalar>> for $scalar {
            type Output = NamedJet2<$scalar>;
            #[inline]
            fn sub(self, rhs: &NamedJet2<$scalar>) -> NamedJet2<$scalar> {
                NamedJet2 {
                    inner: self - rhs.inner,
                    seeded: rhs.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: rhs.gen_id,
                }
            }
        }
        impl ::core::ops::Mul<NamedJet2<$scalar>> for $scalar {
            type Output = NamedJet2<$scalar>;
            #[inline]
            fn mul(self, rhs: NamedJet2<$scalar>) -> NamedJet2<$scalar> {
                NamedJet2 {
                    inner: self * rhs.inner,
                    seeded: rhs.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: rhs.gen_id,
                }
            }
        }
        impl ::core::ops::Mul<&NamedJet2<$scalar>> for $scalar {
            type Output = NamedJet2<$scalar>;
            #[inline]
            fn mul(self, rhs: &NamedJet2<$scalar>) -> NamedJet2<$scalar> {
                NamedJet2 {
                    inner: self * rhs.inner,
                    seeded: rhs.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: rhs.gen_id,
                }
            }
        }
        impl ::core::ops::Div<NamedJet2<$scalar>> for $scalar {
            type Output = NamedJet2<$scalar>;
            #[inline]
            fn div(self, rhs: NamedJet2<$scalar>) -> NamedJet2<$scalar> {
                NamedJet2 {
                    inner: self / rhs.inner,
                    seeded: rhs.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: rhs.gen_id,
                }
            }
        }
        impl ::core::ops::Div<&NamedJet2<$scalar>> for $scalar {
            type Output = NamedJet2<$scalar>;
            #[inline]
            fn div(self, rhs: &NamedJet2<$scalar>) -> NamedJet2<$scalar> {
                NamedJet2 {
                    inner: self / rhs.inner,
                    seeded: rhs.seeded,
                    #[cfg(debug_assertions)]
                    gen_id: rhs.gen_id,
                }
            }
        }
    };
}

__named_d2_scalar_lhs!(f64);
__named_d2_scalar_lhs!(f32);

// ============================================================================
// impl Real for Jet2<f64>
// ============================================================================
//
// `Jet2<T>` is `Copy` (when `T: Copy`) and exposes its transcendentals as
// value-taking inherent methods. The `Real` trait's `&self`-taking methods
// dereference through `Copy` to reach them. `Real::powf(&self, exponent: Self)`
// requires a two-`Jet2` power, which is implemented as
// `exp(exponent * ln(self))` — both factors propagate first and second
// derivatives through the existing primitives. `Real::powi` delegates to the
// inherent `Jet2::powf(self, T)` with a lossless `i32` → `f64` cast.

impl crate::real::Real for Jet2<f64> {
    type Passive = f64;

    #[inline]
    fn value(&self) -> f64 {
        self.value
    }
    #[inline]
    fn ln(&self) -> Self {
        Jet2::ln(*self)
    }
    #[inline]
    fn exp(&self) -> Self {
        Jet2::exp(*self)
    }
    #[inline]
    fn sqrt(&self) -> Self {
        Jet2::sqrt(*self)
    }
    #[inline]
    fn powf(&self, exponent: Self) -> Self {
        // u^v = exp(v * ln(u)), composed via existing Jet2 primitives.
        Jet2::exp(exponent * Jet2::ln(*self))
    }
    #[inline]
    fn powi(&self, exponent: i32) -> Self {
        Jet2::powf(*self, exponent as f64)
    }
    #[inline]
    fn sin(&self) -> Self {
        Jet2::sin(*self)
    }
    #[inline]
    fn cos(&self) -> Self {
        Jet2::cos(*self)
    }
}