cmpa 0.1.6

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

//! Implementation of multiprecision integer multiplication primitives.

use super::limb::{
    LIMB_BITS, LimbChoice, LimbType, black_box_l, ct_add_l_l, ct_mul_add_l_l_l_c, ct_mul_l_l,
};
#[cfg(test)]
use super::limbs::MpMutUIntSlice;
use super::limbs::{MpMutUInt, MpUIntCommon, ct_mp_nlimbs};

/// Conditionally multiply two multiprecision integers.
///
/// The two multiprecision integers may be of different representations.
///
/// If `cond` is true, replace `op0` by the product of `op0`'s value, as given
/// by its `op0_in_len` least significant bytes, and `op1`. If the product
/// doesn't fit into `op0`'s [length in units of bytes](MpUIntCommon::len), its
/// most significant part gets truncated. That is, only the product's
/// `op0.len()` least significant bytes are retained. As the product is less
/// than two to the power of the number of bits in `op0_in_len + op1.len()`
/// bytes, truncation is avoided by having `op0.len() >= op0_in_len +
/// op1.len()`.
///
/// Otherwise, if `cond` is false, `op0`'s value as given by its `op0_in_len`
/// least significant bytes remains unchanged, but any of `op0`'s bytes above
/// those get cleared nonetheless.
///
/// Execution time depends only on `op0_in_len` and on the multiprecision
/// integers' [lengths in units of bytes](MpUIntCommon::len) and is otherwise
/// constant as far as branching is concerned. In particular, it doesn't depend
/// on `cond`. The computation cost is O(m · n), where m is `op0_in_len` and n
/// is `op0`'s [length](MpUIntCommon::len). See [`ct_mul_trunc_mp_mp()`] for
/// the unconditional variant, [`ct_square_trunc_mp()`] for squaring and
/// [`ct_mul_trunc_mp_l()`] for multiplying by a single [`LimbType`].
///
/// # Arguments
///
/// * `op0` - The first factor, given by its `op0_in_len` least significant
///   bytes. Any higher bytes are ignored. It will be overwritten by the
///   resulting, possibly truncated product if `cond` is set.
/// * `op0_in_len` - The first factor's length in units of bytes within `op0`.
///   Must be at most `op0`'s [length in units of bytes](MpUIntCommon::len). As
///   `op0` also receives the product, it usually needs to be allocated larger
///   than the first factor alone would require, hence the latter's length must
///   be specified separately.
/// * `op1` - The second factor. Its [length in units of
///   bytes](MpUIntCommon::len) isn't restricted.
/// * `cond` - Whether or not to replace `op0` by the product. Intended for
///   constant-time implementations of algorithms conditionally multiplying,
///   like e.g. binary exponentiation.
pub fn ct_mul_trunc_cond_mp_mp<T0: MpMutUInt, T1: MpUIntCommon>(
    op0: &mut T0,
    op0_in_len: usize,
    op1: &T1,
    cond: LimbChoice,
) {
    debug_assert!(op0_in_len <= op0.len());
    op0.clear_bytes_above(op0_in_len);
    // If op1's length is zero, interpret that as a zero.
    if op1.is_empty() {
        let cond_mask = cond.select(0, !0);
        for j in 0..op0.nlimbs() {
            op0.store_l(j, op0.load_l(j) & !cond_mask);
        }
        return;
    }
    let op1_nlimbs = op1.nlimbs();

    let result_high_mask = op0.partial_high_mask();
    let op0_nlimbs = op0.nlimbs();
    let op0_in_nlimbs = ct_mp_nlimbs(op0_in_len);

    let mut j = op0_in_nlimbs;
    while j > 0 {
        j -= 1;

        let op0_val = op0.load_l(j);
        op0.store_l(j, 0);

        let mut carry = 0;
        let result_nlimbs = op0_nlimbs - j;
        let mut cond_unit = cond.select(1, 0);
        for k in 0..op1_nlimbs.min(result_nlimbs) {
            let op1_val = cond.select(0, op1.load_l(k)) | cond_unit;
            cond_unit = 0;

            let mut result_val = op0.load_l(j + k);
            (carry, result_val) = ct_mul_add_l_l_l_c(result_val, op0_val, op1_val, carry);

            if k != result_nlimbs - 1 || !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
                op0.store_l_full(j + k, result_val);
            } else {
                op0.store_l(j + k, result_val & result_high_mask);
            }
        }
        // Propagate the carry all the way up.
        for k in op1_nlimbs..result_nlimbs {
            let mut result_val = op0.load_l(j + k);
            (carry, result_val) = ct_add_l_l(result_val, carry);
            if k != result_nlimbs - 1 || !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
                op0.store_l_full(j + k, result_val);
            } else {
                op0.store_l(j + k, result_val & result_high_mask);
            }
        }
    }
}

#[cfg(test)]
fn test_ct_mul_trunc_cond_mp_mp<T0: MpMutUIntSlice, T1: MpMutUIntSlice>() {
    use super::limb::LIMB_BYTES;

    let mut op0 = tst_mk_mp_backing_vec!(T0, 5 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    let mut op1 = tst_mk_mp_backing_vec!(T1, 2 * LIMB_BYTES);
    let mut op1 = T1::from_slice(&mut op1).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    op1.store_l(0, !0);
    op1.store_l(1, !0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(0));
    assert_eq!(op0.load_l(0), !0);
    assert_eq!(op0.load_l(1), !0);
    assert_eq!(op0.load_l(2), 0);
    assert_eq!(op0.load_l(3), 0);
    assert_eq!(op0.load_l(4), 0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(1));
    assert_eq!(op0.load_l(0), 1);
    assert_eq!(op0.load_l(1), 0);
    assert_eq!(op0.load_l(2), !1);
    assert_eq!(op0.load_l(3), !0);
    assert_eq!(op0.load_l(4), 0);

    let mut op0 = tst_mk_mp_backing_vec!(T0, 3 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    let mut op1 = tst_mk_mp_backing_vec!(T1, 2 * LIMB_BYTES);
    let mut op1 = T1::from_slice(&mut op1).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    op1.store_l(0, !0);
    op1.store_l(1, !0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(0));
    assert_eq!(op0.load_l(0), !0);
    assert_eq!(op0.load_l(1), !0);
    assert_eq!(op0.load_l(2), 0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(1));
    assert_eq!(op0.load_l(0), 1);
    assert_eq!(op0.load_l(1), 0);
    assert_eq!(op0.load_l(2), !1);

    let mut op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    let mut op1 = tst_mk_mp_backing_vec!(T1, 2 * LIMB_BYTES);
    let mut op1 = T1::from_slice(&mut op1).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    op1.store_l(0, !0);
    op1.store_l(1, !0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(0));
    assert_eq!(op0.load_l(0), !0);
    assert_eq!(op0.load_l(1), !0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(1));
    assert_eq!(op0.load_l(0), 1);
    assert_eq!(op0.load_l(1), 0);

    // An empty op1 is interpreted as a zero. If cond is unset, op0's bytes above
    // op0_in_len still get cleared, consistent with the non-empty op1 case.
    let mut op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    let mut op1 = tst_mk_mp_backing_vec!(T1, 0);
    let op1 = T1::from_slice(&mut op1).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    ct_mul_trunc_cond_mp_mp(&mut op0, LIMB_BYTES, &op1, LimbChoice::from(0));
    assert_eq!(op0.load_l(0), !0);
    assert_eq!(op0.load_l(1), 0);
    op0.store_l(1, !0);
    ct_mul_trunc_cond_mp_mp(&mut op0, LIMB_BYTES, &op1, LimbChoice::from(1));
    assert_eq!(op0.load_l(0), 0);
    assert_eq!(op0.load_l(1), 0);

    if !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
        return;
    }

    let mut op0 = tst_mk_mp_backing_vec!(T0, 4 * LIMB_BYTES - 1);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    let mut op1 = tst_mk_mp_backing_vec!(T1, 2 * LIMB_BYTES);
    let mut op1 = T1::from_slice(&mut op1).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    op1.store_l(0, !0);
    op1.store_l(1, !0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(0));
    assert_eq!(op0.load_l(0), !0);
    assert_eq!(op0.load_l(1), !0);
    assert_eq!(op0.load_l(2), 0);
    assert_eq!(op0.load_l(3), 0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(1));
    assert_eq!(op0.load_l(0), 1);
    assert_eq!(op0.load_l(1), 0);
    assert_eq!(op0.load_l(2), !1);
    assert_eq!(op0.load_l(3), !0 >> 8);

    let mut op0 = tst_mk_mp_backing_vec!(T0, 3 * LIMB_BYTES - 1);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    let mut op1 = tst_mk_mp_backing_vec!(T1, 2 * LIMB_BYTES);
    let mut op1 = T1::from_slice(&mut op1).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    op1.store_l(0, !0);
    op1.store_l(1, !0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(0));
    assert_eq!(op0.load_l(0), !0);
    assert_eq!(op0.load_l(1), !0);
    assert_eq!(op0.load_l(2), 0);
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES, &op1, LimbChoice::from(1));
    assert_eq!(op0.load_l(0), 1);
    assert_eq!(op0.load_l(1), 0);
    assert_eq!(op0.load_l(2), (!0 >> 8) ^ 1);

    let mut op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES - 1);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    let mut op1 = tst_mk_mp_backing_vec!(T1, 2 * LIMB_BYTES);
    let mut op1 = T1::from_slice(&mut op1).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0 >> (2 * 8));
    op1.store_l(0, !0);
    op1.store_l(1, !0 >> (2 * 8));
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES - 1, &op1, LimbChoice::from(0));
    assert_eq!(op0.load_l(0), !0);
    assert_eq!(op0.load_l(1), !0 >> (2 * 8));
    ct_mul_trunc_cond_mp_mp(&mut op0, 2 * LIMB_BYTES - 1, &op1, LimbChoice::from(1));
    assert_eq!(op0.load_l(0), 1);
    assert_eq!(op0.load_l(1), 0xfe << (8 * (LIMB_BYTES - 2)));
}

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

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

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

/// Multiply two multiprecision integers.
///
/// The two multiprecision integers may be of different representations.
///
/// Replace `op0` by the product of `op0`'s value, as given by its `op0_in_len`
/// least significant bytes, and `op1`. If the product doesn't fit into `op0`'s
/// [length in units of bytes](MpUIntCommon::len), its most significant part
/// gets truncated. That is, only the product's `op0.len()` least significant
/// bytes are retained. As the product is less than two to the power of the
/// number of bits in `op0_in_len + op1.len()` bytes, truncation is avoided by
/// having `op0.len() >= op0_in_len + op1.len()`.
///
/// Execution time depends only on `op0_in_len` and 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 `op0_in_len` and n is `op0`'s
/// [length](MpUIntCommon::len). See [`ct_mul_trunc_cond_mp_mp()`] for the
/// conditional variant, [`ct_square_trunc_mp()`] for squaring and
/// [`ct_mul_trunc_mp_l()`] for multiplying by a single [`LimbType`].
///
/// # Arguments
///
/// * `op0` - The first factor, given by its `op0_in_len` least significant
///   bytes. Any higher bytes are ignored. It will be overwritten by the
///   resulting, possibly truncated product.
/// * `op0_in_len` - The first factor's length in units of bytes within `op0`.
///   Must be at most `op0`'s [length in units of bytes](MpUIntCommon::len). As
///   `op0` also receives the product, it usually needs to be allocated larger
///   than the first factor alone would require, hence the latter's length must
///   be specified separately.
/// * `op1` - The second factor. Its [length in units of
///   bytes](MpUIntCommon::len) isn't restricted.
pub fn ct_mul_trunc_mp_mp<T0: MpMutUInt, T1: MpUIntCommon>(
    op0: &mut T0,
    op0_in_len: usize,
    op1: &T1,
) {
    ct_mul_trunc_cond_mp_mp(op0, op0_in_len, op1, LimbChoice::from(1))
}

/// Square a multiprecision integer.
///
/// Replace `op0` by the square of its value, as given by its `op0_in_len` least
/// significant bytes. If the square doesn't fit into `op0`'s [length in units
/// of bytes](MpUIntCommon::len), its most significant part gets truncated.
/// That is, only the square's `op0.len()` least significant bytes are retained.
/// As the square is less than two to the power of the number of bits in
/// `2 * op0_in_len` bytes, truncation is avoided by having
/// `op0.len() >= 2 * op0_in_len`.
///
/// Execution time depends only on `op0_in_len` and on `op0`'s [length 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 `op0_in_len` and
/// n is `op0`'s [length](MpUIntCommon::len). This is a specialized variant of
/// [`ct_mul_trunc_mp_mp()`] for squaring, which saves about half of the limb
/// multiplications.
///
/// # Arguments
///
/// * `op0` - The operand to square, given by its `op0_in_len` least significant
///   bytes. Any higher bytes are ignored. It will be overwritten by the
///   resulting, possibly truncated square.
/// * `op0_in_len` - The operand's length in units of bytes within `op0`. Must
///   be at most `op0`'s [length in units of bytes](MpUIntCommon::len). As
///   `op0` also receives the square, it usually needs to be allocated larger
///   than the operand alone would require, hence the latter's length must be
///   specified separately.
pub fn ct_square_trunc_mp<T0: MpMutUInt>(op0: &mut T0, op0_in_len: usize) {
    debug_assert!(op0_in_len <= op0.len());
    let result_high_mask = op0.partial_high_mask();
    let op0_nlimbs = op0.nlimbs();
    op0.clear_bytes_above(op0_in_len);
    let op0_in_nlimbs = ct_mp_nlimbs(op0_in_len);

    let mut j = op0_in_nlimbs;
    while j > 0 {
        j -= 1;

        let op0_val = op0.load_l(j);
        op0.store_l(j, 0);

        let mut last_prod_high: LimbType = 0;
        let mut carry = 0;
        let result_nlimbs = op0_nlimbs - j;
        // (u[j] + u[j-1:0])^2 = u[j]^2 + 2 * u[j] * u[j-1:0]) + u[j-1:0]^2.
        // Account for the middle term first.
        for k in 0..j.min(result_nlimbs) {
            // Being the upper half of a multiplication result,
            // last_prod_high <= !1 always holds trivially.
            // As will be shown below, even the invariant
            // carry <= 2 || (2 < carry <= 4 && last_prod_high <= !2)
            // holds at loop entry.
            let op1_val = op0.load_l(k);
            let prod = ct_mul_l_l(op0_val, op1_val);
            let mut result_val = op0.load_l(j + k);

            // Multiply last_prod_high, the upper half of the last iteration's
            // multiplication, by two.
            let carry0 = black_box_l(last_prod_high >> (LIMB_BITS - 1));
            last_prod_high = last_prod_high.wrapping_mul(2);
            // From the loop invariant, it follows that
            // - if carry <= 2, then last_prod_high <= !3 and
            // - if 2 < carry <= 4, then last_prod_high <= !5.
            // In either case, addition of carry to last_prod_high
            // does not overflow and the sum is <= !1.
            debug_assert!(last_prod_high <= !3);
            debug_assert!(carry <= 2 || last_prod_high <= !5);
            last_prod_high += carry;
            // If the sum below wraps around, then by virtue of the fact that
            // last_prod_high <= !1, the result will be <= !2. That is
            // if carry1 != 0, then result_val <= !2
            let carry1;
            (carry1, result_val) = ct_add_l_l(result_val, last_prod_high);

            // Done with last_prod_high, store away the half part to account for in the next
            // iteration. Note that as a basic property of the multiplication,
            // prod.high() is always <= LimbType::MAX - 1. Moreover, either
            // prod.high() is even <= LimbType::MAX - 2 or prod.low() <= 1.
            last_prod_high = prod.high();

            // Multiply the lower part of the current iteration's multiplication by two.
            // If prod.high() == !1, i.e. at the maximum possible value, then
            // the scaled lower half <= 2 * 1 == 2, c.f. the remark above. That is,
            // prod.high() <= !2 || (carry2 == 0 && prod_low <= 2).
            let carry2 = black_box_l(prod.low() >> (LIMB_BITS - 1));
            let prod_low = prod.low().wrapping_mul(2);

            // Add the scaled prod_low to the result.
            // To prove that the loop invariant does indeed hold, assume
            // that last_prod_high (== prod.high()) > !2. It needs to be
            // shown that the carry calculated below as
            // carry = carry0 + carry1 + carry2 + carry3 is <= 2.
            //
            // From the remark right above the preceeding step, the
            // assumption prod.high() > !2 implies that
            // carry2 == 0 && prod_low <= 2.
            //
            // Consider the following two cases:
            // a.) carry1 != 0:
            //     If carry1 != 0, then by the remark right before
            //     the computation of (carry1, result_val) above,
            //     result_val <= !2 at this point. As prod_low <= 2,
            //     the sum of result_val and prod_low below will not
            //     overflow, from which it follows that the associated
            //     carry3 == 0. Thus, carry2 == 0 && carry3 == 0 and
            //     the sum carry = carry0 + carry1 + carry2 + carry3 <= 2.
            // b.) carry1 == 0:
            //     Both carry1 and carry2 are zero, from which it follows
            //     trivially that
            //     carry = carry0 + carry1 + carry2 + carry3 <= 2.
            let carry3;
            (carry3, result_val) = ct_add_l_l(result_val, prod_low);
            carry = carry0 + carry1 + carry2 + carry3;
            // Confirm the loop invariant.
            debug_assert!(last_prod_high <= !1);
            debug_assert!(carry <= 2 || last_prod_high <= !2);

            if k != result_nlimbs - 1 || !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
                op0.store_l_full(j + k, result_val);
            } else {
                op0.store_l(j + k, result_val & result_high_mask);
            }
        }

        // Now handle the u[j]^2 part of the quadratic expansion.
        if j >= result_nlimbs {
            continue;
        }
        let prod = ct_mul_l_l(op0_val, op0_val);
        let mut result_val = op0.load_l(2 * j);
        // Multiply last_prod_high from the previous loop's last iteration by two.
        let carry0 = black_box_l(last_prod_high >> (LIMB_BITS - 1));
        last_prod_high = last_prod_high.wrapping_mul(2);
        // From the previous loop's invariant, it again follows that the addition of
        // carry to last_prod_high does not overflow the sum.
        debug_assert!(last_prod_high <= !3);
        debug_assert!(carry <= 2 || last_prod_high <= !5);
        last_prod_high += carry;
        let carry1;
        (carry1, result_val) = ct_add_l_l(result_val, last_prod_high);
        last_prod_high = prod.high();
        let carry2;
        (carry2, result_val) = ct_add_l_l(result_val, prod.low());
        carry = carry0 + carry1 + carry2;
        if j != result_nlimbs - 1 {
            op0.store_l_full(2 * j, result_val);
        } else {
            op0.store_l(2 * j, result_val & result_high_mask);
        }

        // Propagate the carry all the way up. The first iteration will also account for
        // the previous multiplications upper limb.
        for k in j + 1..result_nlimbs {
            let mut result_val = op0.load_l(j + k);
            let carry0;
            (carry0, result_val) = ct_add_l_l(result_val, last_prod_high);
            last_prod_high = 0;
            let carry1;
            (carry1, result_val) = ct_add_l_l(result_val, carry);
            carry = carry0 + carry1;
            if k != result_nlimbs - 1 || !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
                op0.store_l_full(j + k, result_val);
            } else {
                op0.store_l(j + k, result_val & result_high_mask);
            }
        }
    }
}

#[cfg(test)]
fn test_ct_square_trunc_mp<T0: MpMutUIntSlice>() {
    extern crate alloc;
    use super::limb::LIMB_BYTES;
    use alloc::vec::Vec;

    fn square_by_mul<T0: MpMutUIntSlice>(
        op0: &[T0::BackingSliceElementType],
        op0_in_len: usize,
    ) -> Vec<T0::BackingSliceElementType> {
        let mut _result = Vec::from(op0);
        let mut result = T0::from_slice(&mut _result).unwrap();
        let mut op0 = Vec::from(op0);
        let mut op0 = T0::from_slice(&mut op0).unwrap();
        let op0 = op0.shrink_to(op0_in_len);
        ct_mul_trunc_mp_mp(&mut result, op0_in_len, &op0);
        drop(result);
        _result
    }

    let mut _op0 = tst_mk_mp_backing_vec!(T0, 5 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    drop(op0);
    let expected = square_by_mul::<T0>(&_op0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    ct_square_trunc_mp(&mut op0, 2 * LIMB_BYTES);
    drop(op0);
    assert_eq!(_op0, expected);

    let mut _op0 = tst_mk_mp_backing_vec!(T0, 3 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    drop(op0);
    let expected = square_by_mul::<T0>(&_op0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    ct_square_trunc_mp(&mut op0, 2 * LIMB_BYTES);
    drop(op0);
    assert_eq!(_op0, expected);

    let mut _op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    drop(op0);
    let expected = square_by_mul::<T0>(&_op0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    ct_square_trunc_mp(&mut op0, 2 * LIMB_BYTES);
    drop(op0);
    assert_eq!(_op0, expected);

    if !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
        return;
    }

    let mut _op0 = tst_mk_mp_backing_vec!(T0, 4 * LIMB_BYTES - 1);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    drop(op0);
    let expected = square_by_mul::<T0>(&_op0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    ct_square_trunc_mp(&mut op0, 2 * LIMB_BYTES);
    drop(op0);
    assert_eq!(_op0, expected);

    let mut _op0 = tst_mk_mp_backing_vec!(T0, 3 * LIMB_BYTES - 1);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    drop(op0);
    let expected = square_by_mul::<T0>(&_op0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    ct_square_trunc_mp(&mut op0, 2 * LIMB_BYTES);
    drop(op0);
    assert_eq!(_op0, expected);

    let mut _op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES - 1);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0 >> (2 * 8));
    drop(op0);
    let expected = square_by_mul::<T0>(&_op0, 2 * LIMB_BYTES - 1);
    let mut op0 = T0::from_slice(&mut _op0).unwrap();
    ct_square_trunc_mp(&mut op0, 2 * LIMB_BYTES - 1);
    drop(op0);
    assert_eq!(_op0, expected);
}

#[test]
fn test_ct_square_trunc_be() {
    use super::limbs::MpMutBigEndianUIntByteSlice;
    test_ct_square_trunc_mp::<MpMutBigEndianUIntByteSlice>()
}

#[test]
fn test_ct_square_trunc_le() {
    use super::limbs::MpMutLittleEndianUIntByteSlice;
    test_ct_square_trunc_mp::<MpMutLittleEndianUIntByteSlice>()
}

#[test]
fn test_ct_square_trunc_ne() {
    use super::limbs::MpMutNativeEndianUIntLimbsSlice;
    test_ct_square_trunc_mp::<MpMutNativeEndianUIntLimbsSlice>()
}

/// Multiply a multiprecision integer by a [`LimbType`].
///
/// Replace `op0` by the product of `op0`'s value, as given by its `op0_in_len`
/// least significant bytes, and `op1`. If the product doesn't fit into `op0`'s
/// [length in units of bytes](MpUIntCommon::len), its most significant part
/// gets truncated. That is, only the product's `op0.len()` least significant
/// bytes are retained. As the product is less than two to the power of the
/// number of bits in `op0_in_len` bytes plus one [`LimbType`], truncation is
/// avoided by having `op0.len()` exceed `op0_in_len` by at least the
/// [`LimbType`] size in bytes.
///
/// Execution time depends only on `op0_in_len` and on `op0`'s [length in units
/// of bytes](MpUIntCommon::len) and is otherwise constant as far as branching
/// is concerned. This is a specialized variant of [`ct_mul_trunc_mp_mp()`] and
/// [`ct_mul_trunc_cond_mp_mp()`] for a single [`LimbType`] second factor.
///
/// # Arguments
///
/// * `op0` - The first factor, given by its `op0_in_len` least significant
///   bytes. Any higher bytes are ignored. It will be overwritten by the
///   resulting, possibly truncated product.
/// * `op0_in_len` - The first factor's length in units of bytes within `op0`.
///   Must be at most `op0`'s [length in units of bytes](MpUIntCommon::len). As
///   `op0` also receives the product, it usually needs to be allocated larger
///   than the first factor alone would require, hence the latter's length must
///   be specified separately.
/// * `op1` - The second factor.
pub fn ct_mul_trunc_mp_l<T0: MpMutUInt>(op0: &mut T0, op0_in_len: usize, op1: LimbType) {
    debug_assert!(op0_in_len <= op0.len());
    let result_high_mask = op0.partial_high_mask();
    let op0_nlimbs = op0.nlimbs();
    op0.clear_bytes_above(op0_in_len);
    let op0_in_nlimbs = ct_mp_nlimbs(op0_in_len);

    if op0_in_len == 0 {
        return;
    }

    let mut carry = 0;
    for j in 0..op0_in_nlimbs {
        let op0_val = op0.load_l(j);
        let result_val;
        (carry, result_val) = ct_mul_add_l_l_l_c(0, op0_val, op1, carry);

        if j != op0_nlimbs - 1 || !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
            op0.store_l_full(j, result_val);
        } else {
            op0.store_l(j, result_val & result_high_mask);
        }
    }

    if op0_in_nlimbs != op0_nlimbs {
        if op0_in_nlimbs != op0_nlimbs - 1 || !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
            op0.store_l_full(op0_in_nlimbs, carry);
        } else {
            op0.store_l(op0_in_nlimbs, carry & result_high_mask);
        }
    }
}

#[cfg(test)]
fn test_ct_mul_trunc_mp_l<T0: MpMutUIntSlice>() {
    use super::limb::LIMB_BYTES;

    let mut op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    let op1 = 0;
    ct_mul_trunc_mp_l(&mut op0, 2 * LIMB_BYTES, op1);
    assert_eq!(op0.load_l(0), 0);
    assert_eq!(op0.load_l(1), 0);

    let mut op0 = tst_mk_mp_backing_vec!(T0, 3 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    let op1 = 2;
    ct_mul_trunc_mp_l(&mut op0, 2 * LIMB_BYTES, op1);
    assert_eq!(op0.load_l(0), !1);
    assert_eq!(op0.load_l(1), !0);
    assert_eq!(op0.load_l(2), 1);

    let mut op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    let op1 = 2;
    ct_mul_trunc_mp_l(&mut op0, 2 * LIMB_BYTES, op1);
    assert_eq!(op0.load_l(0), !1);
    assert_eq!(op0.load_l(1), !0);

    if !T0::SUPPORTS_UNALIGNED_BUFFER_LENGTHS {
        return;
    }

    let mut op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES + 1);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0);
    let op1 = 2;
    ct_mul_trunc_mp_l(&mut op0, 2 * LIMB_BYTES, op1);
    assert_eq!(op0.load_l(0), !1);
    assert_eq!(op0.load_l(1), !0);
    assert_eq!(op0.load_l(2), 1);

    let mut op0 = tst_mk_mp_backing_vec!(T0, 2 * LIMB_BYTES - 1);
    let mut op0 = T0::from_slice(&mut op0).unwrap();
    op0.store_l(0, !0);
    op0.store_l(1, !0 >> 8);
    let op1 = 2;
    ct_mul_trunc_mp_l(&mut op0, 2 * LIMB_BYTES - 1, op1);
    assert_eq!(op0.load_l(0), !1);
    assert_eq!(op0.load_l(1), !0 >> 8);
}

#[test]
fn test_ct_mul_trunc_be_l() {
    use super::limbs::MpMutBigEndianUIntByteSlice;
    test_ct_mul_trunc_mp_l::<MpMutBigEndianUIntByteSlice>()
}

#[test]
fn test_ct_mul_trunc_le_l() {
    use super::limbs::MpMutLittleEndianUIntByteSlice;
    test_ct_mul_trunc_mp_l::<MpMutLittleEndianUIntByteSlice>()
}

#[test]
fn test_ct_mul_trunc_ne_l() {
    use super::limbs::MpMutNativeEndianUIntLimbsSlice;
    test_ct_mul_trunc_mp_l::<MpMutNativeEndianUIntLimbsSlice>()
}