echidna 0.9.0

//! `std::ops` implementations for `Laurent<F, K>`.

use std::ops::{
    Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign,
};

use crate::float::Float;
use crate::laurent::Laurent;
use crate::taylor_ops;

// ══════════════════════════════════════════════
//  Laurent<F, K> ↔ Laurent<F, K>
// ══════════════════════════════════════════════

/// # Truncation behavior
///
/// When the pole-order gap between the two operands exceeds `K - 1`, the
/// lower-order coefficients of the smaller-magnitude operand are silently
/// truncated to fit the fixed-size `[F; K]` array. This is inherent to the
/// fixed-width Laurent representation and cannot be avoided without dynamic
/// allocation.
impl<F: Float, const K: usize> Add for Laurent<F, K> {
    type Output = Self;
    #[inline]
    fn add(self, rhs: Self) -> Self {
        // `Laurent::zero()` is the additive identity, but when the other
        // operand has a large-magnitude pole_order the shift-alignment
        // assertion below fires before the algebraic identity has a chance to
        // apply. Short-circuit on either operand being all-zero so the `Zero`
        // trait contract (and generic `.sum()` code) works for every
        // pole_order that fits in `i32`.
        if rhs.is_all_zero_pub() {
            return self;
        }
        if self.is_all_zero_pub() {
            return rhs;
        }
        let p1 = self.pole_order();
        let p2 = rhs.pole_order();
        let p_out = p1.min(p2);
        // Align both to p_out by shifting
        let shift1 = (p1 - p_out) as usize;
        let shift2 = (p2 - p_out) as usize;
        assert!(
            shift1 < K && shift2 < K,
            "Laurent Add: pole-order gap ({}) exceeds K-1 ({}), coefficients would be silently truncated",
            shift1.max(shift2),
            K - 1,
        );

        let a: [F; K] = std::array::from_fn(|i| {
            if i >= shift1 && i - shift1 < K {
                self.coeff(p_out + i as i32)
            } else {
                F::zero()
            }
        });
        let b: [F; K] = std::array::from_fn(|i| {
            if i >= shift2 && i - shift2 < K {
                rhs.coeff(p_out + i as i32)
            } else {
                F::zero()
            }
        });

        let mut c = [F::zero(); K];
        taylor_ops::taylor_add(&a, &b, &mut c);
        Laurent::new(c, p_out)
    }
}

impl<F: Float, const K: usize> Sub for Laurent<F, K> {
    type Output = Self;
    #[inline]
    fn sub(self, rhs: Self) -> Self {
        // Short-circuit on zero operand (see `Add` impl for rationale).
        if rhs.is_all_zero_pub() {
            return self;
        }
        if self.is_all_zero_pub() {
            return -rhs;
        }
        let p1 = self.pole_order();
        let p2 = rhs.pole_order();
        let p_out = p1.min(p2);
        let shift1 = (p1 - p_out) as usize;
        let shift2 = (p2 - p_out) as usize;
        assert!(
            shift1 < K && shift2 < K,
            "Laurent Sub: pole-order gap ({}) exceeds K-1 ({}), coefficients would be silently truncated",
            shift1.max(shift2),
            K - 1,
        );

        let a: [F; K] = std::array::from_fn(|i| {
            if i >= shift1 && i - shift1 < K {
                self.coeff(p_out + i as i32)
            } else {
                F::zero()
            }
        });
        let b: [F; K] = std::array::from_fn(|i| {
            if i >= shift2 && i - shift2 < K {
                rhs.coeff(p_out + i as i32)
            } else {
                F::zero()
            }
        });

        let mut c = [F::zero(); K];
        taylor_ops::taylor_sub(&a, &b, &mut c);
        Laurent::new(c, p_out)
    }
}

// Truncated Laurent series multiplication uses the Cauchy product, which accumulates
// terms via addition — clippy flags the + inside a Mul impl, but this is correct for
// power-series coefficient propagation.
#[allow(clippy::suspicious_arithmetic_impl)]
impl<F: Float, const K: usize> Mul for Laurent<F, K> {
    type Output = Self;
    #[inline]
    fn mul(self, rhs: Self) -> Self {
        let mut c = [F::zero(); K];
        taylor_ops::taylor_mul(&self.leading_coeffs(), &rhs.leading_coeffs(), &mut c);
        // Pole orders can legitimately be large on near-singular inputs; use
        // checked arithmetic so `i32::MAX + 1` folds to a degenerate NaN
        // Laurent instead of silently wrapping into a negative order.
        let p = match self.pole_order().checked_add(rhs.pole_order()) {
            Some(p) => p,
            None => return Laurent::nan_pub(),
        };
        Laurent::new(c, p)
    }
}

// Truncated Laurent series division computes coefficients via recurrence that uses
// multiplication internally — clippy flags the * inside a Div impl, but this is correct
// for power-series coefficient propagation.
#[allow(clippy::suspicious_arithmetic_impl)]
impl<F: Float, const K: usize> Div for Laurent<F, K> {
    type Output = Self;
    #[inline]
    fn div(self, rhs: Self) -> Self {
        if rhs.is_all_zero_pub() {
            return Laurent::nan_pub();
        }
        let mut c = [F::zero(); K];
        taylor_ops::taylor_div(&self.leading_coeffs(), &rhs.leading_coeffs(), &mut c);
        let p = match self.pole_order().checked_sub(rhs.pole_order()) {
            Some(p) => p,
            None => return Laurent::nan_pub(),
        };
        Laurent::new(c, p)
    }
}

impl<F: Float, const K: usize> Neg for Laurent<F, K> {
    type Output = Self;
    #[inline]
    fn neg(self) -> Self {
        let mut c = [F::zero(); K];
        taylor_ops::taylor_neg(&self.leading_coeffs(), &mut c);
        Laurent::new(c, self.pole_order())
    }
}

impl<F: Float, const K: usize> Rem for Laurent<F, K> {
    type Output = Self;
    #[inline]
    fn rem(self, rhs: Self) -> Self {
        // a % b = a - (a / b).trunc() * b
        let quotient = (self / rhs).trunc();
        self - quotient * rhs
    }
}

impl<F: Float, const K: usize> AddAssign for Laurent<F, K> {
    #[inline]
    fn add_assign(&mut self, rhs: Self) {
        *self = *self + rhs;
    }
}

impl<F: Float, const K: usize> SubAssign for Laurent<F, K> {
    #[inline]
    fn sub_assign(&mut self, rhs: Self) {
        *self = *self - rhs;
    }
}

impl<F: Float, const K: usize> MulAssign for Laurent<F, K> {
    #[inline]
    fn mul_assign(&mut self, rhs: Self) {
        *self = *self * rhs;
    }
}

impl<F: Float, const K: usize> DivAssign for Laurent<F, K> {
    #[inline]
    fn div_assign(&mut self, rhs: Self) {
        *self = *self / rhs;
    }
}

impl<F: Float, const K: usize> RemAssign for Laurent<F, K> {
    #[inline]
    fn rem_assign(&mut self, rhs: Self) {
        *self = *self % rhs;
    }
}

// Mixed ops: Laurent<F, K> with primitive floats.
macro_rules! impl_laurent_scalar_ops {
    ($f:ty) => {
        impl<const K: usize> Add<$f> for Laurent<$f, K> {
            type Output = Laurent<$f, K>;
            #[inline]
            fn add(self, rhs: $f) -> Laurent<$f, K> {
                self + Laurent::constant(rhs)
            }
        }

        impl<const K: usize> Add<Laurent<$f, K>> for $f {
            type Output = Laurent<$f, K>;
            #[inline]
            fn add(self, rhs: Laurent<$f, K>) -> Laurent<$f, K> {
                Laurent::constant(self) + rhs
            }
        }

        impl<const K: usize> Sub<$f> for Laurent<$f, K> {
            type Output = Laurent<$f, K>;
            #[inline]
            fn sub(self, rhs: $f) -> Laurent<$f, K> {
                self - Laurent::constant(rhs)
            }
        }

        impl<const K: usize> Sub<Laurent<$f, K>> for $f {
            type Output = Laurent<$f, K>;
            #[inline]
            fn sub(self, rhs: Laurent<$f, K>) -> Laurent<$f, K> {
                Laurent::constant(self) - rhs
            }
        }

        impl<const K: usize> Mul<$f> for Laurent<$f, K> {
            type Output = Laurent<$f, K>;
            #[inline]
            fn mul(self, rhs: $f) -> Laurent<$f, K> {
                self * Laurent::constant(rhs)
            }
        }

        impl<const K: usize> Mul<Laurent<$f, K>> for $f {
            type Output = Laurent<$f, K>;
            #[inline]
            fn mul(self, rhs: Laurent<$f, K>) -> Laurent<$f, K> {
                Laurent::constant(self) * rhs
            }
        }

        // Scalar division delegates to full Laurent Div to reuse pole-order arithmetic.
        // Clippy flags the / inside a Div impl (self-referential dispatch), but this is
        // intentional — the scalar is promoted to a Laurent series first.
        #[allow(clippy::suspicious_arithmetic_impl)]
        impl<const K: usize> Div<$f> for Laurent<$f, K> {
            type Output = Laurent<$f, K>;
            #[inline]
            fn div(self, rhs: $f) -> Laurent<$f, K> {
                self / Laurent::constant(rhs)
            }
        }

        impl<const K: usize> Div<Laurent<$f, K>> for $f {
            type Output = Laurent<$f, K>;
            #[inline]
            fn div(self, rhs: Laurent<$f, K>) -> Laurent<$f, K> {
                Laurent::constant(self) / rhs
            }
        }

        impl<const K: usize> Rem<$f> for Laurent<$f, K> {
            type Output = Laurent<$f, K>;
            #[inline]
            fn rem(self, rhs: $f) -> Laurent<$f, K> {
                self % Laurent::constant(rhs)
            }
        }

        impl<const K: usize> Rem<Laurent<$f, K>> for $f {
            type Output = Laurent<$f, K>;
            #[inline]
            fn rem(self, rhs: Laurent<$f, K>) -> Laurent<$f, K> {
                Laurent::constant(self) % rhs
            }
        }
    };
}

impl_laurent_scalar_ops!(f32);
impl_laurent_scalar_ops!(f64);

impl<F: Float, const K: usize> PartialEq for Laurent<F, K> {
    #[inline]
    fn eq(&self, other: &Self) -> bool {
        self.value() == other.value()
    }
}

impl<F: Float, const K: usize> PartialOrd for Laurent<F, K> {
    #[inline]
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        self.value().partial_cmp(&other.value())
    }
}