equanetwork-math 0.0.4

//! Inventory skew: dead zone + polynomial intensity, path-averaged over a fill.

#[cfg(feature = "wasm")]
use equanetwork_macros::wasm_expose;

use super::consts::{BPS_DENOMINATOR, MAX_VAULT_SKEW_PER_1M, PER_1M_DENOMINATOR};
use super::error::{
    CoreError, AMOUNT_EXCEEDS_MAX_I32, ARITHMETIC_OVERFLOW, DIVISION_BY_ZERO, INVALID_SKEW,
};

/// Polynomial degree for inventory skew outside the dead zone.
#[derive(Default, Debug, Clone, Copy, Eq, PartialEq)]
#[cfg_attr(feature = "wasm", wasm_expose)]
pub enum SkewExponent {
    #[default]
    Linear,
    Quadratic,
    Cubic,
}

impl SkewExponent {
    pub fn value(&self) -> u32 {
        match self {
            Self::Linear => 1,
            Self::Quadratic => 2,
            Self::Cubic => 3,
        }
    }
}

/// Band mid and per-side dead-zone edges (token amounts).
#[derive(Debug, Clone, Copy, Eq, PartialEq)]
struct SkewBand {
    mid: u64,
    lower_start: u64,
    upper_start: u64,
    limit_lower: u64,
    limit_upper: u64,
}

impl SkewBand {
    fn new(
        limit_lower: u64,
        limit_upper: u64,
        lower_skew_offset_bps: u16,
        upper_skew_offset_bps: u16,
    ) -> Result<Option<Self>, CoreError> {
        if lower_skew_offset_bps as u64 > BPS_DENOMINATOR
            || upper_skew_offset_bps as u64 > BPS_DENOMINATOR
        {
            return Err(INVALID_SKEW);
        }
        let band = limit_upper.saturating_sub(limit_lower);
        if band == 0 {
            return Ok(None);
        }
        let half = band / 2;
        if half == 0 {
            return Ok(None);
        }
        let mid = limit_lower.saturating_add(half);
        let lower_span = mid.saturating_sub(limit_lower);
        let upper_span = limit_upper.saturating_sub(mid);
        let lower_offset = (lower_span as u128)
            .checked_mul(lower_skew_offset_bps as u128)
            .ok_or(ARITHMETIC_OVERFLOW)?
            .checked_div(BPS_DENOMINATOR as u128)
            .ok_or(DIVISION_BY_ZERO)?;
        let upper_offset = (upper_span as u128)
            .checked_mul(upper_skew_offset_bps as u128)
            .ok_or(ARITHMETIC_OVERFLOW)?
            .checked_div(BPS_DENOMINATOR as u128)
            .ok_or(DIVISION_BY_ZERO)?;
        let lower_start = mid.saturating_sub(lower_offset as u64);
        let upper_start = mid.saturating_add(upper_offset as u64);
        Ok(Some(Self {
            mid,
            lower_start,
            upper_start,
            limit_lower,
            limit_upper,
        }))
    }
}

/// Map inventory balance to signed deviation in ppm (−1e6..=+1e6).
/// Dead zone → 0; outer limits → ±1e6.
#[cfg_attr(feature = "wasm", wasm_expose)]
pub fn deviation_per_m(
    balance: u64,
    limit_lower: u64,
    limit_upper: u64,
    lower_skew_offset_bps: u16,
    upper_skew_offset_bps: u16,
) -> Result<i32, CoreError> {
    let Some(band) = SkewBand::new(
        limit_lower,
        limit_upper,
        lower_skew_offset_bps,
        upper_skew_offset_bps,
    )?
    else {
        return Ok(0);
    };
    deviation_per_m_in_band(balance, &band)
}

fn deviation_per_m_in_band(balance: u64, band: &SkewBand) -> Result<i32, CoreError> {
    if balance > band.upper_start {
        let span = band.limit_upper.saturating_sub(band.upper_start);
        if span == 0 {
            return Ok(0);
        }
        let dist = balance.saturating_sub(band.upper_start).min(span);
        let ppm = (dist as u128)
            .checked_mul(PER_1M_DENOMINATOR as u128)
            .ok_or(ARITHMETIC_OVERFLOW)?
            .checked_div(span as u128)
            .ok_or(DIVISION_BY_ZERO)?;
        let ppm = core::cmp::min(ppm, PER_1M_DENOMINATOR as u128);
        return i32::try_from(ppm).map_err(|_| AMOUNT_EXCEEDS_MAX_I32);
    }
    if balance < band.lower_start {
        let span = band.lower_start.saturating_sub(band.limit_lower);
        if span == 0 {
            return Ok(0);
        }
        let dist = band.lower_start.saturating_sub(balance).min(span);
        let ppm = (dist as u128)
            .checked_mul(PER_1M_DENOMINATOR as u128)
            .ok_or(ARITHMETIC_OVERFLOW)?
            .checked_div(span as u128)
            .ok_or(DIVISION_BY_ZERO)?;
        let ppm = core::cmp::min(ppm, PER_1M_DENOMINATOR as u128);
        let v = i32::try_from(ppm).map_err(|_| AMOUNT_EXCEEDS_MAX_I32)?;
        return Ok(-v);
    }
    Ok(0)
}

fn select_intensity(
    deviation_per_m: i32,
    increasing: bool,
    positive_skew_per_1m: u32,
    negative_skew_per_1m: u32,
) -> u32 {
    match (deviation_per_m >= 0, increasing) {
        (true, true) | (false, false) => positive_skew_per_1m,
        (false, true) | (true, false) => negative_skew_per_1m,
    }
}

/// Instantaneous skew (ppm) at a deviation. Sign follows deviation; round away from zero.
fn poly_skew_magnitude(
    abs_dev: u128,
    intensity: u32,
    exponent: SkewExponent,
) -> Result<u128, CoreError> {
    if intensity == 0 || abs_dev == 0 {
        return Ok(0);
    }
    let exp = exponent.value();
    let numerator = abs_dev
        .checked_pow(exp)
        .ok_or(ARITHMETIC_OVERFLOW)?
        .checked_mul(intensity as u128)
        .ok_or(ARITHMETIC_OVERFLOW)?;
    let denominator = (PER_1M_DENOMINATOR as u128)
        .checked_pow(exp)
        .ok_or(ARITHMETIC_OVERFLOW)?;
    let quotient = numerator.checked_div(denominator).ok_or(DIVISION_BY_ZERO)?;
    let remainder = numerator.checked_rem(denominator).ok_or(DIVISION_BY_ZERO)?;
    if remainder > 0 {
        quotient.checked_add(1).ok_or(ARITHMETIC_OVERFLOW)
    } else {
        Ok(quotient)
    }
}

/// Instantaneous vault skew contribution (ppm) at `balance` for a move in `increasing` direction.
#[allow(clippy::too_many_arguments)]
#[cfg_attr(feature = "wasm", wasm_expose)]
pub fn compute_skew_per_m(
    balance: u64,
    limit_lower: u64,
    limit_upper: u64,
    lower_skew_offset_bps: u16,
    upper_skew_offset_bps: u16,
    positive_skew_per_1m: u32,
    negative_skew_per_1m: u32,
    exponent: SkewExponent,
    increasing: bool,
) -> Result<i64, CoreError> {
    if positive_skew_per_1m > MAX_VAULT_SKEW_PER_1M || negative_skew_per_1m > MAX_VAULT_SKEW_PER_1M
    {
        return Err(INVALID_SKEW);
    }
    let deviation = deviation_per_m(
        balance,
        limit_lower,
        limit_upper,
        lower_skew_offset_bps,
        upper_skew_offset_bps,
    )?;
    if deviation == 0 {
        return Ok(0);
    }
    let intensity = select_intensity(
        deviation,
        increasing,
        positive_skew_per_1m,
        negative_skew_per_1m,
    );
    let abs_dev = deviation.unsigned_abs() as u128;
    let mag = poly_skew_magnitude(abs_dev, intensity, exponent)?;
    let mag = i64::try_from(mag).map_err(|_| AMOUNT_EXCEEDS_MAX_I32)?;
    let raw = (deviation.signum() as i64)
        .checked_mul(mag)
        .ok_or(ARITHMETIC_OVERFLOW)?;
    if increasing {
        Ok(raw)
    } else {
        raw.checked_neg().ok_or(ARITHMETIC_OVERFLOW)
    }
}

/// Signed intensity for the contribution rate on the upper / lower active region.
fn region_intensity_signed(
    upper: bool,
    increasing: bool,
    positive_skew_per_1m: u32,
    negative_skew_per_1m: u32,
) -> i64 {
    if upper {
        if increasing {
            positive_skew_per_1m as i64
        } else {
            -(negative_skew_per_1m as i64)
        }
    } else if increasing {
        -(negative_skew_per_1m as i64)
    } else {
        positive_skew_per_1m as i64
    }
}

/// ∫_a^b (t(x)/1e6)^n dx where t maps linearly from span endpoints to ppm.
/// `t_at` maps balance → deviation magnitude in ppm on this side.
fn integrate_poly_over_span(
    a: u64,
    b: u64,
    span_start: u64,
    span_end: u64,
    upper: bool,
    exponent: SkewExponent,
) -> Result<u128, CoreError> {
    if a >= b {
        return Ok(0);
    }
    let span = if upper {
        span_end.saturating_sub(span_start)
    } else {
        span_start.saturating_sub(span_end)
    };
    if span == 0 {
        return Ok(0);
    }
    let t = |bal: u64| -> Result<u128, CoreError> {
        let dist = if upper {
            bal.saturating_sub(span_start).min(span)
        } else {
            span_start.saturating_sub(bal).min(span)
        };
        let ppm = (dist as u128)
            .checked_mul(PER_1M_DENOMINATOR as u128)
            .ok_or(ARITHMETIC_OVERFLOW)?
            .checked_div(span as u128)
            .ok_or(DIVISION_BY_ZERO)?;
        Ok(core::cmp::min(ppm, PER_1M_DENOMINATOR as u128))
    };
    let t0 = t(a)?;
    let t1 = t(b)?;
    let n = exponent.value();
    let n1 = n.checked_add(1).ok_or(ARITHMETIC_OVERFLOW)?;
    // ∫ (t/PER)^n db = span / PER^{n+1} * (t1^{n+1} - t0^{n+1}) / (n+1)
    // For lower side, as b increases t decreases, so use |t0^{n+1} - t1^{n+1}|.
    let p0 = t0.checked_pow(n1).ok_or(ARITHMETIC_OVERFLOW)?;
    let p1 = t1.checked_pow(n1).ok_or(ARITHMETIC_OVERFLOW)?;
    let delta_p = p0.abs_diff(p1);
    let numer = (span as u128)
        .checked_mul(delta_p)
        .ok_or(ARITHMETIC_OVERFLOW)?;
    let denom = (PER_1M_DENOMINATOR as u128)
        .checked_pow(n1)
        .ok_or(ARITHMETIC_OVERFLOW)?
        .checked_mul(n1 as u128)
        .ok_or(ARITHMETIC_OVERFLOW)?;
    numer.checked_div(denom).ok_or(DIVISION_BY_ZERO)
}

fn integral_rate_segment(
    a: u64,
    b: u64,
    band: &SkewBand,
    increasing: bool,
    positive_skew_per_1m: u32,
    negative_skew_per_1m: u32,
    exponent: SkewExponent,
) -> Result<i128, CoreError> {
    if a >= b {
        return Ok(0);
    }

    let mut sum: i128 = 0;

    // Upper active region [upper_start, limit_upper] (t=0 at upper_start)
    let u0 = core::cmp::max(a, band.upper_start);
    let u1 = core::cmp::min(b, band.limit_upper);
    if u0 < u1 && band.limit_upper > band.upper_start {
        let poly =
            integrate_poly_over_span(u0, u1, band.upper_start, band.limit_upper, true, exponent)?;
        let intensity =
            region_intensity_signed(true, increasing, positive_skew_per_1m, negative_skew_per_1m);
        sum = sum
            .checked_add(
                (intensity as i128)
                    .checked_mul(poly as i128)
                    .ok_or(ARITHMETIC_OVERFLOW)?,
            )
            .ok_or(ARITHMETIC_OVERFLOW)?;
    }

    // Lower active region [limit_lower, lower_start] (t=0 at lower_start)
    let l0 = core::cmp::max(a, band.limit_lower);
    let l1 = core::cmp::min(b, band.lower_start);
    if l0 < l1 && band.lower_start > band.limit_lower {
        let poly =
            integrate_poly_over_span(l0, l1, band.lower_start, band.limit_lower, false, exponent)?;
        let intensity = region_intensity_signed(
            false,
            increasing,
            positive_skew_per_1m,
            negative_skew_per_1m,
        );
        sum = sum
            .checked_add(
                (intensity as i128)
                    .checked_mul(poly as i128)
                    .ok_or(ARITHMETIC_OVERFLOW)?,
            )
            .ok_or(ARITHMETIC_OVERFLOW)?;
    }

    Ok(sum)
}

/// Path-average skew (ppm) over `[balance_pre, balance_post]`.
#[allow(clippy::too_many_arguments)]
#[cfg_attr(feature = "wasm", wasm_expose)]
pub fn average_skew_per_m(
    balance_pre: u64,
    balance_post: u64,
    limit_lower: u64,
    limit_upper: u64,
    lower_skew_offset_bps: u16,
    upper_skew_offset_bps: u16,
    positive_skew_per_1m: u32,
    negative_skew_per_1m: u32,
    exponent: SkewExponent,
) -> Result<i64, CoreError> {
    if positive_skew_per_1m > MAX_VAULT_SKEW_PER_1M || negative_skew_per_1m > MAX_VAULT_SKEW_PER_1M
    {
        return Err(INVALID_SKEW);
    }
    if balance_pre == balance_post {
        return Ok(0);
    }
    let Some(band) = SkewBand::new(
        limit_lower,
        limit_upper,
        lower_skew_offset_bps,
        upper_skew_offset_bps,
    )?
    else {
        return Ok(0);
    };

    let increasing = balance_post > balance_pre;
    let lo = core::cmp::min(balance_pre, balance_post);
    let hi = core::cmp::max(balance_pre, balance_post);
    let delta = (hi as u128).saturating_sub(lo as u128);
    if delta == 0 {
        return Ok(0);
    }

    // Split at dead-zone edges for numerical stability.
    let mut points = [lo, band.lower_start, band.upper_start, hi];
    points.sort_unstable();
    let mut integral: i128 = 0;
    for w in points.windows(2) {
        let a = core::cmp::max(w[0], lo);
        let b = core::cmp::min(w[1], hi);
        if a >= b {
            continue;
        }
        integral = integral
            .checked_add(integral_rate_segment(
                a,
                b,
                &band,
                increasing,
                positive_skew_per_1m,
                negative_skew_per_1m,
                exponent,
            )?)
            .ok_or(ARITHMETIC_OVERFLOW)?;
    }

    // avg = integral / delta, round away from zero
    let abs_int = integral.unsigned_abs();
    let quot = abs_int.checked_div(delta).ok_or(DIVISION_BY_ZERO)?;
    let rem = abs_int.checked_rem(delta).ok_or(DIVISION_BY_ZERO)?;
    let abs_avg = if rem > 0 {
        quot.checked_add(1).ok_or(ARITHMETIC_OVERFLOW)?
    } else {
        quot
    };
    let abs_avg = i64::try_from(abs_avg).map_err(|_| AMOUNT_EXCEEDS_MAX_I32)?;
    if integral < 0 {
        abs_avg.checked_neg().ok_or(ARITHMETIC_OVERFLOW)
    } else {
        Ok(abs_avg)
    }
}

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

    #[test]
    fn dead_zone_zero_deviation() {
        // mid=500, offset 50% → lower_start=250, upper_start=750
        assert_eq!(deviation_per_m(500, 0, 1000, 5_000, 5_000).unwrap(), 0);
        assert_eq!(deviation_per_m(250, 0, 1000, 5_000, 5_000).unwrap(), 0);
        assert_eq!(deviation_per_m(750, 0, 1000, 5_000, 5_000).unwrap(), 0);
    }

    #[test]
    fn deviation_at_outer_limits() {
        assert_eq!(
            deviation_per_m(1000, 0, 1000, 5_000, 5_000).unwrap(),
            1_000_000
        );
        assert_eq!(
            deviation_per_m(0, 0, 1000, 5_000, 5_000).unwrap(),
            -1_000_000
        );
    }

    #[test]
    fn offset_full_disables_skew() {
        assert_eq!(deviation_per_m(0, 0, 1000, 10_000, 10_000).unwrap(), 0);
        assert_eq!(deviation_per_m(1000, 0, 1000, 10_000, 10_000).unwrap(), 0);
    }

    #[rstest]
    #[case(SkewExponent::Linear, 200_000, 1_000_000, 200_000)]
    #[case(SkewExponent::Quadratic, 200_000, 1_000_000, 40_000)]
    #[case(SkewExponent::Cubic, 200_000, 1_000_000, 8_000)]
    fn poly_point_values(
        #[case] exp: SkewExponent,
        #[case] abs_dev: u128,
        #[case] intensity: u32,
        #[case] expected: u128,
    ) {
        assert_eq!(
            poly_skew_magnitude(abs_dev, intensity, exp).unwrap(),
            expected
        );
    }

    #[test]
    fn worsening_above_mid_positive() {
        // offset 0: skew from mid. Move up from 800→900 on [0,1000], mid=500.
        let skew =
            average_skew_per_m(800, 900, 0, 1000, 0, 0, 10_000, 0, SkewExponent::Linear).unwrap();
        assert!(skew > 0, "worsening should be positive: {skew}");
    }

    #[test]
    fn rebalancing_above_mid_negative() {
        let skew =
            average_skew_per_m(900, 800, 0, 1000, 0, 0, 0, 10_000, SkewExponent::Linear).unwrap();
        assert!(skew < 0, "rebalancing should be negative: {skew}");
    }

    #[test]
    fn linear_average_matches_midpoint() {
        // Linear: average of endpoints equals path average (no dead zone).
        let pre = 600u64;
        let post = 800u64;
        let intensity = 50_000u32;
        let avg = average_skew_per_m(pre, post, 0, 1000, 0, 0, intensity, 0, SkewExponent::Linear)
            .unwrap();
        let s0 = compute_skew_per_m(pre, 0, 1000, 0, 0, intensity, 0, SkewExponent::Linear, true)
            .unwrap();
        let s1 = compute_skew_per_m(
            post,
            0,
            1000,
            0,
            0,
            intensity,
            0,
            SkewExponent::Linear,
            true,
        )
        .unwrap();
        // Continuous average of linear is midpoint; discrete endpoint mean is close.
        let mid = (s0 + s1) / 2;
        assert!((avg - mid).abs() <= 1, "avg={avg} mid={mid}");
    }

    #[test]
    fn quadratic_average_not_endpoint_mean() {
        let pre = 600u64;
        let post = 1000u64;
        let intensity = 50_000u32;
        let avg = average_skew_per_m(
            pre,
            post,
            0,
            1000,
            0,
            0,
            intensity,
            0,
            SkewExponent::Quadratic,
        )
        .unwrap();
        let s0 = compute_skew_per_m(
            pre,
            0,
            1000,
            0,
            0,
            intensity,
            0,
            SkewExponent::Quadratic,
            true,
        )
        .unwrap();
        let s1 = compute_skew_per_m(
            post,
            0,
            1000,
            0,
            0,
            intensity,
            0,
            SkewExponent::Quadratic,
            true,
        )
        .unwrap();
        let endpoint_mean = (s0 + s1) / 2;
        assert_ne!(avg, endpoint_mean);
        // For convex t^2, integral average is below the endpoint mean.
        assert!(
            avg < endpoint_mean,
            "avg={avg} endpoint_mean={endpoint_mean}"
        );
    }

    #[test]
    fn path_through_dead_zone() {
        // 50% dead zone: active only outside [250,750]
        let skew = average_skew_per_m(
            400,
            600,
            0,
            1000,
            5_000,
            5_000,
            10_000,
            10_000,
            SkewExponent::Linear,
        )
        .unwrap();
        assert_eq!(skew, 0);
    }
}