vsrg 0.3.0

use crate::math::easing::eval::easing_value;
use crate::math::{AdvancedEasing, Easable};

pub fn evaluate_integration_reciprocal<T, const N: usize>(
    ease: AdvancedEasing,
    from: T,
    to: T,
    t: T,
    samples: usize,
) -> Option<T>
where
    T: Easable<N>,
{
    let from = from.components();
    let to = to.components();
    let t = t.components();

    let mut out = [0.0f64; N];
    for i in 0..N {
        out[i] = integrate_reciprocal_scalar(ease, from[i], to[i], t[i], samples)?;
    }
    Some(T::compose(out))
}

fn integrate_reciprocal_scalar(
    ease: AdvancedEasing,
    from: f64,
    to: f64,
    t: f64,
    samples: usize,
) -> Option<f64> {
    if !t.is_finite() || !from.is_finite() || !to.is_finite() {
        return None;
    }

    if t < 0.0 {
        // Below the easing range the eased value is constant at `from`.
        div_checked(t, from)
    } else if t > 1.0 {
        // Integrate [0,1] then extend with the constant `to` tail.
        let inner = integrate_reciprocal_in_range(ease, from, to, 1.0, samples)?;
        let tail = div_checked(t - 1.0, to)?;
        Some(inner + tail)
    } else {
        integrate_reciprocal_in_range(ease, from, to, t, samples)
    }
}

/// Integrate `1 / (from*(1-e(s)) + to*e(s))` over `[0, t]` for `t` in `[0, 1]`.
///
/// The denominator is written as `from*(1-e(s)) + to*e(s)` to avoid computing
/// `to - from`, which can overflow to infinity for large same-signed endpoints.
fn integrate_reciprocal_in_range(
    ease: AdvancedEasing,
    from: f64,
    to: f64,
    t: f64,
    samples: usize,
) -> Option<f64> {
    if t <= 0.0 {
        return Some(0.0);
    }

    match ease {
        AdvancedEasing::Linear => {
            // e(s) = s, denom = from*(1-s) + to*s.
            if from == to {
                div_checked(t, from)
            } else if from * to <= 0.0 {
                // from and to have opposite signs (or one is zero) → the
                // denominator crosses zero in (0, 1) → pole.
                None
            } else {
                let denom_t = from * (1.0 - t) + to * t;
                Some(f64::ln(denom_t / from) / (to - from))
            }
        }

        AdvancedEasing::InConst => {
            // e(s) = 0 for s < 1, so the denominator is `from` on [0, 1).
            div_checked(t, from)
        }
        AdvancedEasing::OutConst => {
            // e(s) = 1 for s > 0, so the denominator is `to` on (0, 1].
            div_checked(t, to)
        }
        AdvancedEasing::InOutConst => {
            // e(s) = 0 for s < 0.5, = 1 for s >= 0.5.
            if t <= 0.5 {
                div_checked(t, from)
            } else {
                let first = div_checked(0.5, from)?;
                let second = div_checked(t - 0.5, to)?;
                Some(first + second)
            }
        }

        _ => {
            if has_pole(ease, from, to) {
                None
            } else {
                Some(simpson_reciprocal(ease, from, to, t, samples))
            }
        }
    }
}

/// Whether the easing's value `e(s)` exits `[0, 1]` for some `s ∈ [0, 1]`.
///
/// Only Back and Elastic variants overshoot. All other easings — including
/// Bounce — stay within `[0, 1]`.
#[allow(unused_variables)]
pub fn overshoots(ease: AdvancedEasing) -> bool {
    {
        matches!(
            ease,
            AdvancedEasing::InBack
                | AdvancedEasing::OutBack
                | AdvancedEasing::InOutBack
                | AdvancedEasing::InElastic
                | AdvancedEasing::OutElastic
                | AdvancedEasing::InOutElastic
        )
    }
}

/// The actual `(e_min, e_max)` range of `easing_value(ease, s)` for `s ∈ [0, 1]`,
/// for each easing variant. Pre-computed with 1,000,000 samples and widened
/// by 1e-6 to absorb floating-point noise at the boundaries.
///
/// All non-overshooting easings stay within `[0, 1]`. Overshooting easings
/// (Back, Elastic) have pre-computed bounds that are wider than `[0, 1]`.
pub fn overshoot_bounds(ease: AdvancedEasing) -> (f64, f64) {
    match ease {
        // e(0) = 0, dips to ≈ -0.10 around s ≈ 0.42, then rises to 1.
        AdvancedEasing::InBack => (-0.100005, 1.000001),
        // e(0) = 0, rises to ≈ 1.10 around s ≈ 0.58, then settles to 1.
        AdvancedEasing::OutBack => (-0.000001, 1.100005),
        // Mirrors InBack on the first half and OutBack on the second.
        AdvancedEasing::InOutBack => (-0.100152, 1.100152),
        // Oscillates with decaying amplitude; first trough reaches ≈ -1.0.
        AdvancedEasing::InElastic => (-1.000000, 1.000001),
        // Oscillates with decaying amplitude; first peak reaches ≈ 1.37.
        AdvancedEasing::OutElastic => (-0.000001, 1.373099),
        // Mirrors InElastic on the first half and OutElastic on the second.
        AdvancedEasing::InOutElastic => (-0.118349, 1.118349),
        _ => (0.0, 1.0),
    }
}

/// Check whether `1 / (from*(1-e(s)) + to*e(s))` has a pole for `s ∈ [0, 1]`.
///
/// The denominator `D(s) = from*(1-e(s)) + to*e(s)` is zero when
/// `e(s) = from / (from - to)` (the critical easing value). The function:
///
/// 1. If `from * to <= 0` (opposite signs or one is zero), the denominator
///    is `from` at `s=0` and `to` at `s=1`; by the intermediate value theorem
///    it crosses zero somewhere in `[0, 1]` → pole.
///
/// 2. If `from * to > 0` and the easing is bounded in `[0, 1]`, the
///    denominator is a convex combination of two same-signed values →
///    never zero → no pole.
///
/// 3. If `from * to > 0` and the easing overshoots, numerically compute the
///    easing's actual `[e_min, e_max]` range and check whether the critical
///    value `e_crit = from / (from - to)` falls within it.
fn has_pole(ease: AdvancedEasing, from: f64, to: f64) -> bool {
    if from * to <= 0.0 {
        return true;
    }
    if !overshoots(ease) {
        return false;
    }
    let (e_min, e_max) = overshoot_bounds(ease);
    let e_crit = from / (from - to);
    e_crit >= e_min && e_crit <= e_max
}

fn simpson_reciprocal(ease: AdvancedEasing, from: f64, to: f64, t: f64, samples: usize) -> f64 {
    let denom = |s: f64| from * (1.0 - easing_value(ease, s)) + to * easing_value(ease, s);
    let f = |s: f64| 1.0 / denom(s);

    let n = if samples < 2 {
        2
    } else if samples % 2 == 1 {
        samples + 1
    } else {
        samples
    };
    let h = t / n as f64;

    let mut sum = f(0.0) + f(t);
    for k in 1..n {
        let s = k as f64 * h;
        let w = if k % 2 == 0 { 2.0 } else { 4.0 };
        sum += w * f(s);
    }
    sum * h / 3.0
}

/// `a / b`, returning `None` if `b` is zero.
#[inline]
fn div_checked(a: f64, b: f64) -> Option<f64> {
    if b.abs() <= f64::EPSILON {
        None
    } else {
        Some(a / b)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::math::easing::eval::evaluate_easing;

    const FROM: f64 = 1.0;
    const TO: f64 = 100.0;
    const SAMPLES: usize = 2000;
    const H: f64 = 1e-4;
    const SAMPLE_TS: &[f64] = &[0.1, 0.25, 0.5, 0.75, 0.9];

    const CLOSED_FORM_EASINGS: &[AdvancedEasing] = &[
        AdvancedEasing::Linear,
        AdvancedEasing::InConst,
        AdvancedEasing::OutConst,
        AdvancedEasing::InOutConst,
    ];

    // Easings whose eased value stays within [0, 1] on [0, 1], so the
    // denominator `from*(1-e(s)) + to*e(s)` never vanishes when 0 < from < to.
    const WELL_BEHAVED_EASINGS: &[AdvancedEasing] = &[
        AdvancedEasing::Linear,
        AdvancedEasing::InSine,
        AdvancedEasing::OutSine,
        AdvancedEasing::InOutSine,
        AdvancedEasing::InQuad,
        AdvancedEasing::OutQuad,
        AdvancedEasing::InOutQuad,
        AdvancedEasing::InCubic,
        AdvancedEasing::OutCubic,
        AdvancedEasing::InOutCubic,
        AdvancedEasing::InQuart,
        AdvancedEasing::OutQuart,
        AdvancedEasing::InOutQuart,
        AdvancedEasing::InQuint,
        AdvancedEasing::OutQuint,
        AdvancedEasing::InOutQuint,
        AdvancedEasing::InExpo,
        AdvancedEasing::OutExpo,
        AdvancedEasing::InOutExpo,
        AdvancedEasing::InCirc,
        AdvancedEasing::OutCirc,
        AdvancedEasing::InOutCirc,
    ];

    fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() <= tol
    }

    #[test]
    fn t0_returns_zero() {
        for &ease in AdvancedEasing::ALL {
            let result = evaluate_integration_reciprocal::<f64, 1>(ease, FROM, TO, 0.0, SAMPLES);
            assert!(
                approx_eq(result.unwrap(), 0.0, 1e-6),
                "{ease:?}: expected reciprocal(0)=0, got {result:?}"
            );
        }
    }

    #[test]
    fn negative_t_returns_t_over_from() {
        for &ease in AdvancedEasing::ALL {
            let t = -0.3;
            let result = evaluate_integration_reciprocal::<f64, 1>(ease, FROM, TO, t, SAMPLES);
            assert!(
                approx_eq(result.unwrap(), t / FROM, 1e-6),
                "{ease:?}: expected t/from for t<0, got {result:?}"
            );
        }
    }

    #[test]
    fn t_above_1_extends_with_to() {
        for &ease in WELL_BEHAVED_EASINGS {
            let t = 1.4;
            let inner =
                evaluate_integration_reciprocal::<f64, 1>(ease, FROM, TO, 1.0, SAMPLES).unwrap();
            let result =
                evaluate_integration_reciprocal::<f64, 1>(ease, FROM, TO, t, SAMPLES).unwrap();
            assert!(
                approx_eq(result, inner + (t - 1.0) / TO, 1e-5),
                "{ease:?}: t>1 should be inner + (t-1)/to, got {result}"
            );
        }
    }

    #[test]
    fn derivative_matches_reciprocal_of_ease() {
        // d/dt R(t) = 1 / ease(from, to, t)
        let mut failures: Vec<(AdvancedEasing, f64, f64, f64)> = Vec::new();
        for &ease in WELL_BEHAVED_EASINGS {
            if matches!(
                ease,
                AdvancedEasing::InConst | AdvancedEasing::OutConst | AdvancedEasing::InOutConst
            ) {
                continue;
            }
            for &t in SAMPLE_TS {
                let vp = evaluate_integration_reciprocal::<f64, 1>(ease, FROM, TO, t + H, SAMPLES)
                    .unwrap();
                let vm = evaluate_integration_reciprocal::<f64, 1>(ease, FROM, TO, t - H, SAMPLES)
                    .unwrap();
                let numerical = (vp - vm) / (2.0 * H);
                let expected = 1.0 / evaluate_easing::<f64, 1>(ease, FROM, TO, t);
                if !approx_eq(numerical, expected, 1e-2) {
                    failures.push((ease, t, numerical, expected));
                }
            }
        }
        assert!(
            failures.is_empty(),
            "d/dt reciprocal did not match 1/ease: {}",
            failures
                .iter()
                .map(|(e, t, n, v)| format!("{e:?}@t={t}: d/dt={n}, expected={v}"))
                .collect::<Vec<_>>()
                .join("; ")
        );
    }

    #[test]
    fn linear_matches_closed_form() {
        // R(t) = (ln(from + delta*t) - ln(from)) / delta
        let delta = TO - FROM;
        for &t in SAMPLE_TS {
            let result = evaluate_integration_reciprocal::<f64, 1>(
                AdvancedEasing::Linear,
                FROM,
                TO,
                t,
                SAMPLES,
            )
            .unwrap();
            let expected = (f64::ln(FROM + delta * t) - f64::ln(FROM)) / delta;
            assert!(
                approx_eq(result, expected, 1e-4),
                "Linear: t={t}, expected {expected}, got {result}"
            );
        }
    }

    #[test]
    fn const_variants_match_closed_form() {
        for &t in SAMPLE_TS {
            let in_const = evaluate_integration_reciprocal::<f64, 1>(
                AdvancedEasing::InConst,
                FROM,
                TO,
                t,
                SAMPLES,
            )
            .unwrap();
            assert!(
                approx_eq(in_const, t / FROM, 1e-6),
                "InConst t={t}: {in_const}"
            );

            let out_const = evaluate_integration_reciprocal::<f64, 1>(
                AdvancedEasing::OutConst,
                FROM,
                TO,
                t,
                SAMPLES,
            )
            .unwrap();
            assert!(
                approx_eq(out_const, t / TO, 1e-6),
                "OutConst t={t}: {out_const}"
            );

            let inout_const = evaluate_integration_reciprocal::<f64, 1>(
                AdvancedEasing::InOutConst,
                FROM,
                TO,
                t,
                SAMPLES,
            )
            .unwrap();
            let expected = if t <= 0.5 {
                t / FROM
            } else {
                0.5 / FROM + (t - 0.5) / TO
            };
            assert!(
                approx_eq(inout_const, expected, 1e-6),
                "InOutConst t={t}: {inout_const}"
            );
        }
    }

    #[test]
    fn samples_count_does_not_break_closed_forms() {
        // Closed-form easings must not depend on `samples`.
        for &ease in CLOSED_FORM_EASINGS {
            for &t in SAMPLE_TS {
                let low = evaluate_integration_reciprocal::<f64, 1>(ease, FROM, TO, t, 2).unwrap();
                let high =
                    evaluate_integration_reciprocal::<f64, 1>(ease, FROM, TO, t, 10000).unwrap();
                assert!(
                    approx_eq(low, high, 1e-6),
                    "{ease:?}: closed form should be sample-invariant, {low} vs {high}"
                );
            }
        }
    }

    #[test]
    fn zero_from_returns_none() {
        // A zero lower bound makes the integrand diverge at s=0 for every
        // easing whose value is 0 at s=0. OutConst is the exception: its eased
        // value jumps to `to` immediately, so the integral stays finite.
        for &ease in AdvancedEasing::ALL {
            if matches!(ease, AdvancedEasing::OutConst) {
                continue;
            }
            let result = evaluate_integration_reciprocal::<f64, 1>(ease, 0.0, TO, 0.5, SAMPLES);
            assert!(
                result.is_none(),
                "{ease:?}: from=0 should diverge and return None, got {result:?}"
            );
        }
    }

    #[test]
    fn out_const_with_zero_from_stays_finite() {
        // OutConst's eased value is `to` on (0, 1]; from=0 doesn't cause a pole.
        let result = evaluate_integration_reciprocal::<f64, 1>(
            AdvancedEasing::OutConst,
            0.0,
            TO,
            0.5,
            SAMPLES,
        );
        assert!(
            approx_eq(result.unwrap(), 0.5 / TO, 1e-6),
            "OutConst with from=0 should be t/to, got {result:?}"
        );
    }

    #[test]
    fn zero_to_returns_none_for_t_above_1() {
        // The constant tail (t-1)/to diverges when to == 0.
        for &ease in AdvancedEasing::ALL {
            let result = evaluate_integration_reciprocal::<f64, 1>(ease, FROM, 0.0, 1.5, SAMPLES);
            assert!(
                result.is_none(),
                "{ease:?}: to=0 with t>1 should diverge and return None, got {result:?}"
            );
        }
    }

    #[test]
    fn same_sign_no_pole_for_bounded_easings() {
        // from*to > 0 and bounded easing → no pole, even with extreme ratios.
        let from: f64 = 1e-10;
        let to: f64 = 1e10;
        for &ease in WELL_BEHAVED_EASINGS {
            let result = evaluate_integration_reciprocal::<f64, 1>(ease, from, to, 0.5, SAMPLES);
            assert!(
                result.is_some(),
                "{ease:?}: from*to>0 and bounded → should be finite, got {result:?}"
            );
        }
    }

    #[test]
    fn back_overshoot_same_sign_safe_when_far() {
        // InBack overshoots to ≈ -0.1. With from=100, to=200, the critical
        // e value is from/(from-to) = 100/(100-200) = -1.0, well below -0.1,
        // so the overshoot never reaches the pole → safe.
        let result = evaluate_integration_reciprocal::<f64, 1>(
            AdvancedEasing::InBack,
            100.0,
            200.0,
            0.5,
            SAMPLES,
        );
        assert!(
            result.is_some(),
            "InBack with from=100, to=200: e_crit=-1.0 < e_min≈-0.1 → safe, got {result:?}"
        );
    }

    #[test]
    fn back_overshoot_same_sign_pole_when_close() {
        // InBack overshoots to ≈ -0.1. With from=1, to=1.05, the critical
        // e value is 1/(1-1.05) = 1/(-0.05) = -20, which is below -0.1 → safe.
        // But with from=1, to=10, e_crit = 1/(1-10) = -1/9 ≈ -0.111, which
        // is below e_min ≈ -0.1 → safe.
        // To actually hit the pole, we need e_crit ∈ [-0.1, 1.0], i.e.
        // from/(from-to) ∈ [-0.1, 1.0]. For from, to > 0 and to > from:
        //   e_crit = from/(from-to) < 0. Need e_crit >= -0.1.
        //   from/(from-to) >= -0.1  ⟺  from <= -0.1*(from-to) = 0.1*(to-from)
        //   ⟺  from <= 0.1*to - 0.1*from  ⟺  1.1*from <= 0.1*to  ⟺  to >= 11*from.
        // So with from=1, to=11: e_crit = 1/(1-11) = -0.1, right at the edge.
        // With from=1, to=12: e_crit = 1/(1-12) = -1/11 ≈ -0.0909, which is
        // inside [-0.1, 1.0] → pole!
        let result = evaluate_integration_reciprocal::<f64, 1>(
            AdvancedEasing::InBack,
            1.0,
            12.0,
            0.5,
            SAMPLES,
        );
        assert!(
            result.is_none(),
            "InBack with from=1, to=12: e_crit≈-0.09 ∈ [-0.1, 1.0] → pole, got {result:?}"
        );
    }

    #[test]
    fn hardcoded_overshoot_bounds_match_numeric() {
        // Verify the hardcoded constants match a fresh high-resolution scan.
        for &ease in AdvancedEasing::ALL {
            if !overshoots(ease) {
                continue;
            }
            const N: usize = 1_000_000;
            let mut e_min = f64::INFINITY;
            let mut e_max = f64::NEG_INFINITY;
            for k in 0..=N {
                let s = k as f64 / N as f64;
                let e = crate::math::easing::eval::easing_value(ease, s);
                if e < e_min {
                    e_min = e;
                }
                if e > e_max {
                    e_max = e;
                }
            }
            let (hardcoded_min, hardcoded_max) = overshoot_bounds(ease);
            // Hardcoded bounds must be at least as wide as the scan.
            assert!(
                hardcoded_min <= e_min + 1e-10,
                "{ease:?}: hardcoded e_min {hardcoded_min} > scanned {e_min}"
            );
            assert!(
                hardcoded_max >= e_max - 1e-10,
                "{ease:?}: hardcoded e_max {hardcoded_max} < scanned {e_max}"
            );
        }
    }

    #[test]
    fn opposite_sign_returns_none() {
        // from*to < 0 → denominator crosses zero → pole, for all non-Const easings.
        for &ease in AdvancedEasing::ALL {
            if matches!(
                ease,
                AdvancedEasing::InConst | AdvancedEasing::OutConst | AdvancedEasing::InOutConst
            ) {
                continue;
            }
            let result = evaluate_integration_reciprocal::<f64, 1>(ease, -1.0, 1.0, 0.5, SAMPLES);
            assert!(
                result.is_none(),
                "{ease:?}: from*to<0 → pole, got {result:?}"
            );
        }
    }
}