//! Modular arithmetic — exponentiation, inversion, and the u128 mulmod helper.

use crate::number::extended_gcd;

/// Modular exponentiation `base^exp mod modulus`, by square-and-multiply: walk
/// the bits of `exp` from low to high, folding the running base into the
/// accumulator on a set bit and squaring it for the next bit — `O(log exp)`
/// multiplications instead of `exp` of them. The squaring is guarded so the
/// final bit does not square a value that is never used.
///
/// Total and panic-free over every `u64` input. Each product is formed in `u128`
/// — two values below `2^64` multiply below `2^128` — so nothing overflows, and
/// the reduction keeps both operands below `modulus`. `modulus <= 1` returns `0`:
/// every value is `0 mod 1`, and the degenerate `mod 0` is mapped to `0` rather
/// than trapping. Callers wanting the ring `Z/mZ` pass `modulus >= 2`.
///
/// # Examples
///
/// ```
/// use coprime::number::modpow;
///
/// assert_eq!(modpow(2, 10, 1000), 24); // 1024 mod 1000
/// assert_eq!(modpow(4, 13, 497), 445); // a textbook RSA encryption step
/// assert_eq!(modpow(7, 0, 13), 1); // x^0 == 1
/// assert_eq!(modpow(5, 3, 1), 0); // everything is 0 mod 1
/// ```
#[must_use]
#[allow(clippy::cast_possible_truncation)] // result < modulus <= u64::MAX fits u64
pub fn modpow(base: u64, exp: u64, modulus: u64) -> u64 {
    if modulus <= 1 {
        return 0;
    }
    let m = u128::from(modulus);
    let mut result: u128 = 1;
    let mut b = u128::from(base) % m;
    let mut e = exp;
    while e != 0 {
        if e & 1 == 1 {
            result = result * b % m;
        }
        e >>= 1;
        if e != 0 {
            b = b * b % m;
        }
    }
    result as u64
}

/// The modular inverse of `a` modulo `modulus`: the unique `x` in `[0, modulus)`
/// with `a * x ≡ 1 (mod modulus)`, or `None` when no inverse exists — that is,
/// when `gcd(a, modulus) != 1`. `modulus <= 1` is the degenerate ring and has no
/// units, so it returns `None`.
///
/// Built on [`extended_gcd`]: `extended_gcd(a mod modulus, modulus)` returns
/// `(d, s, _)` with `s*(a mod modulus) + _*modulus = d`. When `d == 1`, `s` is an
/// inverse of `a` modulo `modulus`; it is a signed Bézout cofactor, so it is
/// reduced into the canonical residue `[0, modulus)`.
///
/// # Examples
///
/// ```
/// use coprime::number::modinv;
///
/// assert_eq!(modinv(3, 11), Some(4)); // 3 * 4 == 12 ≡ 1 (mod 11)
/// assert_eq!(modinv(2, 4), None); // gcd(2, 4) != 1, so no inverse
/// assert_eq!(modinv(5, 1), None); // the degenerate ring has no units
/// ```
#[must_use]
#[allow(
    clippy::cast_possible_truncation, // inv in [0, modulus) fits u64
    clippy::cast_sign_loss            // inv >= 0 after reduction
)]
pub fn modinv(a: u64, modulus: u64) -> Option<u64> {
    if modulus <= 1 {
        return None;
    }
    let (d, s, _) = extended_gcd(a % modulus, modulus);
    if d != 1 {
        return None;
    }
    let m = i128::from(modulus);
    let inv = ((s % m) + m) % m; // canonical residue in [0, modulus)
    Some(inv as u64)
}

/// Modular multiplication `a * b mod m`, the product formed in `u128` so it
/// never overflows: two values below `2^64` multiply below `2^128`, and the
/// reduction brings the result back under `m`. This is the widening idiom the
/// whole product path rests on — `sqrtmod` and the Pollard–Brent splitter call
/// this helper, while [`modpow`] deliberately inlines the same shape (its
/// operands are already `u128`-widened, so routing through here would round-trip
/// through `u64`); the Stage Two Montgomery swap is thus a two-site change
/// (here and in `modpow`).
///
/// Total and panic-free over every `u64` input. `m <= 1` returns `0`: every
/// value is `0 mod 1`, and the degenerate `mod 0` is mapped to `0` rather than
/// trapping — the same convention as [`modpow`]. Callers wanting the ring
/// `Z/mZ` pass `m >= 2`.
///
/// # Examples
///
/// ```
/// use coprime::number::mulmod;
///
/// assert_eq!(mulmod(7, 8, 5), 1); // 56 mod 5
/// assert_eq!(mulmod(u64::MAX, u64::MAX, 1_000_000_007), 114_944_269); // no overflow
/// assert_eq!(mulmod(3, 4, 1), 0); // everything is 0 mod 1
/// ```
#[must_use]
#[allow(clippy::cast_possible_truncation)] // a*b % m < m <= u64::MAX fits u64
pub fn mulmod(a: u64, b: u64, m: u64) -> u64 {
    if m <= 1 {
        return 0;
    }
    (u128::from(a) * u128::from(b) % u128::from(m)) as u64
}
