vsrg 0.3.0

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

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

    for i in 0..N {
        // Monotonic only when there is no sign flip and no zero endpoint (a zero
        // endpoint would make the reciprocal integrand diverge).
        if from[i] * to[i] <= 0.0 {
            return None;
        }
        if !from[i].is_finite() || !to[i].is_finite() || !val[i].is_finite() {
            return None;
        }
    }

    Some(T::compose(std::array::from_fn(|i| {
        let from = from[i];
        let to = to[i];
        let val = val[i];
        // R(-from, -to, p) = -R(from, to, p), so solving R(-from, -to, p) = -val
        // is equivalent to solving R(from, to, p) = val.
        if from < 0.0 || to < 0.0 {
            solve_integration_reciprocal(ease, -from, -to, -val, samples)
        } else {
            solve_integration_reciprocal(ease, from, to, val, samples)
        }
    })))
}

/// Solves `R(p) = val` for `p ∈ [0, 1]`, where
/// `R(p) = integrate_reciprocal(from, to, p) = ∫₀ᵖ 1/(from*(1-e(s)) + to*e(s)) ds`.
///
/// `from` and `to` are assumed to be same-signed and non-zero (the caller
/// rejects `from * to <= 0`), so the denominator is a convex combination of two
/// same-signed values and never vanishes; `R` is therefore strictly monotonic.
///
/// # Algorithm
/// `Linear` and the `*Const` variants admit closed-form inverses and are
/// evaluated analytically. All other easings are solved by a **single forward
/// walk**: step `s` over `[0, 1]` on a fixed grid of `samples` panels,
/// accumulate the reciprocal integral one Simpson panel-pair at a time
/// (`1/ease(from, to, s)` is the integrand, by the FTC), and stop as soon as
/// the running sum reaches `val`. The crossing within the final pair is found
/// by linearly interpolating the accumulated antiderivative. This is a single
/// O(`samples`) pass (with early termination), in contrast to a Newton/bisection
/// loop that would re-evaluate the whole integral on every iteration.
///
/// # Accuracy
/// The completed panel-pairs use Simpson's rule (4th-order), but the crossing
/// within the final partial pair uses linear interpolation (2nd-order), so the
/// overall error is O(1/`samples`²) — roughly `1e-6` at the default 1024
/// samples.
///
/// This walk integrates on a **fixed `[0, 1]` grid**, whereas
/// [`Easing::integrate_reciprocal`] uses a grid compressed to `[0, p]`
/// (`h = p/samples`). The two quadratures therefore differ slightly, so solving
/// and then re-integrating does not round-trip bit-exactly: the discrepancy is
/// of the same order as the quadrature error (≈ `1e-6` at 1024 samples). This
/// is the accepted tradeoff for the ~24× speedup over an exact Newton/bisection
/// inversion.
fn solve_integration_reciprocal(ease: Easing, from: f64, to: f64, val: f64, samples: usize) -> f64 {
    if val <= 0.0 {
        return 0.0;
    }

    match ease {
        Easing::Linear => {
            // R(p) = ln(denom_p / from) / (to - from), denom_p = from*(1-p) + to*p.
            let delta = to - from;
            if delta.abs() <= 1e-12 {
                // from == to: R(p) = p / from.
                return (val * from).clamp(0.0, 1.0);
            }
            let denom_p = from * (val * delta).exp();
            let p = (denom_p - from) / delta;
            p.clamp(0.0, 1.0)
        }
        Easing::InConst => {
            // R(p) = p / from.
            (val * from).clamp(0.0, 1.0)
        }
        Easing::OutConst => {
            // R(p) = p / to.
            (val * to).clamp(0.0, 1.0)
        }
        Easing::InOutConst => {
            // R(p) = p/from for p <= 0.5; 0.5/from + (p-0.5)/to for p > 0.5.
            let half_r = 0.5 / from;
            if val <= half_r {
                (val * from).clamp(0.0, 0.5)
            } else {
                0.5 + ((val - half_r) * to).clamp(0.0, 0.5)
            }
        }

        // All remaining easings: single-pass Simpson walk.
        _ => walk_solve_reciprocal(ease.as_easing(), from, to, val, samples),
    }
}

/// Single-pass solver: accumulate `∫ 1/ease(from, to, s) ds` one Simpson
/// panel-pair at a time over a fixed `[0, 1]` grid, stopping at the pair that
/// contains the crossing of `target`. Returns the `s ∈ [0, 1]` where the
/// running integral reaches `target`.
///
/// See [`solve_integration_reciprocal`] for accuracy and round-trip caveats.
fn walk_solve_reciprocal(
    ease: AdvancedEasing,
    from: f64,
    to: f64,
    target: f64,
    samples: usize,
) -> f64 {
    // Simpson needs an even number of panels; pad to at least 2.
    let n = if samples < 2 {
        2
    } else if samples % 2 == 1 {
        samples + 1
    } else {
        samples
    };
    let h = 1.0 / n as f64;
    // The integrand of the reciprocal integral is 1 / ease(from, to, s).
    let integrand = |s: f64| {
        let e = easing_value(ease, s);
        1.0 / (from * (1.0 - e) + to * e)
    };

    let npairs = n / 2;
    let mut cum = 0.0_f64; // Simpson integral from 0 to the start of the current pair.
    let mut f_left = integrand(0.0);

    for k in 0..npairs {
        let f_mid = integrand(((2 * k + 1) as f64) * h);
        let f_right = integrand(((2 * k + 2) as f64) * h);
        // Simpson over [s0, s2] = [2k*h, (2k+2)*h], width 2h.
        let pair = (h / 3.0) * (f_left + 4.0 * f_mid + f_right);

        if cum + pair >= target {
            // The crossing lies within this pair. Linearly interpolate the
            // accumulated antiderivative between (s0, cum) and (s2, cum+pair).
            let frac = if pair > 0.0 && pair.is_finite() {
                ((target - cum) / pair).clamp(0.0, 1.0)
            } else {
                0.0
            };
            let s0 = (2 * k) as f64 * h;
            let s2 = (2 * k + 2) as f64 * h;
            return (s0 + frac * (s2 - s0)).clamp(0.0, 1.0);
        }

        cum += pair;
        f_left = f_right;
    }

    // `target` is at or past R(1); clamp to the end.
    1.0
}