//! Single divstep transition for [`DivstepMatrix`].
use super::DivstepMatrix;

impl DivstepMatrix {
    /// Performs one divstep on the low bits `f0` and `g0`, updating the
    /// transition matrix and returning the new delta.
    ///
    /// The Bernstein-Yang divstep rule:
    /// - If `delta > 0` and `g` is odd: swap `(f, g)`, negate delta, then
    ///   `g ← (g + f) / 2`, `delta ← delta + 1`
    /// - Otherwise: `g ← (g + f * (g & 1)) / 2`, `delta ← delta + 1`
    ///   (when `g` is even, this is just `g ← g / 2`)
    ///
    /// All operations are branchless via arithmetic masks. The matrix rows
    /// track how `f` and `g` evolve so the same transform can be applied
    /// to the full-width values later.
    #[inline]
    pub(crate) fn step(&mut self, delta: i64, f0: &mut i64, g0: &mut i64) -> i64 {
        // cond = -1 if (delta > 0 && g is odd), else 0
        let cond = -(((delta > 0) as i64) & (*g0 & 1));

        // Conditional swap: if cond, swap (f, g) rows and negate delta
        let new_delta = (delta ^ cond) - cond; // if cond=-1: -delta, else delta

        // Conditionally swap f0/g0
        let old_f0 = *f0;
        let old_g0 = *g0;
        *f0 = (old_f0 & !cond) | (old_g0 & cond);
        *g0 = (old_g0 & !cond) | (old_f0 & cond);
        // Negate g0 if swapped (since after swap, g ← -old_f)
        // Actually in Bernstein-Yang: after swap, g is negated
        *g0 = (*g0 ^ cond) - cond;

        // Conditionally swap matrix rows
        let (old_u, old_v) = (self.u, self.v);
        let (old_q, old_r) = (self.q, self.r);
        self.u = (old_u & !cond) | (old_q & cond);
        self.v = (old_v & !cond) | (old_r & cond);
        self.q = (old_q & !cond) | (old_u & cond);
        self.r = (old_r & !cond) | (old_v & cond);
        // Negate q, r row if swapped
        self.q = (self.q ^ cond) - cond;
        self.r = (self.r ^ cond) - cond;

        // Now do the elimination step: g ← (g + f * (g & 1)) / 2
        let g_is_odd = -(*g0 & 1); // mask: -1 if g odd, 0 if even
        *g0 += *f0 & g_is_odd;
        self.q += self.u & g_is_odd;
        self.r += self.v & g_is_odd;

        // Divide g by 2 (arithmetic right shift since g is now even)
        *g0 >>= 1;
        // Scale u, v row by 2 instead of dividing q, r (keeps integers)
        self.u <<= 1;
        self.v <<= 1;

        new_delta + 1
    }
}

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

    /// Single step with delta=1, f=1, g=1 triggers swap path.
    #[test]
    fn step_swap_path() {
        let mut m = DivstepMatrix::identity();
        let mut f0 = 1i64;
        let mut g0 = 1i64;
        let new_delta = m.step(1, &mut f0, &mut g0);
        // After swap: f=g=1, g=-f+g=0, then g/=2=0, delta=-1+1=0
        assert_eq!(new_delta, 0);
        assert_eq!(g0, 0);
    }

    /// Single step with delta=0, f=1, g=3 does not swap.
    #[test]
    fn step_no_swap() {
        let mut m = DivstepMatrix::identity();
        let mut f0 = 1i64;
        let mut g0 = 3i64;
        let new_delta = m.step(0, &mut f0, &mut g0);
        // No swap (delta=0, not > 0). g is odd so g += f → 4, then g/=2 → 2
        assert_eq!(new_delta, 1);
        assert_eq!(f0, 1);
        assert_eq!(g0, 2);
    }

    /// Even g just halves without adding f.
    #[test]
    fn step_even_g() {
        let mut m = DivstepMatrix::identity();
        let mut f0 = 3i64;
        let mut g0 = 4i64;
        let new_delta = m.step(1, &mut f0, &mut g0);
        // g is even, no swap triggered. g/=2 → 2.
        assert_eq!(new_delta, 2);
        assert_eq!(f0, 3);
        assert_eq!(g0, 2);
    }

    /// 62 steps on small coprime values produce GCD=1.
    #[test]
    fn sixty_two_steps_small() {
        let mut m = DivstepMatrix::identity();
        let mut f0 = 7i64;
        let mut g0 = 5i64;
        let mut delta = 1i64;
        for _ in 0..62 {
            delta = m.step(delta, &mut f0, &mut g0);
        }
        // After enough steps, g should converge to 0 for small coprimes
        assert_eq!(g0, 0);
    }
}
