//! Arithmetic modulo the order `L` of the prime-order subgroup generated by the Ed25519 base
//! point, `L = 2^252 + 27742317777372353535851937790883648493`.

use subtle::{Choice, ConditionallySelectable};
use zeroize::Zeroize;

/// `L`, as four little-endian 64-bit limbs.
const L: [u64; 4] = [
    0x5812631a5cf5d3ed,
    0x14def9dea2f79cd6,
    0x0000000000000000,
    0x1000000000000000,
];

/// `floor(2^512 / L)`, the Barrett reduction constant for reducing a 512-bit value modulo `L`
/// (see [`barrett_reduce`]). `L` is just over `2^252`, so this is just under `2^260`, five
/// little-endian 64-bit limbs.
const MU: [u64; 5] = [
    0xed9ce5a30a2c131b,
    0x2106215d086329a7,
    0xffffffffffffffeb,
    0xffffffffffffffff,
    0x000000000000000f,
];

/// An integer modulo `L`, always canonically reduced (`< L`).
#[derive(Copy, Clone, Zeroize)]
pub struct Scalar([u64; 4]);

/// Returns `true` if `a < b`, comparing as 256-bit unsigned integers (little-endian limbs).
fn limbs_lt(a: &[u64; 4], b: &[u64; 4]) -> bool {
    for i in (0..4).rev() {
        if a[i] != b[i] {
            return a[i] < b[i];
        }
    }
    false
}

/// Returns `a - b` modulo `2^(64*N)` and whether the subtraction underflowed.
fn limbs_sub_with_borrow<const N: usize>(a: &[u64; N], b: &[u64; N]) -> ([u64; N], bool) {
    let mut out = [0u64; N];
    let mut borrow = false;
    for i in 0..N {
        let (d1, b1) = a[i].overflowing_sub(b[i]);
        let (d2, b2) = d1.overflowing_sub(borrow as u64);
        out[i] = d2;
        borrow = b1 | b2;
    }
    (out, borrow)
}

/// Returns `a - b`, assuming `a >= b`.
fn limbs_sub(a: &[u64; 4], b: &[u64; 4]) -> [u64; 4] {
    limbs_sub_with_borrow(a, b).0
}

/// Returns `a + b`, assuming the true sum is `< 2^256` (no carry out of the top limb).
fn limbs_add(a: &[u64; 4], b: &[u64; 4]) -> [u64; 4] {
    let mut out = [0u64; 4];
    let mut carry = false;
    for i in 0..4 {
        let (sum, c1) = a[i].overflowing_add(b[i]);
        let (sum, c2) = sum.overflowing_add(carry as u64);
        out[i] = sum;
        carry = c1 | c2;
    }
    out
}

/// Selects `b` when `select_b` is set, in constant time via `subtle`'s optimization barrier.
fn limbs_select<const N: usize>(a: &[u64; N], b: &[u64; N], select_b: Choice) -> [u64; N] {
    core::array::from_fn(|i| u64::conditional_select(&a[i], &b[i], select_b))
}

/// Subtracts `b` from `a` when `a >= b`, leaving `a` unchanged otherwise.
fn limbs_conditional_sub<const N: usize>(a: &[u64; N], b: &[u64; N]) -> [u64; N] {
    let (difference, underflow) = limbs_sub_with_borrow(a, b);
    limbs_select(&difference, a, Choice::from(underflow as u8))
}

/// Returns `(a - b) mod 2^320`, as five little-endian 64-bit limbs: the final borrow (if any) is
/// discarded, which is exactly right both when `a >= b` (ordinary subtraction, no borrow) and,
/// in [`barrett_reduce`], when the true difference is negative (the discarded borrow is
/// equivalent to adding `2^320` back in).
fn limbs_sub5(a: &[u64; 5], b: &[u64; 5]) -> [u64; 5] {
    limbs_sub_with_borrow(a, b).0
}

/// Returns the full product `a * b` as little-endian 64-bit limbs, via the standard
/// schoolbook multiply-accumulate-with-carry ("Comba") method.
fn limbs_mul_wide<const A: usize, const B: usize, const OUT: usize>(
    a: &[u64; A],
    b: &[u64; B],
) -> [u64; OUT] {
    const { assert!(OUT == A + B) };
    let mut t = [0u64; OUT];
    for i in 0..A {
        let mut carry = 0u128;
        for j in 0..B {
            let sum = t[i + j] as u128 + (a[i] as u128) * (b[j] as u128) + carry;
            t[i + j] = sum as u64;
            carry = sum >> 64;
        }
        t[i + B] = carry as u64;
    }
    t
}

/// Reduces a 512-bit value (eight little-endian 64-bit limbs) modulo `L` with a fixed number of
/// limb multiplications using Barrett reduction (Handbook of Applied Cryptography, Algorithm
/// 14.42).
///
/// With word size `b = 2^64` and `L` fitting in `k = 4` words, this reduces any `x < b^(2k) =
/// 2^512`: `q = floor(floor(x / b^(k-1)) * MU / b^(k+1))` approximates `floor(x / L)` (using the
/// precomputed `MU = floor(b^(2k) / L)`), accurate enough that `x - q*L` is guaranteed `< 3L`, so
/// at most two trial subtractions of `L` remain to reach the canonical residue.
fn barrett_reduce(x: [u64; 8]) -> Scalar {
    let q1: [u64; 5] = x[3..8].try_into().expect("slice has 5 elements");
    let q2 = limbs_mul_wide::<5, 5, 10>(&q1, &MU);
    let q3: [u64; 5] = q2[5..10].try_into().expect("slice has 5 elements");

    let r1: [u64; 5] = x[0..5].try_into().expect("slice has 5 elements");
    let r2: [u64; 5] = limbs_mul_wide::<5, 4, 9>(&q3, &L)[0..5]
        .try_into()
        .expect("slice has 5 elements");
    let mut r = limbs_sub5(&r1, &r2);

    let l5 = [L[0], L[1], L[2], L[3], 0];
    // The Barrett quotient estimate leaves a residue below `3L`, so two fixed conditional
    // subtractions are sufficient and avoid data-dependent control flow for secret scalars.
    r = limbs_conditional_sub(&r, &l5);
    r = limbs_conditional_sub(&r, &l5);
    debug_assert_eq!(r[4], 0, "residue must fit in L's 4 limbs after reduction");
    Scalar([r[0], r[1], r[2], r[3]])
}

impl Scalar {
    pub const ZERO: Self = Self([0, 0, 0, 0]);

    /// Returns the additive inverse modulo `L`.
    pub fn neg_mod_l(&self) -> Self {
        Self(limbs_conditional_sub(&limbs_sub(&L, &self.0), &L))
    }

    /// Interprets `bytes` as a little-endian integer and rejects it unless it is already the
    /// canonical representative (`< L`), as required for the `s` component of a signature.
    pub fn from_canonical_bytes(bytes: &[u8; 32]) -> Option<Self> {
        let mut limbs = [0u64; 4];
        for (i, limb) in limbs.iter_mut().enumerate() {
            let mut chunk = [0u8; 8];
            chunk.copy_from_slice(&bytes[i * 8..i * 8 + 8]);
            *limb = u64::from_le_bytes(chunk);
        }
        limbs_lt(&limbs, &L).then_some(Self(limbs))
    }

    /// Reduces a 64-byte little-endian integer (e.g. a SHA-512 digest) modulo `L`, via
    /// [`barrett_reduce`].
    pub fn from_bytes_mod_order_wide(bytes: &[u8; 64]) -> Self {
        let limbs = core::array::from_fn(|i| {
            let mut chunk = [0u8; 8];
            chunk.copy_from_slice(&bytes[i * 8..i * 8 + 8]);
            u64::from_le_bytes(chunk)
        });
        barrett_reduce(limbs)
    }

    /// Returns the bits of this scalar's canonical representative, most significant first.
    #[cfg(test)]
    pub fn bits_be(&self) -> impl Iterator<Item = bool> + '_ {
        (0..256)
            .rev()
            .map(move |i| (self.0[i / 64] >> (i % 64)) & 1 == 1)
    }

    /// Returns the canonical little-endian encoding of this scalar.
    pub fn to_bytes(self) -> [u8; 32] {
        let mut bytes = [0u8; 32];
        for (chunk, limb) in bytes.as_chunks_mut::<8>().0.iter_mut().zip(self.0) {
            *chunk = limb.to_le_bytes();
        }
        bytes
    }

    /// Constructs a scalar from a little-endian 128-bit value. Always canonical, since every
    /// 128-bit value is `< L` (`L > 2^252`).
    pub const fn from_u128(value: u128) -> Self {
        Self([value as u64, (value >> 64) as u64, 0, 0])
    }

    /// Returns the base-`2^width` digit at position `index`, i.e. bits `[index*width,
    /// index*width+width)` of this scalar's canonical representative, as an unsigned integer.
    ///
    /// # Panics
    ///
    /// Panics if `width` is 32 or more.
    pub const fn window(&self, index: usize, width: u32) -> u32 {
        assert!(width < 32);
        let bit_start = index * width as usize;
        if bit_start >= 256 {
            return 0;
        }
        let limb_index = bit_start / 64;
        let bit_offset = (bit_start % 64) as u32;
        let mut digit = self.0[limb_index] >> bit_offset;
        if bit_offset + width > 64 && limb_index + 1 < 4 {
            digit |= self.0[limb_index + 1] << (64 - bit_offset);
        }
        (digit & ((1u64 << width) - 1)) as u32
    }

    /// Recodes this scalar into `N` signed, base-`2^width` digits (each in
    /// `[-2^(width-1), 2^(width-1) - 1]`), least-significant first. Used by the MSM's bucket
    /// method (see [`super::msm`]): a signed digit only ever needs a bucket for its *magnitude*
    /// (`1..=2^(width-1)`, half as many as the `1..2^width` an unsigned digit needs), at the cost
    /// of negating the affine point's X and T coordinates when its digit is negative.
    ///
    /// This is exactly the standard signed-digit recoding: process [`Scalar::window`]'s unsigned
    /// digits from least to most significant, and whenever one is `>= 2^(width-1)` (the upper half
    /// of its range), replace it with `digit - 2^width` (negative, same residue) and carry a `+1`
    /// into the next digit, since `digit * 2^(width*i) = (digit - 2^width) * 2^(width*i) +
    /// 2^width * 2^(width*i)`, and that `2^width` term is exactly one unit of the next digit's
    /// weight. `N` must leave room for the final carry (at most `1`). A scalar is below
    /// `L < 2^253`, so for any `width` of at least 2 the top digit of `256usize.div_ceil(width)`
    /// digits absorbs it.
    pub fn signed_digits<const N: usize>(&self, width: u32) -> [i32; N] {
        let half = 1i64 << (width - 1);
        let mut digits = [0i32; N];
        let mut carry = 0i64;
        let mut i = 0;
        while i < N {
            // `raw` is at most `2^width`, so the carry is 1 exactly when `raw >= half`. Random
            // digits make a branch on that comparison unpredictable.
            let raw = i64::from(self.window(i, width)) + carry;
            carry = (raw + half) >> width;
            digits[i] = (raw - (carry << width)) as i32;
            i += 1;
        }
        digits
    }

    /// Returns `(self + rhs) mod L`, assuming both operands are already `< L`.
    pub fn add_mod_l(&self, rhs: &Self) -> Self {
        let sum = limbs_add(&self.0, &rhs.0);
        Self(limbs_conditional_sub(&sum, &L))
    }

    /// Returns `(self * rhs) mod L`.
    ///
    /// Computes the full 512-bit product via schoolbook multiplication, then reduces it via
    /// [`barrett_reduce`].
    pub fn mul_mod_l(&self, rhs: &Self) -> Self {
        barrett_reduce(limbs_mul_wide(&self.0, &rhs.0))
    }

    /// Recodes this scalar into width-`W` non-adjacent form, least significant digit first.
    ///
    /// The digits sum to this scalar with weights `2^i`. Every nonzero digit is odd with
    /// magnitude below `2^(W-1)`, and any `W` consecutive digits hold at most one nonzero digit.
    /// This operation is variable-time, so the scalar must be public.
    pub const fn naf<const W: usize>(&self) -> [i8; 256] {
        const { assert!(2 <= W && W <= 8) };

        // A zero fifth limb lets windows starting near bit 255 read past it. A carry into the
        // window at `position` is a pending `+1` at that bit. Scalars are below `2^253`, so
        // every carry is absorbed by bit 255.
        let limbs = [self.0[0], self.0[1], self.0[2], self.0[3], 0];
        let width = 1u64 << W;
        let mut digits = [0i8; 256];
        let mut carry = 0;
        let mut position = 0;
        while position < 256 {
            // The `W` bits at `position`, which can straddle two limbs, plus the pending carry.
            let (index, bit) = (position / 64, position % 64);
            let mut bits = limbs[index] >> bit;
            if bit + W > 64 {
                bits |= limbs[index + 1] << (64 - bit);
            }
            let window = (bits & (width - 1)) + carry;

            // An even window makes a zero digit, and the carry moves up with the position. An odd
            // window becomes the odd digit of its residue mod `2^W` with magnitude below
            // `2^(W-1)`. That digit clears the window, so the next `W - 1` digits are zero, and a
            // negative digit carries `+1` into bit `position + W`.
            if window & 1 == 0 {
                position += 1;
                continue;
            }
            if window < width / 2 {
                carry = 0;
                digits[position] = window as i8;
            } else {
                carry = 1;
                digits[position] = (window as i16 - width as i16) as i8;
            }
            position += W;
        }
        digits
    }

    /// Returns the low and high 128-bit halves of this scalar's canonical representative.
    pub const fn halves(&self) -> (u128, u128) {
        (
            self.0[0] as u128 | (self.0[1] as u128) << 64,
            self.0[2] as u128 | (self.0[3] as u128) << 64,
        )
    }

    /// Returns `(u, v)` with `v = u*self (mod L)`, `0 < |u| < 2^127`, and `v < 2^126`.
    ///
    /// Runs the extended Euclidean algorithm on `(L, self)` and stops at the first remainder
    /// below `2^126`, the half-size decomposition of Antipa et al., "Accelerated Verification of
    /// ECDSA Signatures" (SAC 2005). Variable-time, so `self` must be public.
    pub fn half_size(&self) -> (i128, u128) {
        // Each remainder `r` is paired with its coefficient `t`, keeping `r = t*self (mod L)`.
        // Coefficients alternate in sign and satisfy `|t1|*r0 <= L` for consecutive pairs, so
        // stopping once `r1 < 2^126 <= r0` gives `|t1| <= L / 2^126 < 2^127`. Every partial
        // coefficient lies between two consecutive ones, so `i128` never overflows.
        let (mut r0, mut t0) = (L, 0i128);
        let (mut r1, mut t1) = (self.0, 1i128);
        while limbs_bits(&r1) > 126 {
            // Reduces `r0` modulo `r1` by subtracting shifted copies of `r1`. Since
            // `r0 < 2^253` and `r1 >= 2^126`, every shift is below 127.
            while !limbs_lt(&r0, &r1) {
                let mut shift = limbs_bits(&r0) - limbs_bits(&r1);
                let mut shifted = limbs_shl(&r1, shift);
                if limbs_lt(&r0, &shifted) {
                    shift -= 1;
                    shifted = limbs_shl(&r1, shift);
                }
                r0 = limbs_sub(&r0, &shifted);
                t0 -= t1 << shift;
            }
            core::mem::swap(&mut r0, &mut r1);
            core::mem::swap(&mut t0, &mut t1);
        }
        (t1, u128::from(r1[0]) | u128::from(r1[1]) << 64)
    }
}

/// Returns the bit length of a 256-bit integer.
fn limbs_bits(a: &[u64; 4]) -> u32 {
    for i in (0..4).rev() {
        if a[i] != 0 {
            return 64 * i as u32 + 64 - a[i].leading_zeros();
        }
    }
    0
}

/// Returns `a << shift` for `shift < 256`, discarding bits above `2^256`.
fn limbs_shl(a: &[u64; 4], shift: u32) -> [u64; 4] {
    // Whole limbs select the source word; a nonzero bit offset also carries the preceding
    // word's high bits into it. The zero-offset case avoids a shift by 64.
    let words = (shift / 64) as usize;
    let bits = shift % 64;
    let mut out = [0u64; 4];
    for i in words..4 {
        out[i] = a[i - words] << bits;
        if bits != 0 && i > words {
            out[i] |= a[i - words - 1] >> (64 - bits);
        }
    }
    out
}

#[cfg(test)]
impl<'a> arbitrary::Arbitrary<'a> for Scalar {
    fn arbitrary(u: &mut arbitrary::Unstructured<'a>) -> arbitrary::Result<Self> {
        Ok(match u.int_in_range(0u8..=5)? {
            0 => Self::ZERO,
            1 => Self::from_u128(1),
            2 => {
                let mut l_minus_one = L;
                l_minus_one[0] -= 1;
                Self(l_minus_one)
            }
            3 => Self::from_u128(u.arbitrary()?),
            _ => Self::from_bytes_mod_order_wide(&u.arbitrary()?),
        })
    }
}

#[cfg(test)]
mod tests {
    use super::{L, Scalar, limbs_lt, limbs_sub};
    use commonware_invariants::minifuzz::Builder;
    use zeroize::Zeroize;

    /// Returns `(2*a + bit) mod 2^256`, the left-shift-by-one step of the bit-serial reduction
    /// [`reduce_wide_naive`] uses.
    fn limbs_shl1(a: &[u64; 4], bit: u64) -> [u64; 4] {
        let mut out = [0u64; 4];
        let mut carry = bit;
        for i in 0..4 {
            out[i] = (a[i] << 1) | carry;
            carry = a[i] >> 63;
        }
        out
    }

    /// A bit-serial double-and-add-mod reduction used as a differential-test oracle for
    /// [`Scalar::from_bytes_mod_order_wide`]. It performs one conditional subtraction per input
    /// bit, making it slow but straightforward to audit.
    fn reduce_wide_naive(bytes: &[u8; 64]) -> Scalar {
        let mut r = [0u64; 4];
        for byte_index in (0..64).rev() {
            let byte = bytes[byte_index];
            for bit_index in (0..8).rev() {
                let bit = (byte >> bit_index) & 1;
                r = limbs_shl1(&r, bit as u64);
                if !limbs_lt(&r, &L) {
                    r = limbs_sub(&r, &L);
                }
            }
        }
        Scalar(r)
    }

    #[test]
    fn from_bytes_mod_order_wide_matches_naive_reference() {
        Builder::default()
            .with_seed(0)
            .with_search_limit(1024)
            .test(|u| {
                let bytes: [u8; 64] = u.arbitrary()?;
                let expected = reduce_wide_naive(&bytes);
                let actual = Scalar::from_bytes_mod_order_wide(&bytes);
                assert_eq!(actual.0, expected.0);
                Ok(())
            });
    }

    #[test]
    fn from_bytes_mod_order_wide_matches_naive_reference_on_edge_cases() {
        // `x = 0`, `x = L - 1`, `x = L`, `x = L + 1`, and `x = 2^512 - 1` (the Barrett quotient
        // approximation's error is largest for inputs near the top of the input range).
        let l_minus_1 = {
            let mut bytes = [0u8; 64];
            for (i, limb) in L.iter().enumerate() {
                bytes[i * 8..i * 8 + 8].copy_from_slice(&(limb - u64::from(i == 0)).to_le_bytes());
            }
            bytes
        };
        let mut cases = vec![[0u8; 64], l_minus_1];
        let mut l_bytes = [0u8; 64];
        for (i, limb) in L.iter().enumerate() {
            l_bytes[i * 8..i * 8 + 8].copy_from_slice(&limb.to_le_bytes());
        }
        cases.push(l_bytes);
        let mut l_plus_1 = l_bytes;
        l_plus_1[0] = l_plus_1[0].wrapping_add(1);
        cases.push(l_plus_1);
        cases.push([0xffu8; 64]);

        for bytes in cases {
            let expected = reduce_wide_naive(&bytes);
            let actual = Scalar::from_bytes_mod_order_wide(&bytes);
            assert_eq!(actual.0, expected.0);
        }
    }

    #[test]
    fn mul_mod_l_matches_naive_reference() {
        Builder::default()
            .with_seed(0)
            .with_search_limit(1024)
            .test(|u| {
                let a: Scalar = u.arbitrary()?;
                let b: Scalar = u.arbitrary()?;

                let wide: [u64; 8] = super::limbs_mul_wide(&a.0, &b.0);
                let mut bytes = [0u8; 64];
                for (i, limb) in wide.iter().enumerate() {
                    bytes[i * 8..i * 8 + 8].copy_from_slice(&limb.to_le_bytes());
                }
                let expected = reduce_wide_naive(&bytes);

                assert_eq!(a.mul_mod_l(&b).0, expected.0);
                Ok(())
            });
    }

    #[test]
    fn neg_mod_l_is_canonical_additive_inverse() {
        for scalar in [
            Scalar::ZERO,
            Scalar::from_u128(1),
            Scalar(limbs_sub(&L, &[1, 0, 0, 0])),
        ] {
            let negative = scalar.neg_mod_l();
            assert!(limbs_lt(&negative.0, &L));
            assert_eq!(scalar.add_mod_l(&negative).0, Scalar::ZERO.0);
            assert_eq!(negative.neg_mod_l().0, scalar.0);
        }

        Builder::default()
            .with_seed(0)
            .with_search_limit(64)
            .test(|u| {
                let scalar: Scalar = u.arbitrary()?;
                let negative = scalar.neg_mod_l();
                assert!(limbs_lt(&negative.0, &L));
                assert_eq!(scalar.add_mod_l(&negative).0, Scalar::ZERO.0);
                assert_eq!(negative.neg_mod_l().0, scalar.0);
                Ok(())
            });
    }

    #[test]
    fn signed_digits_reconstruct_value() {
        const N: usize = 256usize.div_ceil(6) + 1;

        Builder::default()
            .with_seed(0)
            .with_search_limit(64)
            .test(|u| {
                let s: Scalar = u.arbitrary()?;
                for width in 6..=10 {
                    let digits = s.signed_digits::<N>(width);
                    let base = Scalar::from_u128(1u128 << width);
                    let mut reconstructed = Scalar::ZERO;
                    let mut power = Scalar::from_u128(1);
                    for &digit in &digits {
                        let magnitude = Scalar::from_u128(digit.unsigned_abs() as u128);
                        let term = power.mul_mod_l(&magnitude);
                        let term = if digit < 0 { term.neg_mod_l() } else { term };
                        reconstructed = reconstructed.add_mod_l(&term);
                        power = power.mul_mod_l(&base);
                    }
                    assert_eq!(reconstructed.0, s.0, "width={width}");
                }
                Ok(())
            });
    }

    #[test]
    fn signed_digits_are_in_range() {
        const N: usize = 256usize.div_ceil(6) + 1;

        Builder::default()
            .with_seed(0)
            .with_search_limit(64)
            .test(|u| {
                let s: Scalar = u.arbitrary()?;
                for width in 6..=10 {
                    let half = 1i32 << (width - 1);
                    let digits = s.signed_digits::<N>(width);
                    for digit in digits {
                        assert!((-half..half).contains(&digit));
                    }
                    let windows = 256usize.div_ceil(width as usize) + 1;
                    assert!(digits[windows..].iter().all(|&digit| digit == 0));
                }
                Ok(())
            });
    }

    /// Checks that `scalar`'s width-`W` non-adjacent form reconstructs it, and that its nonzero
    /// digits are odd, below `2^(W-1)` in magnitude, and at least `W` positions apart.
    fn check_naf<const W: usize>(scalar: &Scalar) {
        let digits = scalar.naf::<W>();

        // Horner's rule from the most significant digit, so `previous` is the next nonzero digit
        // above `position`.
        let mut reconstructed = Scalar::ZERO;
        let mut previous = None;
        for (position, &digit) in digits.iter().enumerate().rev() {
            reconstructed = reconstructed.add_mod_l(&reconstructed);
            if digit == 0 {
                continue;
            }
            assert_ne!(digit % 2, 0, "W={W} position={position}");
            assert!(
                digit.unsigned_abs() < 1 << (W - 1),
                "W={W} position={position}"
            );
            if let Some(previous) = previous {
                assert!(previous - position >= W, "W={W} position={position}");
            }
            previous = Some(position);

            // Add the signed digit as its residue mod `L`.
            let magnitude = Scalar::from_u128(u128::from(digit.unsigned_abs()));
            let term = if digit < 0 {
                magnitude.neg_mod_l()
            } else {
                magnitude
            };
            reconstructed = reconstructed.add_mod_l(&term);
        }
        assert_eq!(reconstructed.0, scalar.0, "W={W}");
    }

    #[test]
    fn naf_reconstructs_value_with_sparse_odd_digits() {
        Builder::default()
            .with_seed(0)
            .with_search_limit(64)
            .test(|u| {
                let scalar: Scalar = u.arbitrary()?;
                check_naf::<2>(&scalar);
                check_naf::<5>(&scalar);
                check_naf::<8>(&scalar);
                Ok(())
            });
    }

    /// Checks [`Scalar::half_size`]'s bounds and congruence, including inputs that skip the
    /// reduction loop and inputs whose first quotient spans over 100 bits.
    #[test]
    fn half_size_returns_short_multiple() {
        let check = |scalar: Scalar| {
            let (u, v) = scalar.half_size();
            assert_ne!(u, 0);
            assert!(u.unsigned_abs() < 1 << 127);
            assert!(v < 1 << 126);

            // Map the signed `u` to its residue mod `L` to check `u*scalar = v (mod L)`.
            let magnitude = Scalar::from_u128(u.unsigned_abs());
            let u = if u < 0 {
                magnitude.neg_mod_l()
            } else {
                magnitude
            };
            assert_eq!(u.mul_mod_l(&scalar).0, Scalar::from_u128(v).0);
        };

        // Values below `2^126` skip the reduction loop. From `2^126` through `2^128`, the first
        // quotient, `L` divided by the input, has more than 100 bits. `2^251` and `L - 1` reach
        // the top of the range.
        for value in [0, 1, 3, (1 << 126) - 1, 1 << 126, 1 << 127, u128::MAX] {
            check(Scalar::from_u128(value));
        }
        check(Scalar([0, 0, 1, 0]));
        check(Scalar([0, 0, 0, 1 << 59]));
        check(Scalar(limbs_sub(&L, &[1, 0, 0, 0])));
        Builder::default()
            .with_seed(0)
            .with_search_limit(1024)
            .test(|u| {
                check(u.arbitrary()?);
                Ok(())
            });
    }

    #[test]
    fn canonical_decoding_at_order_boundary() {
        for (limbs, valid) in [
            ([0; 4], true),
            (limbs_sub(&L, &[1, 0, 0, 0]), true),
            (L, false),
            (super::limbs_add(&L, &[1, 0, 0, 0]), false),
            ([u64::MAX; 4], false),
        ] {
            let bytes = Scalar(limbs).to_bytes();
            let decoded = Scalar::from_canonical_bytes(&bytes);
            assert_eq!(decoded.is_some(), valid);
            if let Some(decoded) = decoded {
                assert_eq!(decoded.to_bytes(), bytes);
            }
        }
    }

    #[test]
    fn windows_match_bits_across_limb_boundaries() {
        Builder::default()
            .with_seed(0)
            .with_search_limit(64)
            .test(|u| {
                let scalar: Scalar = u.arbitrary()?;
                let bytes = scalar.to_bytes();
                for width in 6..=10 {
                    for index in 0..=256usize.div_ceil(width as usize) {
                        let mut expected = 0;
                        for offset in 0..width as usize {
                            let bit = index * width as usize + offset;
                            if bit < 256 {
                                expected |= u32::from((bytes[bit / 8] >> (bit % 8)) & 1) << offset;
                            }
                        }
                        assert_eq!(scalar.window(index, width), expected);
                    }
                }
                Ok(())
            });
    }

    #[test]
    fn zeroize_clears_scalar() {
        let mut scalar = Scalar([u64::MAX; 4]);
        scalar.zeroize();
        assert_eq!(scalar.0, [0; 4]);
    }
}