//! The Euclidean algorithm — gcd, lcm, and the extended (Bézout) form.

/// The greatest common divisor of `a` and `b`, by Euclid's algorithm: repeatedly
/// replace `(a, b)` with `(b, a mod b)` until the remainder is zero. `gcd(n, 0)`
/// is `n`; `gcd(0, 0)` is `0`. No overflow — every remainder is smaller than its
/// divisor.
///
/// # Examples
///
/// ```
/// use coprime::number::gcd;
///
/// assert_eq!(gcd(48, 36), 12);
/// assert_eq!(gcd(17, 5), 1); // coprime
/// assert_eq!(gcd(7, 0), 7); // gcd(n, 0) == n
/// assert_eq!(gcd(0, 0), 0);
/// ```
#[must_use]
pub fn gcd(mut a: u64, mut b: u64) -> u64 {
    while b != 0 {
        (a, b) = (b, a % b);
    }
    a
}

/// The least common multiple of `a` and `b`: the smallest non-negative integer
/// divisible by both. Computed as `a / gcd(a, b) * b` — dividing *before*
/// multiplying, so the only place the result can overflow is the true answer.
/// `lcm(a, 0) = lcm(0, b) = 0`. `None` when the result exceeds `u64::MAX`.
///
/// # Examples
///
/// ```
/// use coprime::number::lcm;
///
/// assert_eq!(lcm(4, 6), Some(12));
/// assert_eq!(lcm(21, 6), Some(42));
/// assert_eq!(lcm(7, 0), Some(0));
/// assert_eq!(lcm(1u64 << 62, (1u64 << 62) + 2), None); // overflows u64
/// ```
#[must_use]
pub fn lcm(a: u64, b: u64) -> Option<u64> {
    if a == 0 || b == 0 {
        return Some(0);
    }
    // gcd(a, b) divides a, so `a / gcd` is exact; only the product can overflow.
    (a / gcd(a, b)).checked_mul(b)
}

/// Euclid's algorithm extended to the Bézout identity: returns `(d, s, t)` with
/// `d == gcd(a, b)` and `s*a + t*b == d`. The coefficients `s, t` are signed.
///
/// Correct for every `u64` input. The accumulators are `i128`: the gcd is
/// non-negative and fits `u64`, while the signed Bézout cofactors can exceed
/// `i64`. No `i128` product overflows — for `a, b > 1` the returned cofactors
/// satisfy `|s| <= b/2` and `|t| <= a/2` (each grows only as its paired remainder
/// shrinks), so `|s*a|, |t*b| <= a*b/2 < 2^127`; if `a` or `b` is `0` or `1` a
/// cofactor is `0` or `1` and the product is at most `max(a, b) < 2^64`. The two
/// products have opposite signs and sum to `d`. At the `u64` ceiling this bound
/// is *tight*: `(2^64 - 1)^2 / 2` lands within one `u64` step of `i128::MAX`, so
/// the accumulators carry no spare headroom — widening the operand type past
/// `u64` would require re-deriving the bound, not just reusing `i128`.
///
/// # Examples
///
/// ```
/// use coprime::number::extended_gcd;
///
/// let (d, s, t) = extended_gcd(240, 46);
/// assert_eq!(d, 2); // gcd(240, 46)
/// assert_eq!(s * 240 + t * 46, i128::from(d)); // the Bézout identity
/// ```
#[must_use]
#[allow(
    clippy::cast_possible_truncation, // d = gcd(a,b) < 2^64 fits u64
    clippy::cast_sign_loss,           // d >= 0
    clippy::many_single_char_names    // a, b, q, r, s, t are the canonical Bézout names
)]
pub fn extended_gcd(a: u64, b: u64) -> (u64, i128, i128) {
    let (mut old_r, mut r) = (i128::from(a), i128::from(b));
    let (mut old_s, mut s) = (1i128, 0i128);
    let (mut old_t, mut t) = (0i128, 1i128);
    while r != 0 {
        let q = old_r / r;
        (old_r, r) = (r, old_r - q * r);
        (old_s, s) = (s, old_s - q * s);
        (old_t, t) = (t, old_t - q * t);
    }
    // old_r >= 0: a, b >= 0 and every remainder is non-negative; d <= max(a, b) < 2^64.
    (old_r as u64, old_s, old_t)
}
