cmpa 0.1.6

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

//! Primality testing and prime candidate generation.
//!
//! The primitives here support the generation of large (probable) primes,
//! e.g. for cryptographic key generation:
//! - [`PrimeWheelSieve`] enumerates prime candidates not divisible by any of
//!   the first few primes, starting from a given lower bound.
//! - [`ct_composite_test_small_prime_gcd_mp()`] cheaply identifies candidates
//!   with a small prime factor as composite.
//! - [`ct_prime_test_miller_rabin_mp()`] carries out a round of the
//!   Miller-Rabin probabilistic primality test.

use super::cmp::{ct_is_one_mp, ct_leq_mp_l, ct_lt_mp_mp};
use super::div::ct_mod_mp_l;
use super::euclid::ct_gcd_odd_mp_mp;
use super::hexstr;
use super::limb::{
    CtLDivisor, LIMB_BYTES, LimbChoice, LimbType, ct_eq_l_l, ct_geq_l_l, ct_is_zero_l, ct_lt_l_l,
    ct_sub_l_l_b,
};
use super::limbs::{
    MpBigEndianUIntByteSlice, MpMutNativeEndianUIntLimbsSlice, MpMutUInt, MpMutUIntSlice,
    MpUIntCommon, MpUIntSlicePriv as _, ct_find_first_set_bit_mp,
};
use super::montgomery::{
    CtMontgomeryNegN0InvModLMpError, ct_montgomery_mul_mod_mp_mp, ct_montgomery_neg_n0_inv_mod_l_mp,
};
use super::usize_ct_cmp::{ct_eq_usize_usize, ct_lt_usize_usize};

/// Product of first 3 primes > 2, i.e. from 3 to 7 (inclusive).
///
/// Filters ~54% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_8: [u8; 1] = hexstr::bytes_from_hexstr_cnst::<1>("69");

/// Product of first 5 primes > 2, i.e. from 3 to 13 (inclusive).
///
/// Filters ~62% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_16: [u8; 2] = hexstr::bytes_from_hexstr_cnst::<2>("3aa7");

/// Product of first 9 primes > 2, i.e. from 3 to 29 (inclusive).
///
/// Filters ~68% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_32: [u8; 4] = hexstr::bytes_from_hexstr_cnst::<4>("c0cfd797");

/// Product of first 15 primes > 2, i.e. from 3 to 53 (inclusive).
///
/// Filters ~73% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_64: [u8; 8] = hexstr::bytes_from_hexstr_cnst::<8>("e221f97c30e94e1d");

/// Product of first 25 primes > 2, i.e. from 3 to 101 (inclusive).
///
/// Filters ~76% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_128: [u8; 16] =
    hexstr::bytes_from_hexstr_cnst::<16>("5797d47c51681549d734e4fc4c3eaf7f");

/// Product of first 43 primes > 2, i.e. from 3 to 193 (inclusive).
///
/// Filters ~79% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_256: [u8; 32] = hexstr::bytes_from_hexstr_cnst::<32>(
    "dbf05b6f5654b3c0f524355143958688\
     9f155887819aed2ac05b93352be98677",
);

/// Product of first 74 primes > 2, i.e. from 3 to 379 (inclusive).
///
/// Filters ~81% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_512: [u8; 64] = hexstr::bytes_from_hexstr_cnst::<64>(
    "106aa9fb7646fa6eb0813c28c5d5f09f\
     077ec3ba238bfb99c1b631a203e81187\
     233db117cbc384056ef04659a4a11de4\
     9f7ecb29bada8f980decece92e30c48f",
);

/// Product of first 130 primes > 2, i.e. from 3 to 739 (inclusive).
///
/// Filters ~83% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_1024: [u8; 128] = hexstr::bytes_from_hexstr_cnst::<128>(
    "02c85ff870f24be80f62b1ba6c20bd72\
     b837efdf121206d87db56b7d69fa4c02\
     1c107c3ca206fe8fa7080ef576effc82\
     f9b10f5750656b7794b16afd70996e91\
     aef6e0ad15e91b071ac9b24d98b233ad\
     86ee055518e58e56638ef18bac5c74cb\
     35bbb6e5dae2783dd1c0ce7dec4fc70e\
     5186d411df36368f061aa36011f30179",
);

/// Product of first 232 primes > 2, i.e. from 3 to 1471 (inclusive).
///
/// Filters ~85% of odd prime candidates. Big-endian.
const SMALL_ODD_PRIME_PRODUCT_2048: [u8; 256] = hexstr::bytes_from_hexstr_cnst::<256>(
    "2465a7bd85011e1c9e0527929fff268c\
     82ef7efa416863baa5acdb0971dba0cc\
     ac3ee4999345029f2cf810b99e406aac\
     5fce5dd69d1c717daea5d18ab913f456\
     505679bc91c57d46d9888857862b36e2\
     ede2e473c1f0ab359da25271affe15ff\
     240e299d0b04f4cd0e4d7c0e47b1a7ba\
     007de89aae848fd5bdcd7f9815564eb0\
     60ae14f19cb50c291f0bbd8ed1c4c7f8\
     fc5fba51662001939b532d92dac844a8\
     431d400c832d039f5f900b278a75219c\
     2986140c79045d7759540854c31504dc\
     56f1df5eebe7bee447658b917bf696d6\
     927f2e2428fbeb340e515cb9835d6387\
     1be8bbe09cf13445799f2e6778815157\
     1a93b4c1eee55d1b9072e0b2f5c4607f",
);

/// The products of the first few odd primes, ordered by increasing length.
///
/// [`ct_composite_test_small_prime_gcd_mp()`] selects the longest one fitting
/// the candidate's length.
const SMALL_ODD_PRIME_PRODUCTS: [&[u8]; 9] = [
    &SMALL_ODD_PRIME_PRODUCT_8,
    &SMALL_ODD_PRIME_PRODUCT_16,
    &SMALL_ODD_PRIME_PRODUCT_32,
    &SMALL_ODD_PRIME_PRODUCT_64,
    &SMALL_ODD_PRIME_PRODUCT_128,
    &SMALL_ODD_PRIME_PRODUCT_256,
    &SMALL_ODD_PRIME_PRODUCT_512,
    &SMALL_ODD_PRIME_PRODUCT_1024,
    &SMALL_ODD_PRIME_PRODUCT_2048,
];

/// Error returned by [`ct_composite_test_small_prime_gcd_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtCompositeTestSmallPrimeGcdMpError {
    /// Any of the `scratch` buffers has fewer [`LimbType`] elements than `p`
    /// spans limbs.
    InsufficientScratchSpace,
}

/// Test a prime candidate for small prime factors.
///
/// Returns a [`LimbChoice`] representing `true` if `p` is divisible by two or
/// by any of the small odd primes in a precomputed product, and `false`
/// otherwise. A `true` result means
/// that `p` is composite, unless `p` is one of these small primes itself,
/// which are reported as `true` as well, as are zero and an empty `p`. A
/// `false` result doesn't imply primality, the test serves as a cheap
/// prefilter to [`ct_prime_test_miller_rabin_mp()`].
///
/// The divisibility gets tested by computing the greatest common divisor of
/// `p` with the longest of the precomputed small odd prime products fitting
/// `p`'s [length in units of bytes](MpUIntCommon::len), rounded up to a
/// multiple of the [`LimbType`] size, c.f. [`ct_gcd_odd_mp_mp()`]. The
/// longer `p`, the more small primes get tested for.
///
/// Execution time depends only on `p`'s [length 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 `p`'s
/// [length](MpUIntCommon::len).
///
/// # Arguments
///
/// * `p` - The prime candidate.
/// * `scratch` - Two scratch buffers. Each must have at least as many
///   [`LimbType`] elements as `p` spans limbs, c.f.
///   [`nlimbs_for_len()`](crate::MpNativeEndianUIntLimbsSlice::nlimbs_for_len).
///
/// # Errors
///
/// Returns [`CtCompositeTestSmallPrimeGcdMpError::InsufficientScratchSpace`] if
/// any of the `scratch` buffers is too short.
pub fn ct_composite_test_small_prime_gcd_mp<PT: MpUIntCommon>(
    p: &PT,
    scratch: [&mut [LimbType]; 2],
) -> Result<LimbChoice, CtCompositeTestSmallPrimeGcdMpError> {
    if p.is_empty() {
        return Ok(LimbChoice::from(1));
    }

    let p_nlimbs = MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(p.len());
    for s in scratch.iter() {
        if s.len() < p_nlimbs {
            return Err(CtCompositeTestSmallPrimeGcdMpError::InsufficientScratchSpace);
        }
    }
    let p_len_aligned = p_nlimbs * LIMB_BYTES;

    // The GCD runtime depends on the maximum of both operands' bit length.
    // Select the largest small prime product with length <= p.len().
    let mut i = SMALL_ODD_PRIME_PRODUCTS.len();
    while i > 0 {
        i -= 1;
        if SMALL_ODD_PRIME_PRODUCTS[i].len() <= p_len_aligned {
            break;
        }
    }
    let small_prime_product = SMALL_ODD_PRIME_PRODUCTS[i];

    let [scratch0, scratch1] = scratch;

    let mut gcd = MpMutNativeEndianUIntLimbsSlice::from_limbs(scratch0);
    gcd.clear_bytes_above(p_len_aligned);
    let mut gcd = gcd.shrink_to(p_len_aligned);
    gcd.copy_from(&MpBigEndianUIntByteSlice::from_slice(small_prime_product).unwrap());
    let mut p_work_scratch = MpMutNativeEndianUIntLimbsSlice::from_limbs(scratch1);
    p_work_scratch.copy_from(p);
    let mut p_work_scratch = p_work_scratch.shrink_to(p_len_aligned);

    ct_gcd_odd_mp_mp(&mut gcd, &mut p_work_scratch).unwrap();
    let gcd_is_one = ct_is_one_mp(&gcd);

    // The small prime products don't include a factor of two.
    // Test for it separately.
    let p_is_odd = p.test_bit(0);

    Ok(!gcd_is_one | !p_is_odd)
}

#[cfg(test)]
fn test_ct_composite_test_small_prime_gcd_mp<PT: MpMutUIntSlice>() {
    extern crate alloc;
    use super::add::ct_sub_mp_l;
    use alloc::vec;

    for (i, small_odd_prime_product) in SMALL_ODD_PRIME_PRODUCTS.iter().enumerate().take(8) {
        let l = small_odd_prime_product.len();
        let p_len = if l != 1 { l + 1 } else { l };
        let mut p = tst_mk_mp_backing_vec!(PT, p_len);
        let mut p = PT::from_slice(&mut p).unwrap();
        let scratch_nlimbs = MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(p_len);
        let mut scratch0 = vec![0 as LimbType; scratch_nlimbs];
        let mut scratch1 = vec![0 as LimbType; scratch_nlimbs];

        p.clear_bytes_above(0);
        let scratch = [scratch0.as_mut_slice(), scratch1.as_mut_slice()];
        assert!(
            ct_composite_test_small_prime_gcd_mp(&p, scratch)
                .unwrap()
                .unwrap()
                != 0
        );

        p.store_l(0, 2);
        let scratch = [scratch0.as_mut_slice(), scratch1.as_mut_slice()];
        assert!(
            ct_composite_test_small_prime_gcd_mp(&p, scratch)
                .unwrap()
                .unwrap()
                != 0
        );

        p.store_l(0, 3);
        let scratch = [scratch0.as_mut_slice(), scratch1.as_mut_slice()];
        assert!(
            ct_composite_test_small_prime_gcd_mp(&p, scratch)
                .unwrap()
                .unwrap()
                != 0
        );

        let j = if i > 0 { i - 1 } else { 0 };
        p.copy_from(&MpBigEndianUIntByteSlice::from_slice(SMALL_ODD_PRIME_PRODUCTS[j]).unwrap());
        let scratch = [scratch0.as_mut_slice(), scratch1.as_mut_slice()];
        assert!(
            ct_composite_test_small_prime_gcd_mp(&p, scratch)
                .unwrap()
                .unwrap()
                != 0
        );
    }

    // p = 2^255 - 19.
    let p_len = 256 / 8;
    let mut p = tst_mk_mp_backing_vec!(PT, p_len);
    let mut p = PT::from_slice(&mut p).unwrap();
    p.set_bit_to(255, true);
    ct_sub_mp_l(&mut p, 19);
    let scratch_nlimbs = MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(p_len);
    let mut scratch0 = vec![0 as LimbType; scratch_nlimbs];
    let mut scratch1 = vec![0 as LimbType; scratch_nlimbs];
    let scratch = [scratch0.as_mut_slice(), scratch1.as_mut_slice()];
    assert_eq!(
        ct_composite_test_small_prime_gcd_mp(&p, scratch)
            .unwrap()
            .unwrap(),
        0
    );
}

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

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

#[test]
fn test_ct_composite_test_small_prime_gcd_ne() {
    test_ct_composite_test_small_prime_gcd_mp::<MpMutNativeEndianUIntLimbsSlice>()
}

/// Error returned by [`ct_prime_test_miller_rabin_mp()`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum CtPrimeTestMillerRabinMpError {
    /// The candidate `p` is even, including the case of an empty one.
    InvalidCandidate,
    /// Any of the `scratch` buffers has fewer [`LimbType`] elements than `p`
    /// spans limbs.
    InsufficientScratchSpace,
    /// `mg_base`'s [length in units of bytes](MpUIntCommon::len) doesn't [fit
    /// into](MpUIntCommon::len_is_compatible_with) `p`'s.
    InconsistentBaseLength,
    /// `mg_radix_mod_p`'s [length in units of bytes](MpUIntCommon::len) doesn't
    /// [fit into](MpUIntCommon::len_is_compatible_with) `p`'s.
    InconsistentMgRadixModPLength,
}

/// Carry out a round of the Miller-Rabin primality test.
///
/// Test whether the odd prime candidate `p` is a strong probable prime to the
/// base given in [Montgomery form](crate::montgomery#montgomery-form) by
/// `mg_base`. Returns a [`LimbChoice`] representing `false` if `p` is
/// definitely composite, and `true` if `p` passed the test, i.e. is a probable
/// prime. For a composite `p` and a uniformly random base, the probability of
/// the latter is at most 1/4, so the test is to be repeated for a number of
/// independent random bases, after having filtered candidates with small prime
/// factors by means of [`ct_composite_test_small_prime_gcd_mp()`].
///
/// The base's powers are computed by means of Montgomery multiplications, c.f.
/// [`ct_montgomery_mul_mod_mp_mp()`], with the exponent *`p` - 1* scanned over
/// all of `p`'s bit positions, independent of their values.
///
/// 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. In particular, it doesn't depend on whether `p`
/// passes the test. The computation cost is O(n<sup>3</sup>), where n is `p`'s
/// [length](MpUIntCommon::len).
///
/// # Arguments
///
/// * `mg_base` - The base, in [Montgomery
///   form](crate::montgomery#montgomery-form) with respect to `p`. It's usually
///   chosen uniformly at random from *[2, `p` - 2]* before conversion. Must be
///   less than `p`. Its [length in units of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `p`'s.
/// * `p` - The prime candidate. Must be odd and greater than one.
/// * `mg_radix_mod_p` - The [Montgomery
///   form](crate::montgomery#montgomery-form) of one, i.e. *R mod `p`*, with
///   *R* being the [Montgomery radix](crate::ct_montgomery_radix_shift_len), as
///   obtained e.g. by Montgomery reducing the result of
///   [`ct_montgomery_radix2_mod_n_mp()`](crate::ct_montgomery_radix2_mod_n_mp)
///   with [`ct_montgomery_redc_mp()`](crate::ct_montgomery_redc_mp). Must be
///   less than `p`. Its [length in units of bytes](MpUIntCommon::len) must [fit
///   into](MpUIntCommon::len_is_compatible_with) `p`'s.
/// * `scratch` - Two scratch buffers. Each must have at least as many
///   [`LimbType`] elements as `p` spans limbs, c.f.
///   [`nlimbs_for_len()`](crate::MpNativeEndianUIntLimbsSlice::nlimbs_for_len).
///
/// # Errors
///
/// - Returns [`CtPrimeTestMillerRabinMpError::InvalidCandidate`] if `p` is
///   even.
/// - Returns [`CtPrimeTestMillerRabinMpError::InconsistentBaseLength`] if
///   `mg_base`'s length doesn't fit into `p`'s.
/// - Returns [`CtPrimeTestMillerRabinMpError::InconsistentMgRadixModPLength`]
///   if `mg_radix_mod_p`'s length doesn't fit into `p`'s.
/// - Returns [`CtPrimeTestMillerRabinMpError::InsufficientScratchSpace`] if
///   any of the `scratch` buffers is too short.
pub fn ct_prime_test_miller_rabin_mp<BT: MpUIntCommon, PT: MpUIntCommon, RXT: MpUIntCommon>(
    mg_base: &BT,
    p: &PT,
    mg_radix_mod_p: &RXT,
    scratch: [&mut [LimbType]; 2],
) -> Result<LimbChoice, CtPrimeTestMillerRabinMpError> {
    if p.test_bit(0).unwrap() == 0 {
        return Err(CtPrimeTestMillerRabinMpError::InvalidCandidate);
    }
    if !mg_base.len_is_compatible_with(p.len()) {
        return Err(CtPrimeTestMillerRabinMpError::InconsistentBaseLength);
    }
    debug_assert_ne!(ct_lt_mp_mp(mg_base, p).unwrap(), 0);
    if !mg_radix_mod_p.len_is_compatible_with(p.len()) {
        return Err(CtPrimeTestMillerRabinMpError::InconsistentMgRadixModPLength);
    }
    debug_assert_ne!(ct_lt_mp_mp(mg_radix_mod_p, p).unwrap(), 0);

    let p_nlimbs = MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(p.len());
    for s in scratch.iter() {
        if s.len() < p_nlimbs {
            return Err(CtPrimeTestMillerRabinMpError::InsufficientScratchSpace);
        }
    }
    let [scratch0, scratch1] = scratch;

    let mut p_minus_one = MpMutNativeEndianUIntLimbsSlice::from_limbs(scratch0);
    p_minus_one.copy_from(p);
    let mut p_minus_one = p_minus_one.shrink_to(p.len());
    p_minus_one.store_l(0, p_minus_one.load_l(0) ^ 1); // Minus one for odd p.
    let (p_minus_one_is_nonzero, p_minus_one_first_set_bit_pos) =
        ct_find_first_set_bit_mp(&p_minus_one);
    debug_assert!(p_minus_one_is_nonzero.unwrap() != 0);
    debug_assert!(p_minus_one_first_set_bit_pos > 0);

    let mut base_pow = MpMutNativeEndianUIntLimbsSlice::from_limbs(scratch0);
    // Initialize with a 1 in Montgomery form.
    base_pow.copy_from(mg_radix_mod_p);
    let mut base_pow = base_pow.shrink_to(p.len());
    let mut pow_scratch = MpMutNativeEndianUIntLimbsSlice::from_limbs(scratch1);
    pow_scratch.clear_bytes_above(p.len());
    let mut pow_scratch = pow_scratch.shrink_to(p.len());

    let neg_p0_inv_mod_l = ct_montgomery_neg_n0_inv_mod_l_mp(p).map_err(|e| match e {
        CtMontgomeryNegN0InvModLMpError::InvalidModulus => {
            // Should not be possible, p had been validated for being an odd number at the
            // beginning, but play safe.
            debug_assert!(false);
            CtPrimeTestMillerRabinMpError::InvalidCandidate
        }
    })?;

    let mut is_probable_prime = 0;
    let mut exp_bit_pos = 8 * p.len();
    while exp_bit_pos > 1 {
        exp_bit_pos -= 1;
        ct_montgomery_mul_mod_mp_mp(&mut pow_scratch, &base_pow, &base_pow, p, neg_p0_inv_mod_l)
            .unwrap();
        ct_montgomery_mul_mod_mp_mp(&mut base_pow, &pow_scratch, mg_base, p, neg_p0_inv_mod_l)
            .unwrap();
        // If the current exponent bit is zero, "undo" the latter multiplication.
        // The exponent is p - 1, p is odd, so the lsb of p - 1 is zero and there's no
        // borrow into the next bit position. exp_bit_pos != 0 by the loop's
        // condition, so simply load from the oiginal p itself.
        base_pow.copy_from_cond(&pow_scratch, !p.test_bit(exp_bit_pos));

        // Compare the current power of the base against 1 and -1 (in Montgomery form).
        let mut is_one = 1;
        let mut is_minus_one = 1;
        let mut neg_one_borrow = 0;
        for i in 0..p.nlimbs() {
            let one_val = mg_radix_mod_p.load_l(i);
            let base_pow_val = base_pow.load_l(i);
            is_one &= ct_is_zero_l(base_pow_val ^ one_val);
            let p_val = p.load_l(i);
            let minus_one_val;
            (neg_one_borrow, minus_one_val) = ct_sub_l_l_b(p_val, one_val, neg_one_borrow);
            is_minus_one &= ct_is_zero_l(base_pow_val ^ minus_one_val);
        }

        // Evaluate the Miller-Rabin conditions for the
        // exp_bit_pos == p_minus_one_first_set_bit_pos case.
        let bit_pos_eq_first_set = ct_eq_usize_usize(exp_bit_pos, p_minus_one_first_set_bit_pos);
        is_probable_prime |= bit_pos_eq_first_set.select(0, is_one | is_minus_one);

        // Evaluate the Miller-Rabin conditions for the case that
        // 0 < exp_bit_pos < p_minus_one_first_set_bit_pos.
        //
        // Implementations found in literature often break out early with a result of
        // "composite" if the intermediate power of the base, which is a square
        // by the range of exp_bit_pos, equals one at this point. This is an
        // optimization though, because once the value is one, squaring cannot
        // ever yield -1 (in fact, nothing but 1) again. This is a basic property of
        // (pretty much by definition) a unit in a ring, independent of whether p is a
        // prime or not. This optimization is effective only if p - 1 has a
        // significant number of trailing zeroes, which, for random prime
        // candidates is highly unlikely. So don't bother and keep the algorithm
        // constant-time even for non-primes.
        let bit_pos_below_first_set = ct_lt_usize_usize(exp_bit_pos, p_minus_one_first_set_bit_pos);
        is_probable_prime |= bit_pos_below_first_set.select(0, is_minus_one);
    }

    Ok(LimbChoice::from(is_probable_prime))
}

#[cfg(test)]
fn test_ct_prime_test_miller_rabin_mp<
    BT: MpMutUIntSlice,
    PT: MpMutUIntSlice,
    RXT: MpMutUIntSlice,
>() {
    extern crate alloc;
    use super::add::ct_sub_mp_l;
    use super::limbs::ct_mp_nlimbs;
    use super::mul::ct_mul_trunc_mp_l;
    use alloc::vec;

    fn is_probable_prime<BT: MpMutUIntSlice, PT: MpUIntCommon, RXT: MpMutUIntSlice>(
        base: &BT,
        p: &PT,
    ) -> bool {
        use super::div::{CtMpDivisor, ct_mod_mp_mp};
        use super::montgomery::{
            ct_montgomery_radix2_mod_n_mp, ct_montgomery_redc_mp, ct_to_montgomery_form_mp,
        };

        let mut base_mod_p = tst_mk_mp_backing_vec!(BT, base.len());
        let mut base_mod_p = BT::from_slice(&mut base_mod_p).unwrap();
        base_mod_p.copy_from(base);
        ct_mod_mp_mp(None, &mut base_mod_p, &CtMpDivisor::new(p, None).unwrap());

        let mut mg_radix2_mod_p = tst_mk_mp_backing_vec!(RXT, p.len());
        let mut mg_radix2_mod_p = RXT::from_slice(&mut mg_radix2_mod_p).unwrap();
        ct_montgomery_radix2_mod_n_mp(&mut mg_radix2_mod_p, p).unwrap();

        let mut mg_base_mod_p = tst_mk_mp_backing_vec!(RXT, p.len());
        let mut mg_base_mod_p = RXT::from_slice(&mut mg_base_mod_p).unwrap();
        let neg_p0_inv_mod_l = ct_montgomery_neg_n0_inv_mod_l_mp(p).unwrap();
        ct_to_montgomery_form_mp(
            &mut mg_base_mod_p,
            &base_mod_p,
            p,
            neg_p0_inv_mod_l,
            &mg_radix2_mod_p,
        )
        .unwrap();

        let mut mg_radix_mod_p = mg_radix2_mod_p;
        ct_montgomery_redc_mp(&mut mg_radix_mod_p, p, neg_p0_inv_mod_l).unwrap();

        let mut scratch0 =
            vec![0 as LimbType; MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(p.len())];
        let mut scratch1 =
            vec![0 as LimbType; MpMutNativeEndianUIntLimbsSlice::nlimbs_for_len(p.len())];
        ct_prime_test_miller_rabin_mp(
            &mg_base_mod_p,
            p,
            &mg_radix_mod_p,
            [scratch0.as_mut_slice(), scratch1.as_mut_slice()],
        )
        .unwrap()
        .unwrap()
            != 0
    }

    // p = 2^255 - 19.
    let mut p = tst_mk_mp_backing_vec!(PT, 256 / 8);
    let mut p = PT::from_slice(&mut p).unwrap();
    p.set_bit_to(255, true);
    ct_sub_mp_l(&mut p, 19);
    for i in 1..256 {
        let mut base = tst_mk_mp_backing_vec!(BT, 1);
        let mut base = BT::from_slice(&mut base).unwrap();
        base.store_l(0, i);
        assert!(is_probable_prime::<_, _, RXT>(&base, &p));
    }
    for i in 1..256 {
        let mut base = tst_mk_mp_backing_vec!(BT, p.len());
        let mut base = BT::from_slice(&mut base).unwrap();
        for j in 0..ct_mp_nlimbs(base.len()) {
            base.store_l(j, i);
        }
        base.store_l(0, i);
        assert!(is_probable_prime::<_, _, RXT>(&base, &p));
    }

    // c = 251 * (2^255 - 19)
    let mut c = tst_mk_mp_backing_vec!(PT, (256 + 8) / 8);
    let mut c = PT::from_slice(&mut c).unwrap();
    c.copy_from(&p);
    ct_mul_trunc_mp_l(&mut c, (256 + 8) / 8, 251);
    for i in 2..256 {
        let mut base = tst_mk_mp_backing_vec!(BT, 1);
        let mut base = BT::from_slice(&mut base).unwrap();
        base.store_l(0, i);
        assert!(!is_probable_prime::<_, _, RXT>(&base, &c));
    }
    for i in 1..256 {
        let mut base = tst_mk_mp_backing_vec!(BT, p.len());
        let mut base = BT::from_slice(&mut base).unwrap();
        for j in 0..ct_mp_nlimbs(base.len()) {
            base.store_l(j, i);
        }
        base.store_l(0, i);
        assert!(!is_probable_prime::<_, _, RXT>(&base, &c));
    }
}

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

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

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

/// The first few primes, in increasing order.
///
/// The primorials used by the wheel sieve, c.f. [`PrimeWheelSieve`], are
/// products of a prefix of these.
const FIRST_PRIMES: [u8; 15] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47];

/// Find the first element greater than or equal to a given value in constant
/// time.
///
/// Returns the index of the first element in the sorted slice `v` that is
/// greater than or equal to `val`, by means of a binary search. Both halves
/// are kept equal in length in each step, so that the execution time depends
/// only on `v`'s length and not on `val`. See [`cnst_find_first_ge_than()`] for
/// the non-constant-time variant usable in const contexts.
///
/// # Arguments
///
/// * `v` - The slice to search, sorted in increasing order. Its last element
///   must be greater than or equal to `val`.
/// * `val` - The value to search for.
fn ct_find_first_ge_than(v: &[u8], val: u8) -> usize {
    debug_assert!(*v.iter().last().unwrap() >= val);

    // Binary search. Be careful to always keep the left and right
    // halves equal in length for constant-time.
    let mut lb = 0;
    let mut ub = v.len();
    let mut r = ub - lb;
    while r > 1 {
        // m points past the mid element if r is odd, inbetween the two
        // mid elements otherwise.
        // Don't use usize::div_ceil(), which is implemented with a conditional.
        #[allow(clippy::manual_div_ceil)]
        let m = lb + (r + 1) / 2;
        let pivot_is_lt = ct_lt_l_l(v[m - 1] as LimbType, val as LimbType);
        // Keep the halves balanced for constant-time. If r is odd and
        // the right half is taken, extend it at its front by the pivot
        // element, even though it's not needed for the correctness of
        // the search.
        let r_odd = r & 1;
        debug_assert_eq!(ub - m + r_odd, m - lb);
        lb = pivot_is_lt.select_usize(lb, m - r_odd);
        ub = pivot_is_lt.select_usize(m, ub);
        r = ub - lb;
    }

    lb
}

#[test]
fn test_ct_find_first_ge_than() {
    for val in FIRST_PRIMES.iter() {
        let i = ct_find_first_ge_than(&FIRST_PRIMES, *val);
        assert_eq!(FIRST_PRIMES[i], *val);
    }

    for (i, first_prime) in FIRST_PRIMES.iter().enumerate().take(FIRST_PRIMES.len() - 1) {
        let val = first_prime + 1;
        let j = ct_find_first_ge_than(&FIRST_PRIMES, val);
        assert_eq!(i + 1, j);
    }
}

/// Find the first element greater than or equal to a given value in const
/// contexts.
///
/// Like [`ct_find_first_ge_than()`], but usable in const contexts and not
/// constant-time.
///
/// # Arguments
///
/// * `v` - The slice to search, sorted in increasing order. Must be non-empty
///   and its last element must be greater than or equal to `val`.
/// * `val` - The value to search for.
const fn cnst_find_first_ge_than(v: &[u8], val: u8) -> usize {
    debug_assert!(!v.is_empty());
    debug_assert!(v[v.len() - 1] >= val);

    // Binary search.
    let mut lb = 0;
    let mut ub = v.len();
    let mut r = ub - lb;
    while r > 1 {
        // m points past the mid element if r is odd, inbetween the two
        // mid elements otherwise.
        let m = lb + r.div_ceil(2);
        if v[m - 1] < val {
            lb = m;
        } else {
            ub = m;
        }
        r = ub - lb;
    }

    lb
}

#[test]
fn test_cnst_find_first_ge_than() {
    for val in FIRST_PRIMES.iter() {
        let i = cnst_find_first_ge_than(&FIRST_PRIMES, *val);
        assert_eq!(FIRST_PRIMES[i], *val);
    }

    for (i, first_prime) in FIRST_PRIMES.iter().enumerate().take(FIRST_PRIMES.len() - 1) {
        let val = first_prime + 1;
        let j = cnst_find_first_ge_than(&FIRST_PRIMES, val);
        assert_eq!(i + 1, j);
    }
}

/// Determine the maximum number of primorial factors not exceeding a bound.
///
/// Returns the maximum number *n* of [`FIRST_PRIMES`] such that their
/// product, the primorial, is at most `max`, capped to the number of
/// [`FIRST_PRIMES`] available.
///
/// # Arguments
///
/// * `max` - The upper bound on the primorial.
const fn first_primes_primorial_nfactors_for_max(max: LimbType) -> usize {
    let mut prod: LimbType = 1;
    let mut n = 0;
    while n < FIRST_PRIMES.len() {
        let next_p = FIRST_PRIMES[n] as LimbType;
        if max / next_p < prod {
            break;
        }
        prod *= next_p;
        n += 1;
    }

    n
}

/// Compute the primorial of the first few primes.
///
/// Returns the product of the first `n` [`FIRST_PRIMES`].
///
/// # Arguments
///
/// * `n` - The number of factors. The product must fit a [`LimbType`], c.f.
///   [`first_primes_primorial_nfactors_for_max()`].
const fn first_primes_primorial(n: usize) -> LimbType {
    debug_assert!(n <= first_primes_primorial_nfactors_for_max(LimbType::MAX));
    let mut prod: LimbType = 1;
    let mut i = 0;
    while i < n {
        prod *= FIRST_PRIMES[i] as LimbType;
        i += 1;
    }
    prod
}

/// Level 0 wheel sieve.
///
/// Enumerates the numbers not divisible by any of the first
/// [`PRIMORIAL_NFACTORS`](Self::PRIMORIAL_NFACTORS) primes, i.e. coprime to
/// the level 0 primorial *2 · 3 · 5 = 30*, by producing the deltas between
/// subsequent ones. Serves as a generator of offset candidates for the level 1
/// wheel, [`PrimeWheelSieveLvl1`]. The offsets of the wheel within a turn are
/// precomputed in [`OFFSETS`](Self::OFFSETS).
struct PrimeWheelSieveLvl0 {
    /// The index into [`OFFSETS`](Self::OFFSETS) of the offset to produce
    /// next.
    next_index: usize,
    /// The position within the current turn of the wheel produced last, or the
    /// initial lower bound modulo the primorial.
    last_offset: u8,
    /// Additional increment to apply for the next delta if the wheel has just
    /// wrapped around to its next turn, zero otherwise.
    next_offset_delta: u8,
}

impl PrimeWheelSieveLvl0 {
    /// The number of [`FIRST_PRIMES`] making up the level 0 primorial.
    ///
    /// The primorial of the first four primes would still fit an [`u8`].
    /// However, by sticking to the first three factors only, the list of
    /// wheel offsets is _much_ smaller.
    const PRIMORIAL_NFACTORS: usize = 3;
    /// The level 0 primorial, *2 · 3 · 5 = 30*.
    const PRIMORIAL: u8 = first_primes_primorial(Self::PRIMORIAL_NFACTORS) as u8;

    /// Compile-time assertion that all of the wheel's offsets are primes.
    ///
    /// This holds if the square of the next prime not among the primorial's
    /// factors exceeds the primorial, which allows for looking the number of
    /// offsets up in [`FIRST_PRIMES`], c.f. [`NOFFSETS`](Self::NOFFSETS).
    const fn assert_all_offsets_are_first_primes() {
        // By choice of parameters, the next prime not a factor of the level 0 primorial
        // squared is >= the primorial. It follows that all numbers less than
        // the primorial have a square root < that next prime. It follows in
        // turn that any composite number < the primorial has at
        // least one factor < that next prime, i.e. some common factor with the
        // primorial. For the wheel sieve in general, the set of offsets is
        // defined to be the set of numbers < the primorial, that have no factor
        // in common with it. By the preceeding, these cannot be composite, i.e.
        // must be prime.
        assert!(
            (FIRST_PRIMES[Self::PRIMORIAL_NFACTORS] as u64
                * FIRST_PRIMES[Self::PRIMORIAL_NFACTORS] as u64)
                > Self::PRIMORIAL as u64
        );
    }

    /// The number of offsets within a turn of the wheel.
    const NOFFSETS: usize = {
        Self::assert_all_offsets_are_first_primes();
        cnst_find_first_ge_than(&FIRST_PRIMES, Self::PRIMORIAL) - Self::PRIMORIAL_NFACTORS + 1
    };

    /// The offsets within a turn of the wheel, i.e. the numbers less than the
    /// primorial and coprime to it, in increasing order.
    const OFFSETS: [u8; Self::NOFFSETS] = [1, 7, 11, 13, 17, 19, 23, 29];
    /// The largest offset within a turn of the wheel.
    const MAX_OFFSET: u8 = Self::OFFSETS[Self::NOFFSETS - 1];

    /// The distance from the largest offset to the end of the wheel's turn.
    const fn offset_delta_at_wrap() -> u8 {
        Self::PRIMORIAL - Self::MAX_OFFSET
    }

    /// Start the wheel at a lower bound.
    ///
    /// The first delta produced advances `lower_bound` to the first number
    /// greater than or equal to it that is coprime to the primorial.
    ///
    /// Uses the CPU's division instruction for reducing `lower_bound` modulo
    /// the primorial, which might not run in constant time on all platforms.
    ///
    /// # Arguments
    ///
    /// * `lower_bound` - The lower bound to start from.
    fn start_geq_than(lower_bound: LimbType) -> Self {
        // Take the lower bound modulo the primorial to get a lower bound
        // on the first offset.
        let offset_lower_bound = (lower_bound % Self::PRIMORIAL as LimbType) as u8;

        // Lookup the first wheel offset >= the lower bound.
        // The wheel's last offset is always equal to the primorial minus 1,
        // so the search is well-defined.
        debug_assert_eq!(Self::OFFSETS[Self::NOFFSETS - 1], Self::PRIMORIAL - 1);
        let offset_index = ct_find_first_ge_than(&Self::OFFSETS, offset_lower_bound);
        // By initializing last_offset to offset_lower_bound, the first produced
        // delta will advance the input lower_bound to the first offset.
        Self {
            next_index: offset_index,
            last_offset: offset_lower_bound,
            next_offset_delta: 0,
        }
    }

    /// Produce the delta to the next number coprime to the primorial.
    ///
    /// Runs in constant time.
    fn produce_next_delta(&mut self) -> u8 {
        let next_offset = Self::OFFSETS[self.next_index];
        let delta = next_offset + self.next_offset_delta - self.last_offset;
        self.advance();
        delta
    }

    /// Advance to the next offset, wrapping around to the next turn of the
    /// wheel as appropriate.
    ///
    /// Runs in constant time.
    fn advance(&mut self) {
        let last_index = self.next_index;
        self.next_index += 1;
        let wrapped = ct_eq_usize_usize(self.next_index, Self::NOFFSETS);
        self.next_index = wrapped.select_usize(self.next_index, 0);
        self.last_offset = wrapped.select(Self::OFFSETS[last_index] as LimbType, 0) as u8;
        self.next_offset_delta = wrapped.select(0, Self::offset_delta_at_wrap() as LimbType) as u8;
    }
}

#[test]
fn test_prime_wheel_sieve_lvl0() {
    fn advance_candidate(candidate: &mut LimbType, wheel: &mut PrimeWheelSieveLvl0) {
        let candidate_is_zero = *candidate == 0;
        let delta = wheel.produce_next_delta();
        for j in 0..delta {
            *candidate += 1;
            if candidate_is_zero && j == 0 {
                continue;
            }
            let mut is_not_coprime = false;
            for factor in &FIRST_PRIMES[..PrimeWheelSieveLvl0::PRIMORIAL_NFACTORS] {
                let factor = *factor as LimbType;
                let rem = *candidate % factor;
                if rem == 0 {
                    is_not_coprime = true;
                    break;
                }
            }
            if j + 1 != delta {
                assert!(is_not_coprime);
            } else {
                assert!(!is_not_coprime);
            }
        }
    }

    let mut candidate = 0;
    let mut wheel = PrimeWheelSieveLvl0::start_geq_than(candidate);
    advance_candidate(&mut candidate, &mut wheel);
    assert_eq!(candidate, 1);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }

    let mut candidate = 1;
    let mut wheel = PrimeWheelSieveLvl0::start_geq_than(candidate);
    advance_candidate(&mut candidate, &mut wheel);
    assert_eq!(candidate, 1);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }

    let mut candidate = 2;
    let mut wheel = PrimeWheelSieveLvl0::start_geq_than(candidate);
    advance_candidate(&mut candidate, &mut wheel);
    assert_eq!(candidate, 7);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }

    let mut candidate = PrimeWheelSieveLvl0::PRIMORIAL as LimbType - 2;
    let mut wheel = PrimeWheelSieveLvl0::start_geq_than(candidate);
    advance_candidate(&mut candidate, &mut wheel);
    assert_eq!(candidate, PrimeWheelSieveLvl0::PRIMORIAL as LimbType - 1);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }

    let mut candidate = PrimeWheelSieveLvl0::PRIMORIAL as LimbType - 1;
    let mut wheel = PrimeWheelSieveLvl0::start_geq_than(candidate);
    advance_candidate(&mut candidate, &mut wheel);
    assert_eq!(candidate, PrimeWheelSieveLvl0::PRIMORIAL as LimbType - 1);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }

    let mut candidate = PrimeWheelSieveLvl0::PRIMORIAL as LimbType;
    let mut wheel = PrimeWheelSieveLvl0::start_geq_than(candidate);
    advance_candidate(&mut candidate, &mut wheel);
    assert_eq!(candidate, PrimeWheelSieveLvl0::PRIMORIAL as LimbType + 1);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }
}

/// Wheel sieve for prime candidate generation.
///
/// Enumerates prime candidates greater than or equal to a given lower bound,
/// which aren't divisible by any of the first few primes, i.e. which are
/// coprime to a [`LimbType`]-sized primorial, by producing the deltas between
/// subsequent ones. For a [`LimbType`] of 64 bits, the primorial's factors are
/// the primes from 2 to 47, for one of 32 bits those from 2 to 23. Note that
/// these small primes themselves are never produced. Exported as
/// [`PrimeWheelSieve`].
///
/// To be used as follows: initialize a candidate with the lower bound, start
/// the wheel with [`start_geq_than()`](Self::start_geq_than) and add each
/// delta produced by [`produce_next_delta()`](Self::produce_next_delta) to the
/// candidate to obtain the next one. The first delta advances the lower bound
/// to the first candidate.
///
/// The general method follows the "wheel sieve" approach: start out from a
/// multiple of a primorial as a base and add offsets not divisible by any of
/// its factors to it. The larger the primorial, the more effective the
/// preselection, because all potential candidates divisible by any of its
/// factors are filtered. On the other hand, the list of offsets grows about
/// linearly with the primorial's magnitude, and it's not feasible to store a
/// number of precomputed offsets proportional to e.g. *2<sup>64</sup>*. Hence a
/// two-level approach is taken: a level 0 wheel uses the very small primorial
/// *2 · 3 · 5 = 30* with precomputed offsets and serves the sole purpose of
/// generating offset candidates for this level 1 wheel, which filters those not
/// coprime to its larger primorial.
///
/// The wheel sieve is not constant-time: the number of level 0 offset
/// candidates to skip for producing the next delta, and thus the execution
/// time, depends on the current position in the wheel and hence on the lower
/// bound.
pub struct PrimeWheelSieveLvl1 {
    /// The level 0 wheel generating the offset candidates.
    lvl0_wheel: PrimeWheelSieveLvl0,
    /// The offset within the current turn of the wheel produced last, or the
    /// initial lower bound modulo the primorial.
    last_offset: LimbType,
    /// An additional increment to account for in the first delta produced, if
    /// the lower bound had to be advanced internally, zero otherwise.
    last_offset_delta: u8,
}

impl PrimeWheelSieveLvl1 {
    /// The number of [`FIRST_PRIMES`] making up the level 1 primorial, chosen
    /// as large as possible such that the primorial fits a [`LimbType`].
    const PRIMORIAL_NFACTORS: usize = first_primes_primorial_nfactors_for_max(LimbType::MAX);
    /// The level 1 primorial, i.e. the product of the first
    /// [`PRIMORIAL_NFACTORS`](Self::PRIMORIAL_NFACTORS) primes.
    const PRIMORIAL: LimbType = first_primes_primorial(Self::PRIMORIAL_NFACTORS);

    /// Start the wheel at a lower bound.
    ///
    /// The first delta produced advances `lower_bound` to the first number
    /// greater than or equal to it that is coprime to the primorial, but at
    /// least to the first such number greater than one.
    ///
    /// # Arguments
    ///
    /// * `lower_bound` - The lower bound to start from, a multiprecision
    ///   integer of any representation and length.
    pub fn start_geq_than<LBT: MpUIntCommon>(lower_bound: &LBT) -> Self {
        let offset_lower_bound = ct_mod_mp_l(
            lower_bound,
            &CtLDivisor::nonct_new(Self::PRIMORIAL).unwrap(),
        );
        // If the lower bound is <= 1, skip the 1 generated by the level 0.
        // It only got the potential for causing confusion or subtle corner cases.
        // It's important to account for that internal increment when producing
        // the first offset delta to the outside, so remember it in last_offset_delta.
        let lower_bound_leq_one = ct_leq_mp_l(lower_bound, 1);
        let last_offset_delta =
            lower_bound_leq_one.select(0, (2 as LimbType).wrapping_sub(offset_lower_bound));
        let offset_lower_bound = offset_lower_bound + last_offset_delta;
        // Remember that increment and adjust for it when emitting the first produced
        // offset delta.
        Self {
            lvl0_wheel: PrimeWheelSieveLvl0::start_geq_than(offset_lower_bound),
            last_offset: offset_lower_bound,
            last_offset_delta: last_offset_delta as u8,
        }
    }

    /// Produce the delta to the next prime candidate.
    ///
    /// Returns the delta to add to the current candidate, or, for the first
    /// invocation, to the lower bound, for obtaining the next number coprime to
    /// the primorial.
    ///
    /// Not constant-time, c.f. the [type level documentation](Self).
    pub fn produce_next_delta(&mut self) -> LimbType {
        let mut next_offset = self.last_offset;
        let mut next_offset_delta;
        let mut last_offset = self.last_offset;
        loop {
            let lvl0_delta = self.lvl0_wheel.produce_next_delta();
            // The addition does not wrap: next_offset < Self::PRIMORIAL and Self::PRIMORIAL
            // has been chosen such that there's enough room until LimbType::MAX
            // to accomodate the largest delta value the level 0 wheel could
            // ever produce.
            next_offset += lvl0_delta as LimbType;
            let wrapped = ct_geq_l_l(next_offset, Self::PRIMORIAL);
            next_offset = wrapped.select(next_offset, next_offset.wrapping_sub(Self::PRIMORIAL));
            last_offset = wrapped.select(last_offset, 0);
            next_offset_delta = wrapped.select(0, Self::PRIMORIAL - self.last_offset);

            let mut is_not_coprime = LimbChoice::from(0);
            for factor in FIRST_PRIMES
                .iter()
                .take(Self::PRIMORIAL_NFACTORS)
                .skip(PrimeWheelSieveLvl0::PRIMORIAL_NFACTORS)
            {
                let rem = next_offset % *factor as LimbType;
                is_not_coprime |= ct_eq_l_l(rem, 0);
            }
            if is_not_coprime.unwrap() != 0 {
                continue;
            }
            break;
        }

        let last_offset_delta = self.last_offset_delta;
        self.last_offset_delta = 0;
        self.last_offset = next_offset;
        next_offset + next_offset_delta - (last_offset - last_offset_delta as LimbType)
    }
}

#[test]
fn test_prime_wheel_sieve_lvl1() {
    use super::add::ct_add_mp_l;
    use super::cmp::ct_is_zero_mp;

    fn advance_candidate(
        candidate: &mut MpMutNativeEndianUIntLimbsSlice,
        wheel: &mut PrimeWheelSieveLvl1,
    ) {
        let candidate_is_zero = ct_is_zero_mp(candidate).unwrap() != 0;
        let delta = wheel.produce_next_delta();
        for j in 0..delta {
            ct_add_mp_l(candidate, 1);
            if candidate_is_zero && j == 0 {
                continue;
            }
            let mut is_not_coprime = false;
            for factor in &FIRST_PRIMES[..PrimeWheelSieveLvl1::PRIMORIAL_NFACTORS] {
                let factor = *factor as LimbType;
                let factor = CtLDivisor::nonct_new(factor).unwrap();
                let rem = ct_mod_mp_l(candidate, &factor);
                if rem == 0 {
                    is_not_coprime = true;
                    break;
                }
            }
            if j + 1 != delta {
                assert!(is_not_coprime);
            } else {
                assert!(!is_not_coprime);
            }
        }
    }

    let mut candidate = [0 as LimbType; 2];
    let mut candidate = MpMutNativeEndianUIntLimbsSlice::from_limbs(&mut candidate);
    let mut wheel = PrimeWheelSieveLvl1::start_geq_than(&candidate);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }

    let mut candidate = [0 as LimbType; 2];
    let mut candidate = MpMutNativeEndianUIntLimbsSlice::from_limbs(&mut candidate);
    candidate.store_l(0, 1);
    let mut wheel = PrimeWheelSieveLvl1::start_geq_than(&candidate);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }

    let mut candidate = [0 as LimbType; 2];
    let mut candidate = MpMutNativeEndianUIntLimbsSlice::from_limbs(&mut candidate);
    candidate.store_l(0, 2);
    let mut wheel = PrimeWheelSieveLvl1::start_geq_than(&candidate);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }

    let mut candidate = [0 as LimbType; 2];
    let mut candidate = MpMutNativeEndianUIntLimbsSlice::from_limbs(&mut candidate);
    candidate.store_l(0, PrimeWheelSieveLvl1::PRIMORIAL - 2);
    let mut wheel = PrimeWheelSieveLvl1::start_geq_than(&candidate);
    advance_candidate(&mut candidate, &mut wheel);
    assert_eq!(candidate.load_l(0), PrimeWheelSieveLvl1::PRIMORIAL - 1);
    for _i in 0..1024 {
        advance_candidate(&mut candidate, &mut wheel);
    }
}

pub use PrimeWheelSieveLvl1 as PrimeWheelSieve;