cmpa 0.1.6

// SPDX-License-Identifier: Apache-2.0
// Copyright 2023, 2026 SUSE LLC
// Author: Nicolai Stange <nstange@suse.de>

//! Implementation of Montgomery multiplication and related primitives.
//!
//! Montgomery multiplication, as introduced in Peter L. Montgomery, "Modular
//! Multiplication Without Trial Division", Mathematics of Computation, vol. 44,
//! no. 170, pp. 519-521, 1985, computes modular products without any trial
//! divisions, for factors given in [Montgomery form](#montgomery-form). This
//! makes it the building block of choice for algorithms involving many modular
//! multiplications with the same modulus, like modular exponentiation.
//!
//! # Montgomery form
//!
//! For a given odd modulus *N*, the Montgomery form of an integer *a* is
//! *a · R mod N*, where *R* is the [Montgomery
//! radix](ct_montgomery_radix_shift_len), a power of two greater than *N*. For
//! the purpose of this implementation, *R* is set to
//! *2<sup>8 · `ct_montgomery_radix_shift_len(N.len())`</sup>*, where
//! `ct_montgomery_radix_shift_len()` rounds *N*'s [length in units of
//! bytes](MpUIntCommon::len) up to the next multiple of the [`LimbType`] size.
//! That is, *R* is two to the power of the number of bits in the [`LimbType`]
//! limbs spanned by *N*.
//!
//! The following primitives are provided for converting to and from the
//! Montgomery form:
//! - [`ct_to_montgomery_form_direct_mp()`] converts an integer to Montgomery
//!   form directly, by means of a division.
//! - [`ct_to_montgomery_form_mp()`] converts an integer to Montgomery form by a
//!   Montgomery multiplication with the precomputed *R<sup>2</sup> mod N*, c.f.
//!   [`ct_montgomery_radix2_mod_n_mp()`]. This is more efficient for converting
//!   multiple integers for the same modulus.
//! - [`ct_montgomery_redc_mp()`] converts an integer in Montgomery form back,
//!   by means of a Montgomery reduction.

use super::add::ct_sub_cond_mp_mp;
use super::cmp::{ct_geq_mp_mp, ct_lt_mp_mp};
use super::div::{
    CtModLshiftedMpMpError, CtModPow2MpError, CtMpDivisor, CtMpDivisorError, ct_mod_lshifted_mp_mp,
    ct_mod_pow2_mp,
};
use super::limb::{
    LIMB_BITS, LimbChoice, LimbType, ct_add_l_l, ct_eq_l_l, ct_inv_mod_l, ct_lsb_mask_l,
    ct_mul_add_l_l_l_c,
};
use super::limbs::{
    MpMutNativeEndianUIntLimbsSlice, MpMutUInt, MpMutUIntSlice, MpUIntCommon,
    ct_mp_limbs_align_len, ct_mp_nlimbs,
};

/// Determine the Montgomery radix's exponent in units of bytes.
///
/// For a modulus *N*, the Montgomery radix *R*, c.f. the section on the
/// [Montgomery form](self#montgomery-form), is chosen as
/// *2<sup>8 · `ct_montgomery_radix_shift_len(N.len())`</sup>* in this
/// implementation. That is, return the modulus' [length in units of
/// bytes](MpUIntCommon::len), rounded up to the next multiple of the
/// [`LimbType`] size, i.e. the length of the [`LimbType`] limbs spanned by the
/// modulus.
///
/// Runs in constant time.
///
/// # Arguments
///
/// * `n_len` - The modulus' [length in units of bytes](MpUIntCommon::len).
pub fn ct_montgomery_radix_shift_len(n_len: usize) -> usize {
    ct_mp_limbs_align_len(n_len)
}

/// Determine the Montgomery radix's exponent in units of limbs.
///
/// For a modulus *N*, the [Montgomery radix](ct_montgomery_radix_shift_len) is
/// *R = 2<sup>[`LimbType::BITS`] · `ct_montgomery_radix_shift_mp_nlimbs(N.len())`</sup>*.
/// That is, return the number of [`LimbType`] limbs spanned by the modulus. See
/// [`ct_montgomery_radix_shift_len()`] for the exponent in units of bytes.
///
/// Runs in constant time.
///
/// # Arguments
///
/// * `n_len` - The modulus' [length in units of bytes](MpUIntCommon::len).
fn ct_montgomery_radix_shift_mp_nlimbs(n_len: usize) -> usize {
    ct_mp_nlimbs(n_len)
}

/// Error returned by [`ct_montgomery_neg_n0_inv_mod_l_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtMontgomeryNegN0InvModLMpError {
    /// The modulus `n` is empty or even.
    InvalidModulus,
}

/// Compute the negated inverse of a modulus' least significant limb.
///
/// Returns *-`n`<sub>0</sub><sup>-1</sup> mod 2<sup>[`LimbType::BITS`]</sup>*,
/// where *`n`<sub>0</sub>* denotes `n`'s least significant limb. That is, the
/// returned value's product with *`n`<sub>0</sub>* equals *-1* modulo the
/// limb base. This value is needed by the Montgomery reduction and
/// multiplication, c.f. [`ct_montgomery_redc_mp()`],
/// [`ct_montgomery_mul_mod_mp_mp()`] and [`ct_to_montgomery_form_mp()`], and
/// is to be computed once per modulus.
///
/// Runs in constant time.
///
/// # Arguments
///
/// * `n` - The modulus. Must be odd.
///
/// # Errors
///
/// Returns [`CtMontgomeryNegN0InvModLMpError::InvalidModulus`] if `n` is empty
/// or even.
pub fn ct_montgomery_neg_n0_inv_mod_l_mp<NT: MpUIntCommon>(
    n: &NT,
) -> Result<LimbType, CtMontgomeryNegN0InvModLMpError> {
    if n.is_empty() {
        return Err(CtMontgomeryNegN0InvModLMpError::InvalidModulus);
    }
    let n0 = n.load_l(0);
    if ct_eq_l_l(n0 & 1, 0).unwrap() != 0 {
        return Err(CtMontgomeryNegN0InvModLMpError::InvalidModulus);
    }
    let n0_inv_mod_l = ct_inv_mod_l(n0);
    Ok((!n0_inv_mod_l).wrapping_add(1))
}

#[test]
fn test_ct_montgomery_neg_n0_inv_mod_l_mp() {
    use super::limb::LIMB_BYTES;
    use super::limbs::{MpMutBigEndianUIntByteSlice, MpMutUIntSlicePriv as _};

    for n0 in 0 as LimbType..128 {
        let n0 = 2 * n0 + 1;
        for j in 0..2048 {
            const MERSENNE_PRIME_13: LimbType = 8191 as LimbType;
            let v = MERSENNE_PRIME_13.wrapping_mul((511 as LimbType).wrapping_mul(j));
            let v = v << 8;
            let n0 = n0.wrapping_add(v);

            let mut n0_buf: [u8; LIMB_BYTES] = [0; LIMB_BYTES];
            let mut n = MpMutBigEndianUIntByteSlice::from_slice(n0_buf.as_mut_slice()).unwrap();
            n.store_l(0, n0);
            let neg_n0_inv = ct_montgomery_neg_n0_inv_mod_l_mp(&n).unwrap();
            assert_eq!(n0.wrapping_mul(neg_n0_inv), !0);
        }
    }
}

/// Incremental, limb-wise Montgomery reduction step.
///
/// Carries out a single step of Montgomery reduction by
/// *2<sup>e</sup>*, *e* ≤ [`LimbType::BITS`], on a multiprecision integer *t*
/// modulo an odd *n*: compute *(t + m · n) / 2<sup>e</sup>*, with
/// *m = t · (-n<sup>-1</sup>) mod 2<sup>e</sup>* chosen such that the division
/// is exact. The result is congruent to *t · 2<sup>-e</sup>* modulo *n*. A full
/// Montgomery reduction, c.f. [`ct_montgomery_redc_mp()`], consists of as many
/// such steps with *e* = [`LimbType::BITS`] as there are limbs in the
/// modulus.
///
/// The limbs of *t* and *n* are to be fed in from least to most significant,
/// one at a time, and the result limbs are returned likewise. This allows for
/// interleaving the reduction with other computations producing *t*'s limbs,
/// like the multiplication in [`ct_montgomery_mul_mod_mp_mp()`]:
/// - [`start()`](Self::start) consumes *t*'s and *n*'s least significant limbs.
/// - Each subsequent [`update()`](Self::update) consumes the respective next
///   limbs of *t* and *n* and returns the result's next limb.
/// - [`finish()`](Self::finish) consumes the limb above the last one fed to
///   [`update()`](Self::update), if any, and returns the result's most
///   significant limb along with a carry.
///
/// Each step runs in constant time.
pub(super) struct CtMontgomeryRedcKernel {
    /// The right shift distance for moving a reduced limb's bits above
    /// position *e* to the bottom of a result limb, i.e. *e* modulo
    /// [`LimbType::BITS`].
    redc_pow2_rshift: u32,
    /// Mask of the *e* least significant bits.
    redc_pow2_mask: LimbType,
    /// The multiplier
    /// *m = t<sub>0</sub> · (-n<sub>0</sub><sup>-1</sup>) mod 2<sup>e</sup>*,
    /// where *t<sub>0</sub>* and *n<sub>0</sub>* denote *t*'s and *n*'s least
    /// significant limbs.
    m: LimbType,
    /// The bits of the most recently reduced limb above position *e*, already
    /// shifted to the bottom. They make up the lower part of the result limb
    /// returned next.
    last_redced_val: LimbType,
    /// The carry of the running computation of *t + m · n*.
    carry: LimbType,
}

impl CtMontgomeryRedcKernel {
    /// Extract a reduced limb's bits above position *e*, shifted to the bottom.
    ///
    /// They make up the lower part of a result limb.
    ///
    /// # Arguments
    ///
    /// * `redc_pow2_rshift` - c.f.
    ///   [`redc_pow2_rshift`](Self::redc_pow2_rshift).
    /// * `redc_pow2_mask` - c.f. [`redc_pow2_mask`](Self::redc_pow2_mask).
    /// * `val` - The reduced limb.
    const fn redc_rshift_lo(
        redc_pow2_rshift: u32,
        redc_pow2_mask: LimbType,
        val: LimbType,
    ) -> LimbType {
        // There are two possible cases where redc_pow2_rshift == 0:
        // a.) redc_pow2_exp == 0. In this case !redc_pow2_mask == !0.
        // b.) redc_pow2_exp == LIMB_BITS. In this case !redc_pow2_mask == !!0 == 0.
        (val & !redc_pow2_mask) >> redc_pow2_rshift
    }

    /// Extract a reduced limb's *e* least significant bits, shifted to the top.
    ///
    /// They make up the upper part of the preceding result limb.
    ///
    /// # Arguments
    ///
    /// * `redc_pow2_rshift` - c.f.
    ///   [`redc_pow2_rshift`](Self::redc_pow2_rshift).
    /// * `redc_pow2_mask` - c.f. [`redc_pow2_mask`](Self::redc_pow2_mask).
    /// * `val` - The reduced limb.
    const fn redc_lshift_hi(
        redc_pow2_rshift: u32,
        redc_pow2_mask: LimbType,
        val: LimbType,
    ) -> LimbType {
        // There are two possible cases where redc_pow2_rshift == 0:
        // a.) redc_pow2_exp == 0. In this case redc_pow2_mask == 0 as well.
        // b.) redc_pow2_exp == LIMB_BITS. In this case redc_pow2_mask == !0.
        (val & redc_pow2_mask) << ((LIMB_BITS - redc_pow2_rshift) % LIMB_BITS)
    }

    /// Start a reduction step.
    ///
    /// Compute the multiplier *m* and consume *t*'s and *n*'s least
    /// significant limbs.
    ///
    /// # Arguments
    ///
    /// * `redc_pow2_exp` - The exponent *e* of the power of two to reduce by.
    ///   Must be at most [`LimbType::BITS`].
    /// * `t0_val` - *t*'s least significant limb.
    /// * `n0_val` - *n*'s least significant limb. Must be odd.
    /// * `neg_n0_inv_mod_l` - The negated inverse of `n0_val` modulo
    ///   *2<sup>[`LimbType::BITS`]</sup>*, as computed by
    ///   [`ct_montgomery_neg_n0_inv_mod_l_mp()`].
    pub fn start(
        redc_pow2_exp: u32,
        t0_val: LimbType,
        n0_val: LimbType,
        neg_n0_inv_mod_l: LimbType,
    ) -> Self {
        debug_assert!(redc_pow2_exp <= LIMB_BITS);

        // Calculate the shift distance for right shifting a redced limb into its final
        // position. Be careful to keep it < LIMB_BITS for not running into
        // undefined behaviour with the shift.  See the comments in
        // redc_rshift_lo()/redc_rshift_hi() for the mask<->rshift interaction.
        let redc_pow2_rshift = redc_pow2_exp % LIMB_BITS;
        let redc_pow2_mask = ct_lsb_mask_l(redc_pow2_exp);

        // For any i >= j, if n' == -n^{-1} mod 2^i, then n' mod 2^j == -n^{-1} mod 2^j.
        let neg_n0_inv_mod_l = neg_n0_inv_mod_l & redc_pow2_mask;
        debug_assert_eq!(
            neg_n0_inv_mod_l.wrapping_mul(n0_val) & redc_pow2_mask,
            redc_pow2_mask
        );

        let m = t0_val.wrapping_mul(neg_n0_inv_mod_l) & redc_pow2_mask;

        let (carry, redced_t0_val) = ct_mul_add_l_l_l_c(t0_val, m, n0_val, 0);
        debug_assert_eq!(redced_t0_val & redc_pow2_mask, 0);
        // If redc_pow2_exp < LIMB_BITS, the upper bits of the reduced zeroth limb
        // will become the lower bits of the resulting zeroth limb.
        let last_redced_val = Self::redc_rshift_lo(redc_pow2_rshift, redc_pow2_mask, redced_t0_val);

        Self {
            redc_pow2_rshift,
            redc_pow2_mask,
            m,
            last_redced_val,
            carry,
        }
    }

    /// Feed the next pair of limbs and return the next result limb.
    ///
    /// The `j`'th invocation, counting from zero, consumes *t*'s and *n*'s
    /// limbs at index `j + 1` and returns the result's limb at index `j`.
    ///
    /// # Arguments
    ///
    /// * `t_val` - *t*'s next limb.
    /// * `n_val` - *n*'s next limb, or zero if *n* has no more limbs.
    pub fn update(&mut self, t_val: LimbType, n_val: LimbType) -> LimbType {
        let redced_t_val;
        (self.carry, redced_t_val) = ct_mul_add_l_l_l_c(t_val, self.m, n_val, self.carry);

        // If redc_pow2_exp < LIMB_BITS, the lower bits of the reduced current limb
        // correspond to the upper bits of the returned result limb.
        let result_val = self.last_redced_val
            | Self::redc_lshift_hi(self.redc_pow2_rshift, self.redc_pow2_mask, redced_t_val);
        // If redc_pow2_exp < LIMB_BITS, the upper bits of the reduced current limb
        // will become the lower bits of the subsequently returned result limb.
        self.last_redced_val =
            Self::redc_rshift_lo(self.redc_pow2_rshift, self.redc_pow2_mask, redced_t_val);

        result_val
    }

    /// Finish the reduction step.
    ///
    /// Consume *t*'s last limb, i.e. the one above the last one fed to
    /// [`update()`](Self::update), and return a pair of the carry beyond the
    /// result's most significant limb and that limb itself.
    ///
    /// # Arguments
    ///
    /// * `t_val` - *t*'s last limb, which, for a full Montgomery reduction,
    ///   is usually the carry returned from the preceding reduction step.
    pub fn finish(self, t_val: LimbType) -> (LimbType, LimbType) {
        let (carry_hi, redced_t_val) = ct_add_l_l(t_val, self.carry);
        let carry_hi = Self::redc_lshift_hi(self.redc_pow2_rshift, self.redc_pow2_mask, carry_hi);
        let carry_lo =
            Self::redc_rshift_lo(self.redc_pow2_rshift, self.redc_pow2_mask, redced_t_val);
        (
            carry_hi | carry_lo,
            self.last_redced_val
                | Self::redc_lshift_hi(self.redc_pow2_rshift, self.redc_pow2_mask, redced_t_val),
        )
    }
}

/// Error returned by [`ct_montgomery_redc_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtMontgomeryRedcMpError {
    /// The modulus `n` is empty or even.
    InvalidModulus,
    /// The modulus' [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) `t`'s.
    InsufficientResultSpace,
    /// `t`'s [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) twice the modulus'.
    InputValueOutOfRange,
}

/// Montgomery reduction, converting from Montgomery form.
///
/// The multiprecision integers may be of different representations.
///
/// Replace `t` by *`t` · R<sup>-1</sup> mod `n`*, where *R* is the [Montgomery
/// radix](ct_montgomery_radix_shift_len), c.f. the section on the [Montgomery
/// form](self#montgomery-form). The result is less than `n`. In particular, for
/// `t` in [Montgomery form](self#montgomery-form), i.e. *`t` = a · R mod `n`*
/// for some *a*, this converts back to *a mod `n`*. Likewise, applied to the
/// product of two integers in [Montgomery form](self#montgomery-form), it
/// yields their product in [Montgomery form](self#montgomery-form).
///
/// Execution time depends only on the multiprecision integers' [lengths in
/// units of bytes](MpUIntCommon::len) and is otherwise constant as far as
/// branching is concerned. The computation cost is O(m · n), where m is `t`'s
/// and n is `n`'s [length](MpUIntCommon::len). See
/// [`ct_to_montgomery_form_mp()`] and [`ct_to_montgomery_form_direct_mp()`] for
/// converting to [Montgomery form](self#montgomery-form).
///
/// # Arguments
///
/// * `t` - The multiprecision integer to reduce. Its value must be less than
///   *`n` · R*, which holds e.g. for the product of two integers less than `n`.
///   Its [length in units of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) twice `n`'s. It will be
///   overwritten by the result.
/// * `n` - The modulus. Must be odd. Its [length in units of
///   bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `t`'s.
/// * `neg_n0_inv_mod_l` - The negated inverse of `n`'s least significant limb
///   modulo *2<sup>[`LimbType::BITS`]</sup>*, as computed by
///   [`ct_montgomery_neg_n0_inv_mod_l_mp()`].
///
/// # Errors
///
/// - Returns [`CtMontgomeryRedcMpError::InvalidModulus`] if `n` is empty or
///   even.
/// - Returns [`CtMontgomeryRedcMpError::InsufficientResultSpace`] if `n`'s
///   length doesn't fit into `t`'s.
/// - Returns [`CtMontgomeryRedcMpError::InputValueOutOfRange`] if `t`'s length
///   doesn't fit into twice `n`'s.
///
/// Nothing gets modified in any of these cases.
pub fn ct_montgomery_redc_mp<TT: MpMutUInt, NT: MpUIntCommon>(
    t: &mut TT,
    n: &NT,
    neg_n0_inv_mod_l: LimbType,
) -> Result<(), CtMontgomeryRedcMpError> {
    if n.is_empty() {
        return Err(CtMontgomeryRedcMpError::InvalidModulus);
    }
    let n0_val = n.load_l(0);
    if ct_eq_l_l(n0_val & 1, 0).unwrap() != 0 {
        return Err(CtMontgomeryRedcMpError::InvalidModulus);
    }
    if !n.len_is_compatible_with(t.len()) {
        return Err(CtMontgomeryRedcMpError::InsufficientResultSpace);
    }
    if !t.len_is_compatible_with(2 * n.len()) {
        return Err(CtMontgomeryRedcMpError::InputValueOutOfRange);
    }
    let t_nlimbs = t.nlimbs();
    let n_nlimbs = n.nlimbs();
    debug_assert!(n0_val.wrapping_mul(neg_n0_inv_mod_l) == !0);

    let mut reduced_t_carry = 0;
    // t's high limb might be a partial one, do not update directly in the course of
    // reducing in order to avoid overflowing it. Use a shadow instead.
    let mut t_high_shadow = t.load_l(t_nlimbs - 1);
    for _i in 0..ct_montgomery_radix_shift_mp_nlimbs(n.len()) {
        let mut redc_kernel =
            CtMontgomeryRedcKernel::start(LIMB_BITS, t.load_l(0), n0_val, neg_n0_inv_mod_l);
        let mut j = 0;
        while j + 2 < n_nlimbs {
            debug_assert!(j < t_nlimbs - 1);
            t.store_l_full(
                j,
                redc_kernel.update(t.load_l_full(j + 1), n.load_l_full(j + 1)),
            );
            j += 1;
        }

        if n_nlimbs != 1 {
            debug_assert_eq!(j + 2, n_nlimbs);
            // Do not read the potentially partial, stale high limb directly from t, use the
            // t_high_shadow shadow instead.
            let t_val = if j + 2 != t_nlimbs {
                t.load_l_full(j + 1)
            } else {
                t_high_shadow
            };
            t.store_l_full(j, redc_kernel.update(t_val, n.load_l(j + 1)));
            j += 1;
        }

        while j + 2 < t_nlimbs {
            t.store_l_full(j, redc_kernel.update(t.load_l_full(j + 1), 0));
            j += 1;
        }
        if j + 1 == t_nlimbs - 1 {
            t.store_l_full(j, redc_kernel.update(t_high_shadow, 0));
            j += 1;
        }
        debug_assert_eq!(j, t_nlimbs - 1);

        (reduced_t_carry, t_high_shadow) = redc_kernel.finish(reduced_t_carry);
        debug_assert!(reduced_t_carry <= 1);
    }

    // Now apply the high limb shadow back.
    let t_high_shadow_mask = t.partial_high_mask();
    let t_high_shadow_shift = t.partial_high_shift();
    assert!(t_high_shadow_shift == 0 || reduced_t_carry == 0);
    reduced_t_carry |= (t_high_shadow & !t_high_shadow_mask) >> t_high_shadow_shift;
    t_high_shadow &= t_high_shadow_mask;
    t.store_l(t_nlimbs - 1, t_high_shadow);

    ct_sub_cond_mp_mp(t, n, LimbChoice::from(reduced_t_carry) | ct_geq_mp_mp(t, n));
    debug_assert!(ct_geq_mp_mp(t, n).unwrap() == 0);
    Ok(())
}

#[cfg(test)]
fn test_ct_montgomery_redc_mp<TT: MpMutUIntSlice, NT: MpMutUIntSlice>() {
    use super::div::ct_mod_mp_mp;
    use super::limb::LIMB_BYTES;

    for i in 0..64 {
        const MERSENNE_PRIME_13: LimbType = 8191 as LimbType;
        let n_high = MERSENNE_PRIME_13.wrapping_mul((16385 as LimbType).wrapping_mul(i));
        for j in 0..64 {
            const MERSENNE_PRIME_17: LimbType = 131071 as LimbType;
            let n_low = MERSENNE_PRIME_17.wrapping_mul((1023 as LimbType).wrapping_mul(j));
            // Force n_low odd.
            let n_low = n_low | 1;
            let mut n = tst_mk_mp_backing_vec!(NT, 2 * LIMB_BYTES);
            let mut n = NT::from_slice(n.as_mut_slice()).unwrap();
            n.store_l(0, n_low);
            n.store_l(1, n_high);
            let neg_n0_inv = ct_montgomery_neg_n0_inv_mod_l_mp(&n).unwrap();

            for k in 0..8 {
                let t_high = MERSENNE_PRIME_17.wrapping_mul((8191 as LimbType).wrapping_mul(k));
                for l in 0..8 {
                    let t_low =
                        MERSENNE_PRIME_13.wrapping_mul((131087 as LimbType).wrapping_mul(l));

                    let mut t = tst_mk_mp_backing_vec!(TT, 2 * LIMB_BYTES);
                    let mut t = TT::from_slice(t.as_mut_slice()).unwrap();
                    t.store_l(0, t_low);
                    t.store_l(1, t_high);

                    // All montgomery operations are defined mod n, compute t mod n
                    ct_mod_mp_mp(None, &mut t, &CtMpDivisor::new(&n, None).unwrap());
                    let t_low = t.load_l(0);
                    let t_high = t.load_l(1);

                    // To Montgomery form: t * R mod n
                    ct_to_montgomery_form_direct_mp(&mut t, &n).unwrap();

                    // And back to normal: (t * R mod n) / R mod n
                    ct_montgomery_redc_mp(&mut t, &n, neg_n0_inv).unwrap();
                    assert_eq!(t.load_l(0), t_low);
                    assert_eq!(t.load_l(1), t_high);
                }
            }
        }
    }
}

#[test]
fn test_ct_montgomery_redc_be_be() {
    use super::limbs::MpMutBigEndianUIntByteSlice;
    test_ct_montgomery_redc_mp::<MpMutBigEndianUIntByteSlice, MpMutBigEndianUIntByteSlice>()
}

#[test]
fn test_ct_montgomery_redc_le_le() {
    use super::limbs::MpMutLittleEndianUIntByteSlice;
    test_ct_montgomery_redc_mp::<MpMutLittleEndianUIntByteSlice, MpMutLittleEndianUIntByteSlice>()
}

#[test]
fn test_ct_montgomery_redc_ne_ne() {
    use super::limbs::MpMutNativeEndianUIntLimbsSlice;
    test_ct_montgomery_redc_mp::<MpMutNativeEndianUIntLimbsSlice, MpMutNativeEndianUIntLimbsSlice>()
}

/// Error returned by [`ct_montgomery_mul_mod_mp_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtMontgomeryMulModMpMpError {
    /// The modulus `n` is empty or even.
    InvalidModulus,
    /// The modulus' [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) `result`'s.
    InsufficientResultSpace,
    /// `op0`'s or `op1`'s [length in units of bytes](MpUIntCommon::len) doesn't
    /// [fit into](MpUIntCommon::len_is_compatible_with) the modulus'.
    InconsistentOperandLengths,
}

/// Montgomery multiplication.
///
/// The multiprecision integers may be of different representations.
///
/// Compute *`op0` · `op1` · R<sup>-1</sup> mod `n`*, where *R* is the
/// [Montgomery radix](ct_montgomery_radix_shift_len), and store it in `result`,
/// with any of `result`'s bytes above `n`'s length cleared. The result is less
/// than `n`. For factors in [Montgomery form](self#montgomery-form), i.e.
/// *`op0` = a · R mod `n`* and *`op1` = b · R mod `n`*, the result is the
/// product's [Montgomery form](self#montgomery-form), *a · b · R mod `n`*.
///
/// The multiplication and the Montgomery reduction are fused, following the
/// "Finely Integrated Operand Scanning (FIOS)" method from Ç. K. Koç, T. Acar
/// and B. S. Kaliski, "Analyzing and Comparing Montgomery Multiplication
/// Algorithms", IEEE Micro, 16(3):26-33, June 1996.
///
/// Execution time depends only on the multiprecision integers' [lengths in
/// units of bytes](MpUIntCommon::len) and is otherwise constant as far as
/// branching is concerned. The computation cost is O(n<sup>2</sup>), where n
/// is `n`'s [length](MpUIntCommon::len). See [`ct_to_montgomery_form_mp()`]
/// and [`ct_to_montgomery_form_direct_mp()`] for converting the factors to,
/// and [`ct_montgomery_redc_mp()`] for converting the result back from
/// [Montgomery form](self#montgomery-form).
///
/// # Arguments
///
/// * `result` - The output for the result. Its [length in units of
///   bytes](MpUIntCommon::len) must be such that `n`'s [fits
///   into](MpUIntCommon::len_is_compatible_with) it.
/// * `op0` - The first factor, usually in [Montgomery
///   form](self#montgomery-form). Must be less than `n`. Its [length in units
///   of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `n`'s.
/// * `op1` - The second factor, usually in [Montgomery
///   form](self#montgomery-form). Must be less than `n`. Its [length in units
///   of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `n`'s.
/// * `n` - The modulus. Must be odd.
/// * `neg_n0_inv_mod_l` - The negated inverse of `n`'s least significant limb
///   modulo *2<sup>[`LimbType::BITS`]</sup>*, as computed by
///   [`ct_montgomery_neg_n0_inv_mod_l_mp()`].
///
/// # Errors
///
/// - Returns [`CtMontgomeryMulModMpMpError::InvalidModulus`] if `n` is empty
///   or even.
/// - Returns [`CtMontgomeryMulModMpMpError::InsufficientResultSpace`] if `n`'s
///   length doesn't fit into `result`'s.
/// - Returns [`CtMontgomeryMulModMpMpError::InconsistentOperandLengths`] if
///   `op0`'s or `op1`'s length doesn't fit into `n`'s.
///
/// Nothing gets modified in any of these cases.
pub fn ct_montgomery_mul_mod_mp_mp<
    RT: MpMutUIntSlice,
    T0: MpUIntCommon,
    T1: MpUIntCommon,
    NT: MpUIntCommon,
>(
    result: &mut RT,
    op0: &T0,
    op1: &T1,
    n: &NT,
    neg_n0_inv_mod_l: LimbType,
) -> Result<(), CtMontgomeryMulModMpMpError> {
    // This is an implementation of the "Finely Integrated Operand Scanning (FIOS)
    // Method" approach to fused multiplication and Montgomery reduction, as
    // described in "Analyzing and Comparing Montgomery Multiplication
    // Algorithm", IEEE Micro, 16(3):26-33, June 1996.
    if n.is_empty() {
        return Err(CtMontgomeryMulModMpMpError::InvalidModulus);
    }
    let n0_val = n.load_l(0);
    if ct_eq_l_l(n0_val & 1, 0).unwrap() != 0 {
        return Err(CtMontgomeryMulModMpMpError::InvalidModulus);
    }
    if !n.len_is_compatible_with(result.len()) {
        return Err(CtMontgomeryMulModMpMpError::InsufficientResultSpace);
    }
    debug_assert!(n.nlimbs() <= result.nlimbs());
    if !op0.len_is_compatible_with(n.len()) || !op1.len_is_compatible_with(n.len()) {
        return Err(CtMontgomeryMulModMpMpError::InconsistentOperandLengths);
    }
    debug_assert!(op0.nlimbs() <= n.nlimbs());
    debug_assert!(ct_lt_mp_mp(op0, n).unwrap() != 0);
    debug_assert!(op1.nlimbs() <= n.nlimbs());
    debug_assert!(ct_lt_mp_mp(op1, n).unwrap() != 0);

    let op0_nlimbs = op0.nlimbs();
    let op1_nlimbs = op1.nlimbs();
    let n_nlimbs = n.nlimbs();
    debug_assert!(n0_val.wrapping_mul(neg_n0_inv_mod_l) == !0);

    result.clear_bytes_above(0);
    let mut result = result.shrink_to(n.len());
    debug_assert_eq!(result.nlimbs(), n.nlimbs());
    let mut result_carry = 0;
    // result's high limb might be a partial one, do not update directly in the
    // course of reducing in order to avoid overflowing it. Use a shadow
    // instead.
    let mut result_high_shadow = 0;
    for i in 0..op0_nlimbs {
        debug_assert!(result_carry <= 1); // Loop invariant.
        let op0_val = op0.load_l(i);

        // Do not read the potentially partial, stale high limb directly from result,
        // use the result_high_shadow shadow instead.
        let result_val = if n_nlimbs != 1 {
            result.load_l_full(0)
        } else {
            result_high_shadow
        };
        let op1_val = op1.load_l(0);
        let (mut op0_op1_add_carry, result_val) =
            ct_mul_add_l_l_l_c(result_val, op0_val, op1_val, 0);

        let mut redc_kernel =
            CtMontgomeryRedcKernel::start(LIMB_BITS, result_val, n0_val, neg_n0_inv_mod_l);

        let mut j = 0;
        while j + 1 < op1_nlimbs {
            let op1_val = op1.load_l(j + 1);

            // Do not read the potentially partial, stale high limb directly from result,
            // use the result_high_shadow shadow instead.
            let mut result_val = if j + 1 != n_nlimbs - 1 {
                result.load_l_full(j + 1)
            } else {
                result_high_shadow
            };

            (op0_op1_add_carry, result_val) =
                ct_mul_add_l_l_l_c(result_val, op0_val, op1_val, op0_op1_add_carry);

            let n_val = n.load_l(j + 1);
            let result_val = redc_kernel.update(result_val, n_val);
            result.store_l_full(j, result_val);
            j += 1;
        }
        debug_assert_eq!(j + 1, op1_nlimbs);

        // If op1_nlimbs < n_nlimbs, handle the rest by propagating the multiplication
        // carry and continue redcing.
        while j + 1 < n_nlimbs {
            // Do not read the potentially partial, stale high limb directly from result,
            // use the result_high_shadow shadow instead.
            let mut result_val = if j + 1 != n_nlimbs - 1 {
                result.load_l_full(j + 1)
            } else {
                result_high_shadow
            };

            (op0_op1_add_carry, result_val) = ct_add_l_l(result_val, op0_op1_add_carry);

            let n_val = n.load_l(j + 1);
            let result_val = redc_kernel.update(result_val, n_val);
            result.store_l_full(j, result_val);
            j += 1;
        }
        debug_assert_eq!(j + 1, n_nlimbs);

        let mut result_val;
        debug_assert!(result_carry <= 1);
        (result_carry, result_val) = ct_add_l_l(result_carry, op0_op1_add_carry);
        debug_assert!(result_carry <= 1);
        debug_assert!(result_carry == 0 || result_val == 0);

        (result_carry, result_val) = redc_kernel.finish(result_val);
        debug_assert!(result_carry <= 1);
        result_high_shadow = result_val;
    }

    // If op0_nlimbs < the montgomery radix shift distance, handle the rest by
    // REDCing it.
    for _i in op0_nlimbs..ct_montgomery_radix_shift_mp_nlimbs(n.len()) {
        // Do not read the potentially partial, stale high limb directly from result,
        // use the result_high_shadow shadow instead.
        let result_val = if n_nlimbs != 1 {
            result.load_l_full(0)
        } else {
            result_high_shadow
        };

        let mut redc_kernel =
            CtMontgomeryRedcKernel::start(LIMB_BITS, result_val, n0_val, neg_n0_inv_mod_l);

        let mut j = 0;
        while j + 1 < n_nlimbs {
            // Do not read the potentially partial, stale high limb directly from result,
            // use the result_high_shadow shadow instead.
            let result_val = if j + 1 != n_nlimbs - 1 {
                result.load_l_full(j + 1)
            } else {
                result_high_shadow
            };

            let n_val = n.load_l(j + 1);
            let result_val = redc_kernel.update(result_val, n_val);
            result.store_l_full(j, result_val);
            j += 1;
        }
        debug_assert_eq!(j + 1, n_nlimbs);

        (result_carry, result_high_shadow) = redc_kernel.finish(result_carry);
        debug_assert!(result_carry <= 1);
    }

    // Now apply the high limb shadow back.
    debug_assert!(result.nlimbs() == n.nlimbs());
    let result_high_shadow_mask = n.partial_high_mask();
    let result_high_shadow_shift = n.partial_high_shift();
    assert!(result_high_shadow_shift == 0 || result_carry == 0);
    result_carry |= (result_high_shadow & !result_high_shadow_mask) >> result_high_shadow_shift;
    result_high_shadow &= result_high_shadow_mask;
    result.store_l(n_nlimbs - 1, result_high_shadow);

    let result_geq_n = LimbChoice::from(result_carry) | ct_geq_mp_mp(&result, n);
    ct_sub_cond_mp_mp(&mut result, n, result_geq_n);
    debug_assert!(ct_geq_mp_mp(&result, n).unwrap() == 0);
    Ok(())
}

#[cfg(test)]
fn test_ct_montgomery_mul_mod_mp_mp<
    RT: MpMutUIntSlice,
    T0: MpMutUIntSlice,
    T1: MpMutUIntSlice,
    NT: MpMutUIntSlice,
>() {
    use super::div::ct_mod_mp_mp;
    use super::limb::LIMB_BYTES;
    use super::mul::ct_mul_trunc_mp_mp;

    for i in 0..16 {
        const MERSENNE_PRIME_13: LimbType = 8191 as LimbType;
        let n_high = MERSENNE_PRIME_13.wrapping_mul((65543 as LimbType).wrapping_mul(i));
        for j in 0..16 {
            const MERSENNE_PRIME_17: LimbType = 131071 as LimbType;
            let n_low = MERSENNE_PRIME_17.wrapping_mul((4095 as LimbType).wrapping_mul(j));
            // Force n_low odd.
            let n_low = n_low | 1;
            let mut n_buf = tst_mk_mp_backing_vec!(NT, 2 * LIMB_BYTES);
            let mut n = NT::from_slice(n_buf.as_mut_slice()).unwrap();
            n.store_l(0, n_low);
            n.store_l(1, n_high);
            drop(n);
            let n_lengths = if !RT::SUPPORTS_UNALIGNED_BUFFER_LENGTHS
                || !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS
                || !T1::SUPPORTS_UNALIGNED_BUFFER_LENGTHS
                || !NT::SUPPORTS_UNALIGNED_BUFFER_LENGTHS
            {
                [LIMB_BYTES, 2 * LIMB_BYTES]
            } else {
                [2 * LIMB_BYTES - 1, 2 * LIMB_BYTES]
            };
            for n_len in n_lengths {
                let mut n_buf = n_buf.clone();
                let mut n = NT::from_slice(n_buf.as_mut_slice()).unwrap();
                n.clear_bytes_above(n_len);
                let n = n.shrink_to(n_len);
                let neg_n0_inv = ct_montgomery_neg_n0_inv_mod_l_mp(&n).unwrap();

                // r_mod_n = 2^(2 * LIMB_BITS) % n.
                let mut r_mod_n = tst_mk_mp_backing_vec!(RT, 3 * LIMB_BYTES);
                let mut r_mod_n = RT::from_slice(r_mod_n.as_mut_slice()).unwrap();
                r_mod_n.store_l_full(ct_montgomery_radix_shift_mp_nlimbs(n_len), 1);
                ct_mod_mp_mp(None, &mut r_mod_n, &CtMpDivisor::new(&n, None).unwrap());
                let r_mod_n = r_mod_n.shrink_to(n.len());

                for k in 0..4 {
                    let a_high =
                        MERSENNE_PRIME_17.wrapping_mul((16383 as LimbType).wrapping_mul(k));
                    for l in 0..4 {
                        let a_low =
                            MERSENNE_PRIME_13.wrapping_mul((262175 as LimbType).wrapping_mul(l));
                        let mut a_buf = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES);
                        let mut a = T0::from_slice(a_buf.as_mut_slice()).unwrap();
                        a.store_l(0, a_low);
                        a.store_l(1, a_high);
                        // All montgomery operations are defined mod n, compute a mod n
                        ct_mod_mp_mp(None, &mut a, &CtMpDivisor::new(&n, None).unwrap());
                        drop(a);
                        for s in 0..4 {
                            let b_high = MERSENNE_PRIME_13
                                .wrapping_mul((262175 as LimbType).wrapping_mul(s));
                            for t in 0..4 {
                                const MERSENNE_PRIME_19: LimbType = 524287 as LimbType;
                                let b_low = MERSENNE_PRIME_19
                                    .wrapping_mul((4095 as LimbType).wrapping_mul(t));
                                let mut b_buf = tst_mk_mp_backing_vec!(T1, 2 * LIMB_BYTES);
                                let mut b = T1::from_slice(b_buf.as_mut_slice()).unwrap();
                                b.store_l(0, b_low);
                                b.store_l(1, b_high);
                                // All montgomery operations are defined mod n, compute b mod n
                                ct_mod_mp_mp(None, &mut b, &CtMpDivisor::new(&n, None).unwrap());
                                drop(b);

                                for op_len in [0, LIMB_BYTES, n_len] {
                                    let mut a_buf = a_buf.clone();
                                    let mut a = T0::from_slice(a_buf.as_mut_slice()).unwrap();
                                    a.clear_bytes_above(op_len);
                                    let a = a.shrink_to(op_len);
                                    let mut b_buf = b_buf.clone();
                                    let mut b = T1::from_slice(b_buf.as_mut_slice()).unwrap();
                                    b.clear_bytes_above(op_len);
                                    let b = b.shrink_to(op_len);

                                    let mut _result = tst_mk_mp_backing_vec!(RT, 4 * LIMB_BYTES);
                                    let mut result =
                                        RT::from_slice(_result.as_mut_slice()).unwrap();
                                    let mut mg_mul_result = result.shrink_to(n_len);
                                    ct_montgomery_mul_mod_mp_mp(
                                        &mut mg_mul_result,
                                        &a,
                                        &b,
                                        &n,
                                        neg_n0_inv,
                                    )
                                    .unwrap();
                                    drop(mg_mul_result);

                                    // For testing against the expected result computed using the
                                    // "conventional" methods only, multiply by r_mod_n -- this
                                    // avoids having to multiply
                                    // the conventional product by r^-1 mod n, which is
                                    // not known without implementing Euklid's algorithm.
                                    ct_mul_trunc_mp_mp(&mut result, n.len(), &r_mod_n);
                                    ct_mod_mp_mp(
                                        None,
                                        &mut result,
                                        &CtMpDivisor::new(&n, None).unwrap(),
                                    );
                                    drop(result);

                                    let mut _expected = tst_mk_mp_backing_vec!(RT, 4 * LIMB_BYTES);
                                    let mut expected =
                                        RT::from_slice(_expected.as_mut_slice()).unwrap();
                                    expected.copy_from(&a);
                                    ct_mul_trunc_mp_mp(&mut expected, op_len, &b);
                                    ct_mod_mp_mp(
                                        None,
                                        &mut expected,
                                        &CtMpDivisor::new(&n, None).unwrap(),
                                    );
                                    drop(expected);

                                    assert_eq!(_result, _expected);
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

#[test]
fn test_ct_montgomery_mul_mod_be_be_be_be() {
    use super::limbs::MpMutBigEndianUIntByteSlice;
    test_ct_montgomery_mul_mod_mp_mp::<
        MpMutBigEndianUIntByteSlice,
        MpMutBigEndianUIntByteSlice,
        MpMutBigEndianUIntByteSlice,
        MpMutBigEndianUIntByteSlice,
    >()
}

#[test]
fn test_ct_montgomery_mul_mod_le_le_le_le() {
    use super::limbs::MpMutLittleEndianUIntByteSlice;
    test_ct_montgomery_mul_mod_mp_mp::<
        MpMutLittleEndianUIntByteSlice,
        MpMutLittleEndianUIntByteSlice,
        MpMutLittleEndianUIntByteSlice,
        MpMutLittleEndianUIntByteSlice,
    >()
}

#[test]
fn test_ct_montgomery_mul_mod_ne_ne_ne_ne() {
    use super::limbs::MpMutNativeEndianUIntLimbsSlice;
    test_ct_montgomery_mul_mod_mp_mp::<
        MpMutNativeEndianUIntLimbsSlice,
        MpMutNativeEndianUIntLimbsSlice,
        MpMutNativeEndianUIntLimbsSlice,
        MpMutNativeEndianUIntLimbsSlice,
    >()
}

/// Error returned by [`ct_to_montgomery_form_direct_mp()`] and
/// [`ct_montgomery_radix2_mod_n_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtMontgomeryTransformationError {
    /// The modulus `n` is even, including the case of an empty one.
    InvalidModulus,
    /// The modulus' [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) the result's.
    InsufficientResultSpace,
}

/// Convert a multiprecision integer to Montgomery form by division.
///
/// The multiprecision integers may be of different representations.
///
/// Replace `t` by its [Montgomery form](self#montgomery-form),
/// *`t` · R mod `n`*, where *R* is the [Montgomery
/// radix](ct_montgomery_radix_shift_len). The result is less than `n`. `t`'s
/// value isn't restricted, in particular it may exceed `n`.
///
/// The conversion is carried out directly, by means of a division. See
/// [`ct_to_montgomery_form_mp()`] for a variant based on Montgomery
/// multiplication with a precomputed *R<sup>2</sup> mod `n`*, which is more
/// efficient when converting multiple integers for the same modulus, and
/// [`ct_montgomery_redc_mp()`] for converting back.
///
/// Execution time depends only on the multiprecision integers' [lengths in
/// units of bytes](MpUIntCommon::len) as well as on `n`'s length without any
/// leading zero bytes, c.f. [`CtMpDivisor`], and is otherwise constant as far
/// as branching is concerned. The computation cost is O((m + n) · n), where m
/// is `t`'s and n is `n`'s [length](MpUIntCommon::len).
///
/// # Arguments
///
/// * `t` - The multiprecision integer to convert. It will be overwritten by the
///   result.
/// * `n` - The modulus. Must be odd. Its [length in units of
///   bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `t`'s.
///
/// # Errors
///
/// - Returns [`CtMontgomeryTransformationError::InvalidModulus`] if `n` is
///   even.
/// - Returns [`CtMontgomeryTransformationError::InsufficientResultSpace`] if
///   `n`'s length doesn't fit into `t`'s.
///
/// Nothing gets modified in either case.
pub fn ct_to_montgomery_form_direct_mp<TT: MpMutUInt, NT: MpUIntCommon>(
    t: &mut TT,
    n: &NT,
) -> Result<(), CtMontgomeryTransformationError> {
    if n.test_bit(0).unwrap() == 0 {
        return Err(CtMontgomeryTransformationError::InvalidModulus);
    }
    if !n.len_is_compatible_with(t.len()) {
        return Err(CtMontgomeryTransformationError::InsufficientResultSpace);
    }
    debug_assert!(t.nlimbs() >= n.nlimbs());
    let radix_shift_len = ct_montgomery_radix_shift_len(n.len());
    let n = CtMpDivisor::new(n, None).map_err(|e| match e {
        CtMpDivisorError::DivisorIsZero => {
            // n had been checked for being odd above, so should be unreachable, but play
            // safe.
            debug_assert!(false);
            CtMontgomeryTransformationError::InvalidModulus
        }
    })?;
    ct_mod_lshifted_mp_mp(t, t.len(), radix_shift_len, &n).map_err(|e| match e {
        CtModLshiftedMpMpError::InsufficientRemainderSpace => {
            // The result space had been checked at function entry already, but play safe.
            debug_assert!(false);
            CtMontgomeryTransformationError::InsufficientResultSpace
        }
    })?;
    Ok(())
}

/// Compute the square of the Montgomery radix modulo a multiprecision integer.
///
/// The multiprecision integers may be of different representations.
///
/// Compute *R<sup>2</sup> mod `n`*, where *R* is the [Montgomery
/// radix](ct_montgomery_radix_shift_len), and store it in `radix2_mod_n_out`,
/// with any of its bytes above `n`'s length cleared. The result is less than
/// `n`. It's to be computed once per modulus and serves as input to
/// [`ct_to_montgomery_form_mp()`] for converting integers to [Montgomery
/// form](self#montgomery-form).
///
/// Execution time depends only on the multiprecision integers' [lengths in
/// units of bytes](MpUIntCommon::len) as well as on `n`'s length without any
/// leading zero bytes, c.f. [`CtMpDivisor`], and is otherwise constant as far
/// as branching is concerned. The computation cost is O(n<sup>2</sup>), where
/// n is `n`'s [length](MpUIntCommon::len).
///
/// # Arguments
///
/// * `radix2_mod_n_out` - The output for the result. Its [length in units of
///   bytes](MpUIntCommon::len) must be such that `n`'s [fits
///   into](MpUIntCommon::len_is_compatible_with) it.
/// * `n` - The modulus. Must be odd.
///
/// # Errors
///
/// - Returns [`CtMontgomeryTransformationError::InvalidModulus`] if `n` is
///   even.
/// - Returns [`CtMontgomeryTransformationError::InsufficientResultSpace`] if
///   `n`'s length doesn't fit into `radix2_mod_n_out`'s.
///
/// Nothing gets modified in either case.
pub fn ct_montgomery_radix2_mod_n_mp<RX2T: MpMutUIntSlice, NT: MpUIntCommon>(
    radix2_mod_n_out: &mut RX2T,
    n: &NT,
) -> Result<(), CtMontgomeryTransformationError> {
    if n.test_bit(0).unwrap() == 0 {
        return Err(CtMontgomeryTransformationError::InvalidModulus);
    }
    if !n.len_is_compatible_with(radix2_mod_n_out.len()) {
        return Err(CtMontgomeryTransformationError::InsufficientResultSpace);
    }

    radix2_mod_n_out.clear_bytes_above(n.len());
    let mut radix2_mod_n_out = radix2_mod_n_out.shrink_to(n.len());
    debug_assert_eq!(radix2_mod_n_out.nlimbs(), n.nlimbs());

    let radix_shift_len = ct_montgomery_radix_shift_len(n.len());
    let n = CtMpDivisor::new(n, None).map_err(|e| match e {
        CtMpDivisorError::DivisorIsZero => {
            // n had been checked for being odd above, so should be unreachable, but play
            // safe.
            debug_assert!(false);
            CtMontgomeryTransformationError::InvalidModulus
        }
    })?;
    ct_mod_pow2_mp::<_, _>(2 * 8 * radix_shift_len, &mut radix2_mod_n_out, &n).map_err(
        |e| match e {
            CtModPow2MpError::InsufficientRemainderSpace => {
                // The result space had been checked at function entry already, but play safe.
                debug_assert!(false);
                CtMontgomeryTransformationError::InsufficientResultSpace
            }
        },
    )?;
    Ok(())
}

/// Error returned by [`ct_to_montgomery_form_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtToMontgomeryFormMpError {
    /// The modulus `n` is empty or even.
    InvalidModulus,
    /// The modulus' [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) the result's.
    InsufficientResultSpace,
    /// `t`'s [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) the modulus'.
    InconsistentOperandLengths,
    /// `radix2_mod_n`'s [length in units of bytes](MpUIntCommon::len) doesn't
    /// [fit into](MpUIntCommon::len_is_compatible_with) the modulus'.
    InconsistentRadix2ModNLenth,
}

/// Convert a multiprecision integer to Montgomery form by Montgomery
/// multiplication.
///
/// The multiprecision integers may be of different representations.
///
/// Store the [Montgomery form](self#montgomery-form) of `t`, *`t` · R mod `n`*,
/// where *R* is the [Montgomery radix](ct_montgomery_radix_shift_len), in
/// `result`. The result is less than `n`. The conversion is carried out by a
/// Montgomery multiplication of `t` with the precomputed
/// *R<sup>2</sup> mod `n`*, c.f. [`ct_montgomery_radix2_mod_n_mp()`], which is
/// more efficient than [`ct_to_montgomery_form_direct_mp()`] when converting
/// multiple integers for the same modulus. See [`ct_montgomery_redc_mp()`] for
/// converting back.
///
/// Execution time depends only on the multiprecision integers' [lengths in
/// units of bytes](MpUIntCommon::len) and is otherwise constant as far as
/// branching is concerned. The computation cost is O(n<sup>2</sup>), where n is
/// `n`'s [length](MpUIntCommon::len).
///
/// # Arguments
///
/// * `result` - The output for the result. Its [length in units of
///   bytes](MpUIntCommon::len) must be such that `n`'s [fits
///   into](MpUIntCommon::len_is_compatible_with) it.
/// * `t` - The multiprecision integer to convert. Must be less than `n`. Its
///   [length in units of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `n`'s.
/// * `n` - The modulus. Must be odd.
/// * `neg_n0_inv_mod_l` - The negated inverse of `n`'s least significant limb
///   modulo *2<sup>[`LimbType::BITS`]</sup>*, as computed by
///   [`ct_montgomery_neg_n0_inv_mod_l_mp()`].
/// * `radix2_mod_n` - *R<sup>2</sup> mod `n`*, as computed by
///   [`ct_montgomery_radix2_mod_n_mp()`]. Its [length in units of
///   bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `n`'s.
///
/// # Errors
///
/// - Returns [`CtToMontgomeryFormMpError::InvalidModulus`] if `n` is empty or
///   even.
/// - Returns [`CtToMontgomeryFormMpError::InsufficientResultSpace`] if `n`'s
///   length doesn't fit into `result`'s.
/// - Returns [`CtToMontgomeryFormMpError::InconsistentOperandLengths`] if `t`'s
///   length doesn't fit into `n`'s.
/// - Returns [`CtToMontgomeryFormMpError::InconsistentRadix2ModNLenth`] if
///   `radix2_mod_n`'s length doesn't fit into `n`'s.
///
/// Nothing gets modified in any of these cases.
pub fn ct_to_montgomery_form_mp<
    RT: MpMutUIntSlice,
    TT: MpUIntCommon,
    NT: MpUIntCommon,
    RX2T: MpUIntCommon,
>(
    result: &mut RT,
    t: &TT,
    n: &NT,
    neg_n0_inv_mod_l: LimbType,
    radix2_mod_n: &RX2T,
) -> Result<(), CtToMontgomeryFormMpError> {
    // The Montgomery multiplication will do all error checking needed. However, to
    // disambiguate which of the two factors has a length inconsistent with n,
    // if any, check that here.
    if !t.len_is_compatible_with(n.len()) {
        return Err(CtToMontgomeryFormMpError::InconsistentOperandLengths);
    }
    if !radix2_mod_n.len_is_compatible_with(n.len()) {
        return Err(CtToMontgomeryFormMpError::InconsistentRadix2ModNLenth);
    }
    debug_assert!(ct_lt_mp_mp(t, n).unwrap() != 0);
    debug_assert!(ct_lt_mp_mp(radix2_mod_n, n).unwrap() != 0);

    // All input arguments have been validated above, just unwrap().
    ct_montgomery_mul_mod_mp_mp(result, t, radix2_mod_n, n, neg_n0_inv_mod_l).map_err(
        |e| match e {
            CtMontgomeryMulModMpMpError::InsufficientResultSpace => {
                CtToMontgomeryFormMpError::InsufficientResultSpace
            }
            CtMontgomeryMulModMpMpError::InvalidModulus => {
                CtToMontgomeryFormMpError::InvalidModulus
            }
            CtMontgomeryMulModMpMpError::InconsistentOperandLengths => {
                // The multiplication's factors have been validated above, but play safe.
                CtToMontgomeryFormMpError::InconsistentOperandLengths
            }
        },
    )?;

    Ok(())
}

#[cfg(test)]
fn test_ct_to_montgomery_form_mp<TT: MpMutUIntSlice, NT: MpMutUIntSlice, RX2T: MpMutUIntSlice>() {
    use super::cmp::ct_eq_mp_mp;
    use super::div::ct_mod_mp_mp;
    use super::limb::LIMB_BYTES;

    for i in 0..16 {
        const MERSENNE_PRIME_13: LimbType = 8191 as LimbType;
        let n_high = MERSENNE_PRIME_13.wrapping_mul((65543 as LimbType).wrapping_mul(i));
        for j in 0..16 {
            const MERSENNE_PRIME_17: LimbType = 131071 as LimbType;
            let n_low = MERSENNE_PRIME_17.wrapping_mul((4095 as LimbType).wrapping_mul(j));
            // Force n_low odd.
            let n_low = n_low | 1;
            let mut n_buf = tst_mk_mp_backing_vec!(NT, 2 * LIMB_BYTES);
            let mut n = NT::from_slice(n_buf.as_mut_slice()).unwrap();
            n.store_l(0, n_low);
            n.store_l(1, n_high);
            drop(n);
            let n_lengths = if !NT::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
                [LIMB_BYTES, 2 * LIMB_BYTES]
            } else {
                [2 * LIMB_BYTES - 1, 2 * LIMB_BYTES]
            };

            for n_len in n_lengths {
                let mut n_buf = n_buf.clone();
                let mut n = NT::from_slice(n_buf.as_mut_slice()).unwrap();
                n.clear_bytes_above(n_len);
                let n = n.shrink_to(n_len);
                let neg_n0_inv = ct_montgomery_neg_n0_inv_mod_l_mp(&n).unwrap();

                let mut radix2_mod_n = tst_mk_mp_backing_vec!(RX2T, n_len);
                let mut radix2_mod_n = RX2T::from_slice(radix2_mod_n.as_mut_slice()).unwrap();
                ct_montgomery_radix2_mod_n_mp(&mut radix2_mod_n, &n).unwrap();

                for k in 0..4 {
                    let a_high =
                        MERSENNE_PRIME_17.wrapping_mul((16383 as LimbType).wrapping_mul(k));
                    for l in 0..4 {
                        let a_low =
                            MERSENNE_PRIME_13.wrapping_mul((262175 as LimbType).wrapping_mul(l));
                        let mut a = tst_mk_mp_backing_vec!(TT, 2 * LIMB_BYTES);
                        let mut a = TT::from_slice(a.as_mut_slice()).unwrap();
                        a.store_l(0, a_low);
                        a.store_l(1, a_high);
                        // All montgomery operations are defined mod n, compute a mod n
                        ct_mod_mp_mp(None, &mut a, &CtMpDivisor::new(&n, None).unwrap());
                        let mut a = a.shrink_to(n_len);

                        let mut result = tst_mk_mp_backing_vec!(TT, n_len);
                        let mut result = TT::from_slice(result.as_mut_slice()).unwrap();
                        ct_to_montgomery_form_mp(&mut result, &a, &n, neg_n0_inv, &radix2_mod_n)
                            .unwrap();

                        ct_to_montgomery_form_direct_mp(&mut a, &n).unwrap();
                        assert_eq!(ct_eq_mp_mp(&result, &a).unwrap(), 1);
                    }
                }
            }
        }
    }
}

#[test]
fn test_ct_to_montgomery_form_be_be_be() {
    use super::limbs::MpMutBigEndianUIntByteSlice;
    test_ct_to_montgomery_form_mp::<
        MpMutBigEndianUIntByteSlice,
        MpMutBigEndianUIntByteSlice,
        MpMutBigEndianUIntByteSlice,
    >()
}

#[test]
fn test_ct_to_montgomery_form_le_le_le() {
    use super::limbs::MpMutLittleEndianUIntByteSlice;
    test_ct_to_montgomery_form_mp::<
        MpMutLittleEndianUIntByteSlice,
        MpMutLittleEndianUIntByteSlice,
        MpMutLittleEndianUIntByteSlice,
    >()
}

#[test]
fn test_ct_to_montgomery_form_ne_ne_ne() {
    use super::limbs::MpMutNativeEndianUIntLimbsSlice;
    test_ct_to_montgomery_form_mp::<
        MpMutNativeEndianUIntLimbsSlice,
        MpMutNativeEndianUIntLimbsSlice,
        MpMutNativeEndianUIntLimbsSlice,
    >()
}

/// Montgomery exponentiation kernel.
///
/// Compute the [Montgomery form](self#montgomery-form) of `op0`'s value raised
/// to the power of `exponent` modulo `n`, for `op0` given in [Montgomery
/// form](self#montgomery-form), and store it in `result`. `result` must be
/// initialized to the [Montgomery form](self#montgomery-form) of one, i.e. to
/// *R mod `n`*, with *R* being the [Montgomery
/// radix](ct_montgomery_radix_shift_len).
///
/// Only the `exponent_nbits` least significant bits of `exponent` are
/// considered, with `exponent_nbits` capped to `exponent`'s bit length. The
/// computation follows the left-to-right binary method, with the
/// multiplication by `op0` carried out for every bit and undone by a
/// conditional copy for zero bits, so that the execution time doesn't depend
/// on `exponent`'s value.
///
/// Execution time depends only on the multiprecision integers' [lengths in
/// units of bytes](MpUIntCommon::len) and on `exponent_nbits`, and is otherwise
/// constant as far as branching is concerned. The computation cost is
/// O(e · n<sup>2</sup>), where e is `exponent_nbits` and n is `n`'s
/// [length](MpUIntCommon::len). The arguments are validated by the callers,
/// [`ct_montogmery_exp_mod_odd_mp_mp()`] and [`ct_exp_mod_odd_mp_mp()`].
///
/// # Arguments
///
/// * `result` - The output for the result. Must be initialized to the
///   [Montgomery form](self#montgomery-form) of one, *R mod `n`*, and span
///   exactly as many limbs as `n`.
/// * `op0` - The base, in [Montgomery form](self#montgomery-form). Must be less
///   than `n`. Its [length in units of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `n`'s.
/// * `n` - The modulus. Must be odd.
/// * `neg_n0_inv_mod_l` - The negated inverse of `n`'s least significant limb
///   modulo *2<sup>[`LimbType::BITS`]</sup>*, as computed by
///   [`ct_montgomery_neg_n0_inv_mod_l_mp()`].
/// * `exponent` - The exponent.
/// * `exponent_nbits` - The number of `exponent`'s least significant bits to
///   consider.
/// * `scratch` - Scratch buffer. Must have at least as many [`LimbType`]
///   elements as `n` spans limbs, c.f.
///   [`nlimbs_for_len()`](crate::MpNativeEndianUIntLimbsSlice::nlimbs_for_len).
fn _ct_montogmery_exp_mod_mp_mp<
    RT: MpMutUIntSlice,
    T0: MpUIntCommon,
    NT: MpUIntCommon,
    ET: MpUIntCommon,
>(
    result: &mut RT,
    op0: &T0,
    n: &NT,
    neg_n0_inv_mod_l: LimbType,
    exponent: &ET,
    exponent_nbits: usize,
    scratch: &mut [LimbType],
) {
    debug_assert_eq!(result.nlimbs(), n.nlimbs());

    let n_nlimbs = MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(n.len());
    debug_assert!(scratch.len() >= n_nlimbs);
    let mut scratch = MpMutNativeEndianUIntLimbsSlice::from_limbs(scratch);
    scratch.clear_bytes_above(n.len());
    let mut scratch = scratch.shrink_to(n.len());

    let exponent_nbits = exponent_nbits.min(8 * exponent.len());
    for i in 0..exponent_nbits {
        // Input arguments have been validated/setup by callers, just unwrap() the
        // result.
        ct_montgomery_mul_mod_mp_mp(&mut scratch, result, result, n, neg_n0_inv_mod_l).unwrap();
        ct_montgomery_mul_mod_mp_mp(result, &scratch, op0, n, neg_n0_inv_mod_l).unwrap();
        // If the current exponent bit is zero, "undo" the latter multiplication.
        result.copy_from_cond(&scratch, !exponent.test_bit(exponent_nbits - i - 1));
    }
}

/// Error returned by [`ct_montogmery_exp_mod_odd_mp_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtMontgomeryExpModOddMpMpError {
    /// The modulus `n` is even, including the case of an empty one.
    InvalidModulus,
    /// The modulus' [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) `result`'s.
    InsufficientResultSpace,
    /// `scratch` has fewer [`LimbType`] elements than the modulus spans limbs.
    InsufficientScratchSpace,
    /// `op0`'s [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) the modulus'.
    InconsistentOperandLengths,
    /// `radix_mod_n`'s [length in units of bytes](MpUIntCommon::len) doesn't
    /// [fit into](MpUIntCommon::len_is_compatible_with) the modulus'.
    InconsistendRadixModNLengh,
}

/// Modular exponentiation in Montgomery form.
///
/// The multiprecision integers may be of different representations.
///
/// Compute the [Montgomery form](self#montgomery-form) of `op0`'s value raised
/// to the power of `exponent` modulo an odd `n`, for `op0` given in [Montgomery
/// form](self#montgomery-form), and store it in `result`, with any of
/// `result`'s bytes above `n`'s length cleared. That is, for
/// *`op0` = a · R mod `n`*, with *R* being the [Montgomery
/// radix](ct_montgomery_radix_shift_len), compute
/// *a<sup>`exponent`</sup> · R mod `n`*. The result is less than `n`. Only the
/// `exponent_nbits` least significant bits of `exponent` are considered, which
/// allows for not revealing the exponent's actual bit length.
///
/// The exponentiation follows the left-to-right binary method with Montgomery
/// multiplications, c.f. [`ct_montgomery_mul_mod_mp_mp()`], carrying out a
/// multiplication by `op0` for every exponent bit, so that the execution time
/// doesn't depend on `exponent`'s value. See [`ct_exp_mod_odd_mp_mp()`] for a
/// variant operating on, and returning, integers not in [Montgomery
/// form](self#montgomery-form).
///
/// Execution time depends only on the multiprecision integers' [lengths in
/// units of bytes](MpUIntCommon::len), on `scratch`'s length and on
/// `exponent_nbits`, capped to `exponent`'s bit length, and is otherwise
/// constant as far as branching is concerned. The computation cost is O(e ·
/// n<sup>2</sup>), where e is `exponent_nbits` and n is `n`'s
/// [length](MpUIntCommon::len).
///
/// # Arguments
///
/// * `result` - The output for the result. Its [length in units of
///   bytes](MpUIntCommon::len) must be such that `n`'s [fits
///   into](MpUIntCommon::len_is_compatible_with) it.
/// * `op0` - The base, in [Montgomery form](self#montgomery-form). Must be less
///   than `n`. Its [length in units of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `n`'s.
/// * `n` - The modulus. Must be odd.
/// * `neg_n0_inv_mod_l` - The negated inverse of `n`'s least significant limb
///   modulo *2<sup>[`LimbType::BITS`]</sup>*, as computed by
///   [`ct_montgomery_neg_n0_inv_mod_l_mp()`].
/// * `radix_mod_n` - The [Montgomery form](self#montgomery-form) of one, i.e.
///   *R mod `n`*, as obtained e.g. by Montgomery reducing the result of
///   [`ct_montgomery_radix2_mod_n_mp()`] with [`ct_montgomery_redc_mp()`]. Its
///   [length in units of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `n`'s.
/// * `exponent` - The exponent. Its [length in units of
///   bytes](MpUIntCommon::len) isn't restricted.
/// * `exponent_nbits` - The number of `exponent`'s least significant bits to
///   consider. Values exceeding `exponent`'s bit length are capped to the
///   latter.
/// * `scratch` - Scratch buffer. Must have at least as many [`LimbType`]
///   elements as `n` spans limbs, c.f.
///   [`nlimbs_for_len()`](crate::MpNativeEndianUIntLimbsSlice::nlimbs_for_len).
///
/// # Errors
///
/// - Returns [`CtMontgomeryExpModOddMpMpError::InvalidModulus`] if `n` is
///   even.
/// - Returns [`CtMontgomeryExpModOddMpMpError::InsufficientResultSpace`] if
///   `n`'s length doesn't fit into `result`'s.
/// - Returns [`CtMontgomeryExpModOddMpMpError::InsufficientScratchSpace`] if
///   `scratch` is too short.
/// - Returns [`CtMontgomeryExpModOddMpMpError::InconsistentOperandLengths`] if
///   `op0`'s length doesn't fit into `n`'s.
/// - Returns [`CtMontgomeryExpModOddMpMpError::InconsistendRadixModNLengh`] if
///   `radix_mod_n`'s length doesn't fit into `n`'s.
///
/// Nothing gets modified in any of these cases.
#[allow(clippy::too_many_arguments)]
pub fn ct_montogmery_exp_mod_odd_mp_mp<
    RT: MpMutUIntSlice,
    T0: MpUIntCommon,
    NT: MpUIntCommon,
    RXT: MpUIntCommon,
    ET: MpUIntCommon,
>(
    result: &mut RT,
    op0: &T0,
    n: &NT,
    neg_n0_inv_mod_l: LimbType,
    radix_mod_n: &RXT,
    exponent: &ET,
    exponent_nbits: usize,
    scratch: &mut [LimbType],
) -> Result<(), CtMontgomeryExpModOddMpMpError> {
    if n.test_bit(0).unwrap() == 0 {
        return Err(CtMontgomeryExpModOddMpMpError::InvalidModulus);
    }
    if !n.len_is_compatible_with(result.len()) {
        return Err(CtMontgomeryExpModOddMpMpError::InsufficientResultSpace);
    }
    if scratch.len() < MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(n.len()) {
        return Err(CtMontgomeryExpModOddMpMpError::InsufficientScratchSpace);
    }
    if !op0.len_is_compatible_with(n.len()) {
        return Err(CtMontgomeryExpModOddMpMpError::InconsistentOperandLengths);
    }
    debug_assert!(ct_lt_mp_mp(op0, n).unwrap() != 0);
    if !radix_mod_n.len_is_compatible_with(n.len()) {
        return Err(CtMontgomeryExpModOddMpMpError::InconsistendRadixModNLengh);
    }
    debug_assert!(ct_lt_mp_mp(radix_mod_n, n).unwrap() != 0);

    // Initialize the result with a one in Montgomery form.
    result.clear_bytes_above(n.len());
    let mut result = result.shrink_to(n.len());
    debug_assert_eq!(result.nlimbs(), n.nlimbs());
    result.copy_from(radix_mod_n);

    _ct_montogmery_exp_mod_mp_mp(
        &mut result,
        op0,
        n,
        neg_n0_inv_mod_l,
        exponent,
        exponent_nbits,
        scratch,
    );
    Ok(())
}

/// Error returned by [`ct_exp_mod_odd_mp_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtExpModOddMpMpError {
    /// The modulus `n` is empty or even.
    InvalidModulus,
    /// The modulus' [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) `result`'s.
    InsufficientResultSpace,
    /// `scratch` has fewer [`LimbType`] elements than the modulus spans limbs.
    InsufficientScratchSpace,
    /// `op0`'s and the modulus' [lengths in units of bytes](MpUIntCommon::len)
    /// don't [fit into](MpUIntCommon::len_is_compatible_with) each other
    /// mutually.
    InconsistentOperandLengths,
}

/// Modular exponentiation with an odd modulus.
///
/// The multiprecision integers may be of different representations.
///
/// Compute *`op0`<sup>`exponent`</sup> mod `n`* for an odd `n` and store it in
/// `result`, with any of `result`'s bytes above `n`'s length cleared. The
/// result is less than `n`. Only the `exponent_nbits` least significant bits of
/// `exponent` are considered, which allows for not revealing the exponent's
/// actual bit length.
///
/// The exponentiation is carried out by means of Montgomery multiplications:
/// `op0` gets converted to [Montgomery form](self#montgomery-form), c.f.
/// [`ct_to_montgomery_form_mp()`], and the result back from it, c.f.
/// [`ct_montgomery_redc_mp()`]. To save a buffer, `op0` gets overwritten by its
/// [Montgomery form](self#montgomery-form) in the course. See
/// [`ct_montogmery_exp_mod_odd_mp_mp()`] for a variant operating on, and
/// returning, integers in [Montgomery form](self#montgomery-form), which avoids
/// the conversions, e.g. for multiple exponentiations with the same modulus.
///
/// Execution time depends only on the multiprecision integers' [lengths in
/// units of bytes](MpUIntCommon::len), on `scratch`'s length, on
/// `exponent_nbits`, capped to `exponent`'s bit length, as well as on `n`'s
/// length without any leading zero bytes, c.f. [`CtMpDivisor`], and is
/// otherwise constant as far as branching is concerned. In particular, it
/// doesn't depend on `exponent`'s value. The computation cost is
/// O((e + 1) · n<sup>2</sup>), where e is `exponent_nbits` and n is `n`'s
/// [length](MpUIntCommon::len).
///
/// # Arguments
///
/// * `result` - The output for the result. Its [length in units of
///   bytes](MpUIntCommon::len) must be such that `n`'s [fits
///   into](MpUIntCommon::len_is_compatible_with) it.
/// * `op0` - The base. Must be less than `n`. Its [length in units of
///   bytes](MpUIntCommon::len) and `n`'s must [fit
///   into](MpUIntCommon::len_is_compatible_with) each other mutually. It will
///   be overwritten by its [Montgomery form](self#montgomery-form).
/// * `n` - The modulus. Must be odd.
/// * `exponent` - The exponent. Its [length in units of
///   bytes](MpUIntCommon::len) isn't restricted.
/// * `exponent_nbits` - The number of `exponent`'s least significant bits to
///   consider. Values exceeding `exponent`'s bit length are capped to the
///   latter.
/// * `scratch` - Scratch buffer. Must have at least as many [`LimbType`]
///   elements as `n` spans limbs, c.f.
///   [`nlimbs_for_len()`](crate::MpNativeEndianUIntLimbsSlice::nlimbs_for_len).
///
/// # Errors
///
/// - Returns [`CtExpModOddMpMpError::InvalidModulus`] if `n` is empty or even.
/// - Returns [`CtExpModOddMpMpError::InsufficientResultSpace`] if `n`'s length
///   doesn't fit into `result`'s.
/// - Returns [`CtExpModOddMpMpError::InsufficientScratchSpace`] if `scratch`
///   is too short.
/// - Returns [`CtExpModOddMpMpError::InconsistentOperandLengths`] if `op0`'s
///   and `n`'s lengths don't fit into each other mutually.
///
/// Nothing gets modified in any of these cases.
pub fn ct_exp_mod_odd_mp_mp<
    RT: MpMutUIntSlice,
    T0: MpMutUInt,
    NT: MpUIntCommon,
    ET: MpUIntCommon,
>(
    result: &mut RT,
    op0: &mut T0,
    n: &NT,
    exponent: &ET,
    exponent_nbits: usize,
    scratch: &mut [LimbType],
) -> Result<(), CtExpModOddMpMpError> {
    if !n.len_is_compatible_with(result.len()) {
        return Err(CtExpModOddMpMpError::InsufficientResultSpace);
    }
    if !op0.len_is_compatible_with(n.len()) {
        return Err(CtExpModOddMpMpError::InconsistentOperandLengths);
    }
    debug_assert!(ct_lt_mp_mp(op0, n).unwrap() != 0);
    if !n.len_is_compatible_with(op0.len()) {
        // op0 will get transformed in-place into Montgomery form. So the
        // backing byte slice must be large enough.
        return Err(CtExpModOddMpMpError::InconsistentOperandLengths);
    }
    debug_assert_eq!(op0.nlimbs(), n.nlimbs());
    if scratch.len() < MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(n.len()) {
        return Err(CtExpModOddMpMpError::InsufficientScratchSpace);
    }

    // This checks the modulus for validity as a side-effect.
    let neg_n0_inv_mod_l = ct_montgomery_neg_n0_inv_mod_l_mp(n).map_err(|e| match e {
        CtMontgomeryNegN0InvModLMpError::InvalidModulus => CtExpModOddMpMpError::InvalidModulus,
    })?;

    // Shrink result[] to the length of n. It will be used to temporarily hold the
    // radix^2 mod n and ct_to_montgomery_form_mp() below would complain if its
    // length is unexpectedly large.
    result.clear_bytes_above(n.len());
    let mut result = result.shrink_to(n.len());
    debug_assert_eq!(result.nlimbs(), n.nlimbs());

    // The radix squared mod n gets into result[], it will be reduced
    // later on to a one in Montgomery form.
    ct_montgomery_radix2_mod_n_mp(&mut result, n).unwrap();

    // Transform op0 into Montgomery form, the function argument will get
    // overwritten to save an extra scratch buffer.
    let mut mg_op0 = MpMutNativeEndianUIntLimbsSlice::from_limbs(scratch);
    mg_op0.clear_bytes_above(n.len());
    let mut mg_op0 = mg_op0.shrink_to(n.len());
    ct_to_montgomery_form_mp(&mut mg_op0, op0, n, neg_n0_inv_mod_l, &result).unwrap();
    op0.copy_from(&mg_op0);

    // Reduce the radix squared mod n in result[] to the radix mod n,
    // i.e. to a one in Montgomery form.
    ct_montgomery_redc_mp(&mut result, n, neg_n0_inv_mod_l).unwrap();

    // Do the Montgomery exponentiation.
    _ct_montogmery_exp_mod_mp_mp(
        &mut result,
        op0,
        n,
        neg_n0_inv_mod_l,
        exponent,
        exponent_nbits,
        scratch,
    );

    // And transform the result back from Montgomery form.
    ct_montgomery_redc_mp(&mut result, n, neg_n0_inv_mod_l).unwrap();
    Ok(())
}

#[cfg(test)]
fn test_ct_exp_mod_odd_mp_mp<
    RT: MpMutUIntSlice,
    T0: MpMutUIntSlice,
    NT: MpMutUIntSlice,
    ET: MpMutUIntSlice,
>() {
    extern crate alloc;
    use super::limb::LIMB_BYTES;
    use super::mul::ct_mul_trunc_mp_l;
    use super::shift::ct_lshift_mp;
    use alloc::vec;

    fn test_one<RT: MpMutUIntSlice, T0: MpMutUIntSlice, NT: MpUIntCommon, ET: MpUIntCommon>(
        op0: &T0,
        n: &NT,
        exponent: &ET,
    ) {
        use super::cmp::ct_eq_mp_mp;
        use super::div::ct_mod_mp_mp;
        use super::mul::{ct_mul_trunc_mp_mp, ct_square_trunc_mp};

        let n_len = n.len();

        let mut op0_mod_n = tst_mk_mp_backing_vec!(T0, n_len.max(op0.len()));
        let mut op0_mod_n = T0::from_slice(&mut op0_mod_n).unwrap();
        op0_mod_n.copy_from(op0);
        ct_mod_mp_mp(None, &mut op0_mod_n, &CtMpDivisor::new(n, None).unwrap());

        let mut op0_scratch = tst_mk_mp_backing_vec!(T0, n_len);
        let mut op0_scratch = T0::from_slice(&mut op0_scratch).unwrap();
        op0_scratch.copy_from(&op0_mod_n);
        let mut result = tst_mk_mp_backing_vec!(RT, n_len);
        let mut result = RT::from_slice(&mut result).unwrap();
        let mut scratch =
            vec![0 as LimbType; MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(n_len)];
        ct_exp_mod_odd_mp_mp(
            &mut result,
            &mut op0_scratch,
            n,
            exponent,
            8 * exponent.len(),
            &mut scratch,
        )
        .unwrap();

        // Compute the expected value using repeated multiplications/squarings and
        // modular reductions.
        let mut expected = tst_mk_mp_backing_vec!(RT, 2 * n_len);
        let mut expected = RT::from_slice(&mut expected).unwrap();
        expected.clear_bytes_above(0);
        expected.store_l(0, 1);
        for i in 0..8 * exponent.len() {
            ct_square_trunc_mp(&mut expected, n_len);
            ct_mod_mp_mp(None, &mut expected, &CtMpDivisor::new(n, None).unwrap());
            if exponent.test_bit(8 * exponent.len() - i - 1).unwrap() != 0 {
                ct_mul_trunc_mp_mp(&mut expected, n_len, &op0_mod_n);
                ct_mod_mp_mp(None, &mut expected, &CtMpDivisor::new(n, None).unwrap());
            }
        }
        assert_ne!(ct_eq_mp_mp(&result, &expected).unwrap(), 0);
    }

    let exponent_len = LIMB_BYTES + 1;
    let mut e0_buf = tst_mk_mp_backing_vec!(ET, exponent_len);
    let mut e1_buf = tst_mk_mp_backing_vec!(ET, exponent_len);
    let mut e1 = ET::from_slice(&mut e1_buf).unwrap();
    e1.store_l(0, 1);
    drop(e1);
    let mut e2_buf = tst_mk_mp_backing_vec!(ET, exponent_len);
    let mut e2 = ET::from_slice(&mut e2_buf).unwrap();
    e2.store_l(0, 2);
    drop(e2);
    let mut ef_buf = tst_mk_mp_backing_vec!(ET, exponent_len);
    ef_buf.fill(0xefu8.into());

    for n_len in [1, LIMB_BYTES + 1, 2 * LIMB_BYTES - 1, 3 * LIMB_BYTES] {
        let mut n = tst_mk_mp_backing_vec!(NT, n_len);
        let mut n = NT::from_slice(&mut n).unwrap();
        n.store_l(0, 1);
        while n.load_l((n_len - 1) / LIMB_BYTES) >> (8 * ((n_len - 1) % LIMB_BYTES)) == 0 {
            ct_mul_trunc_mp_l(&mut n, n_len, 251);
        }

        for op0_len in 1..n_len {
            let mut op0 = tst_mk_mp_backing_vec!(T0, n_len);
            let mut op0 = T0::from_slice(&mut op0).unwrap();
            op0.store_l(0, 1);
            for _ in 0..op0_len {
                ct_mul_trunc_mp_l(&mut op0, n_len, 241);
            }
            ct_lshift_mp(&mut op0, 8 * (n_len - op0_len));
            for e_buf in [&mut e0_buf, &mut e1_buf, &mut e2_buf, &mut ef_buf] {
                let e = ET::from_slice(e_buf).unwrap();
                test_one::<RT, _, _, _>(&op0, &n, &e);
            }
        }
    }
}

#[test]
fn test_ct_exp_mod_odd_be_be_be_be() {
    use super::limbs::MpMutBigEndianUIntByteSlice;
    test_ct_exp_mod_odd_mp_mp::<
        MpMutBigEndianUIntByteSlice,
        MpMutBigEndianUIntByteSlice,
        MpMutBigEndianUIntByteSlice,
        MpMutBigEndianUIntByteSlice,
    >()
}

#[test]
fn test_ct_exp_mod_odd_le_le_le_le() {
    use super::limbs::MpMutLittleEndianUIntByteSlice;
    test_ct_exp_mod_odd_mp_mp::<
        MpMutLittleEndianUIntByteSlice,
        MpMutLittleEndianUIntByteSlice,
        MpMutLittleEndianUIntByteSlice,
        MpMutLittleEndianUIntByteSlice,
    >()
}

#[test]
fn test_ct_exp_mod_odd_ne_ne_ne_ne() {
    use super::limbs::MpMutNativeEndianUIntLimbsSlice;
    test_ct_exp_mod_odd_mp_mp::<
        MpMutNativeEndianUIntLimbsSlice,
        MpMutNativeEndianUIntLimbsSlice,
        MpMutNativeEndianUIntLimbsSlice,
        MpMutNativeEndianUIntLimbsSlice,
    >()
}