//! A macro to define efficient and constant-time arithmetic for finite fields Fp with
//! p = 3 mod 4 using Montgomery multiplication.
//!
//! # Traits
//!
//! This macro defines a finite field type for GF(p) and also implements the trait Fp
//! by re-exporting all necessary functions.
//!
//! # Authorship and History
//!
//! The majority of this code has been adapted from code written by Thomas Pornin
//! from collaboration in previous projects and several methods which appear in other
//! macros in the cryptographic research library crrl <https://github.com/pornin/crrl>
//!
//! This code has also been used in a handful of isogeny-based cryptography research
//! projects before being rewritten for this crate, including:
//! - <https://github.com/ThetaIsogenies/two-isogenies>
//! - <https://github.com/GiacomoPope/cubical-pairings>
//! - <https://github.com/GiacomoPope/ThetaCGL>

/// A macro to define the finite field Fp, all functions are designed to run in
/// constant time. Assumes that the characteristic is p = 3 mod 4.
///
/// Macro expectations:
/// - A typename for the finite field generated
/// - An array of `N` words which represent the finite field characteristic
///   in base 2^64
#[macro_export]
macro_rules! define_fp_core {
    (
        typename = $typename:ident,
        modulus = $modulus:expr,
    ) => {
        /// A finite field element. Contents are opaque.
        /// All functions are constant-time.
        ///
        /// A field element x is encoded into bytes by using the unsigned
        /// little-endian convention over the unique representant of x in the
        /// [0..(p-1)] range. There is no sign bit.
        #[derive(Clone, Copy, Debug)]
        pub struct $typename([u64; $typename::N]);

        // For very large moduli (such as 4000+ bits) some of the const functions can
        // take a long time. This is solved in two ways.
        //
        // 1. Allow the user to obtain R, R2, TDEC etc as macro parameters
        // 2. Be patient
        //
        // I have gone with option two as it makes the macro input cleaner and most
        // of the time smaller moduli are used and the compile time is not an issue.
        #[allow(long_running_const_eval)]
        impl $typename {
            // IMPLEMENTATION NOTES
            // --------------------
            //
            // Modulus is p. Each element is represented over N limbs, in base
            // 2^64.
            //
            // Let R = 2^(64*N) mod p. A field element x is represented by
            // the integer x*R mod p, in the [0..(p-1)] range. The limbs are in
            // little-endian order (limb 0 is the least significant).
            //
            // Multiplications use Montgomery multiplication: given x and y,
            // the value (x*y)/R mod p is computed. Since our values are in
            // Montgomery representation, what is computed is really
            // (x*R)*(y*R)/R = x*y*R mod p, which is the correct product. The
            // decoding and encoding functions apply the required conversions; the
            // use of Montgomery representation is not visible to other code using
            // this type.

            // Number of words and bit length of the field characteristic
            pub const N: usize = Self::top_word_index() + 1;
            pub const BIT_LENGTH: usize = Self::mod_bitlen();
            pub const MODULUS: [u64; Self::N] = $modulus;

            // Constants used for internal arithmetic
            const P0I: u64 = Self::ninv64(Self::MODULUS[0]);
            const R: Self = Self::pow2mod(Self::N * 64);
            const R2: Self = Self::pow2mod(Self::N * 128);
            const P1: u64 = Self::top_u32();
            const P1DIV_M: u64 =
                1 + ((((((1u64 << 32) - Self::P1) as u128) << 64) / (Self::P1 as u128)) as u64);
            // const NUM1: usize = (2 * Self::BIT_LENGTH - 34) / 31;
            // const NUM2: usize = 2 * Self::BIT_LENGTH - 31 * Self::NUM1 - 2;

            // When the modulus has 32 bits or less, we can set NUM1 = 0
            const NUM1: usize = {
                if Self::BIT_LENGTH <= 32 {
                    0
                } else {
                    (2 * Self::BIT_LENGTH - 34) / 31 // original formula, safe now
                }
            };

            // NUM2 covers the remaining bit-length after NUM1 outer iterations.
            const NUM2: usize = {
                let remaining = (2 * Self::BIT_LENGTH).saturating_sub(31 * Self::NUM1);
                if remaining >= 2 { remaining - 2 } else { 0 }
            };

            // Constant for correction after the binary GCD loop
            const TFIXDIV: Self = Self::const_mmul(
                Self::const_mmul(Self::pow2mod(Self::NUM1 * 33 + 64 - Self::NUM2), Self::R2),
                Self::R2,
            );
            const SQRT_EXP: [u64; Self::N] = Self::const_sqrt_exp();
            const FOURTH_ROOT_EXP: [u64; Self::N] = Self::const_fourth_root_exp();
            pub const SUM_OF_PRODUCTS_ADDITIONAL_SUB: bool = Self::sum_of_products_check();

            /// Encoding length of a field element (in bytes). All elements
            /// always encode into exactly that many bytes. Encoding is
            /// canonical: a given field element has a unique valid encoding,
            /// and the decoding process verifies that this specific encoding
            /// was used.
            pub const ENCODED_LENGTH: usize = (Self::BIT_LENGTH + 7) >> 3;

            // Multiplier for decode_reduce().
            const CLEN: usize = if Self::N == 1 {
                Self::ENCODED_LENGTH - 1
            } else {
                8 * (Self::N - 1)
            };
            const TDEC: Self = Self::pow2mod(Self::N * 64 + 8 * Self::CLEN);

            // Predefined constants used externally
            pub const ZERO: Self = Self([0u64; Self::N]);
            pub const ONE: Self = Self::R;
            pub const TWO: Self = Self::const_small(2);
            pub const THREE: Self = Self::const_small(3);
            pub const FOUR: Self = Self::const_small(4);
            pub const MINUS_ONE: Self = Self::const_neg(Self::R);

            pub const fn new(input: [u64; Self::N]) -> Self {
                return Self(input);
            }

            /// Return 0xFFFFFFFF if this value is zero, or 0x00000000 otherwise.
            #[inline]
            pub fn is_zero(self) -> u32 {
                let mut x = self.0[0];
                for i in 1..Self::N {
                    x |= self.0[i];
                }
                (!$crate::utils64::sgnw(x | x.wrapping_neg())) as u32
            }

            /// Return 0xFFFFFFFF if this value is equal to rhs, or 0x00000000
            /// otherwise.
            #[inline(always)]
            pub fn equals(self, rhs: &Self) -> u32 {
                let mut r = 0u64;
                for i in 0..Self::N {
                    r |= self.0[i] ^ rhs.0[i];
                }
                (((r | r.wrapping_neg()) >> 63) as u32).wrapping_sub(1)
            }

            /// Add `rhs` to this value.
            #[inline]
            fn set_add(&mut self, rhs: &Self) {
                // raw addition.
                let mut cc1 = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::addcarry_u64(self.0[i], rhs.0[i], cc1);
                    self.0[i] = d;
                    cc1 = ee;
                }

                // subtract modulus.
                let mut cc2 = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::subborrow_u64(self.0[i], Self::MODULUS[i], cc2);
                    self.0[i] = d;
                    cc2 = ee;
                }

                // add back modulus if the result was negative, i.e. cc1 - cc2 < 0.
                let mm = (cc1 as u64).wrapping_sub(cc2 as u64);
                let mut cc = 0;
                for i in 0..Self::N {
                    let (d, ee) =
                        $crate::utils64::addcarry_u64(self.0[i], mm & Self::MODULUS[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }
            }

            /// Subtract `rhs` from this value.
            #[inline]
            fn set_sub(&mut self, rhs: &Self) {
                // raw subtraction
                let mut cc = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::subborrow_u64(self.0[i], rhs.0[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }

                // add back modulus if the result was negative
                let mm = (cc as u64).wrapping_neg();
                cc = 0;
                for i in 0..Self::N {
                    let (d, ee) =
                        $crate::utils64::addcarry_u64(self.0[i], mm & Self::MODULUS[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }
            }

            /// Negate this value.
            #[inline]
            pub fn set_neg(&mut self) {
                // subtract from zero
                let mut cc = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::subborrow_u64(0, self.0[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }

                // add back the modulus if needed (i.e. if input was non-zero)
                let mm = (cc as u64).wrapping_neg();
                cc = 0;
                for i in 0..Self::N {
                    let (d, ee) =
                        $crate::utils64::addcarry_u64(self.0[i], mm & Self::MODULUS[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }
            }

            // Perform Montgomery reduction (division by R) on this value.
            // Internal note: if self has proper contents (value less than p), then
            // this necessarily yields a properly reduced value. If self is not
            // properly reduced, then the output is in [0..p] inclusive.
            #[inline]
            fn set_montyred(&mut self) {
                for _ in 0..Self::N {
                    let f = self.0[0].wrapping_mul(Self::P0I);
                    let (_, mut cc) = $crate::utils64::umull_add(f, Self::MODULUS[0], self.0[0]);
                    for i in 1..Self::N {
                        let (d, hi) =
                            $crate::utils64::umull_add2(f, Self::MODULUS[i], self.0[i], cc);
                        self.0[i - 1] = d;
                        cc = hi;
                    }
                    self.0[Self::N - 1] = cc;
                }
            }

            /// Multiply this value by `rhs`, optimised for when N is "small"
            #[inline]
            fn set_mul_small_word_len(&mut self, rhs: &Self) {
                let mut t = Self::ZERO;

                let mut cch: u8 = 0;
                for i in 0..Self::N {
                    let (lo, mut cc1) = $crate::utils64::umull_add(rhs.0[i], self.0[0], t.0[0]);
                    t.0[0] = lo;
                    for j in 1..Self::N {
                        let (d, hi1) =
                            $crate::utils64::umull_add2(rhs.0[i], self.0[j], t.0[j], cc1);
                        cc1 = hi1;
                        t.0[j] = d;
                    }

                    let q = t.0[0].wrapping_mul(Self::P0I);

                    let (_, mut cc2) = $crate::utils64::umull_add(q, Self::MODULUS[0], t.0[0]);
                    for j in 1..Self::N {
                        let (d, hi2) =
                            $crate::utils64::umull_add2(q, Self::MODULUS[j], t.0[j], cc2);
                        cc2 = hi2;
                        t.0[j - 1] = d;
                    }
                    let (d, ee) = $crate::utils64::addcarry_u64(cc1, cc2, cch);
                    t.0[Self::N - 1] = d;
                    cch = ee;
                }

                // final reduction: subtract modulus if necessary
                let mut cc = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::subborrow_u64(t.0[i], Self::MODULUS[i], cc);
                    t.0[i] = d;
                    cc = ee;
                }
                let mask = (cch as u64).wrapping_sub(cc as u64);
                cc = 0;
                for i in 0..Self::N {
                    let (d, ee) =
                        $crate::utils64::addcarry_u64(t.0[i], mask & Self::MODULUS[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }
            }

            /// Multiply this value by `rhs`, optimised for when N is "large"
            #[inline]
            fn set_mul_large_word_len(&mut self, rhs: &Self) {
                let mut t = Self::ZERO;

                // combined muls + reduction
                let mut cch = 0;
                for i in 0..Self::N {
                    let f = rhs.0[i];
                    let (lo, mut cc1) = $crate::utils64::umull_add(f, self.0[0], t.0[0]);
                    let g = lo.wrapping_mul(Self::P0I);
                    let (_, mut cc2) = $crate::utils64::umull_add(g, Self::MODULUS[0], lo);
                    for j in 1..Self::N {
                        let (d, hi1) = $crate::utils64::umull_add2(f, self.0[j], t.0[j], cc1);
                        cc1 = hi1;
                        let (d, hi2) = $crate::utils64::umull_add2(g, Self::MODULUS[j], d, cc2);
                        cc2 = hi2;
                        t.0[j - 1] = d;
                    }
                    let (d, ee) = $crate::utils64::addcarry_u64(cc1, cc2, cch);
                    t.0[Self::N - 1] = d;
                    cch = ee;
                }

                // final reduction: subtract modulus if necessary
                let mut cc = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::subborrow_u64(t.0[i], Self::MODULUS[i], cc);
                    t.0[i] = d;
                    cc = ee;
                }
                let mm = (cch as u64).wrapping_sub(cc as u64);
                cc = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::addcarry_u64(t.0[i], mm & Self::MODULUS[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }
            }

            /// Multiply this value by `rhs`.
            #[inline]
            fn set_mul(&mut self, rhs: &Self) {
                // TODO: what's the best bound here?
                if Self::N < 15 {
                    self.set_mul_small_word_len(rhs);
                } else {
                    self.set_mul_large_word_len(rhs);
                }
            }

            /// Replace this value with its square.
            #[inline]
            pub fn set_square(&mut self) {
                // FIXME: this turns out to be slower than set_mul() on x86_64
                // when N >= 23. This is probably due to the more complicated
                // loop bounds. Full unrolling helps, but can only be done at
                // the crate level:
                //   RUSTFLAGS="-C llvm-args=-unroll-threshold=1200"
                // This impacts all the code in the crate, and is thus
                // probably not a very good idea.

                // TODO handle this in a different way?
                if Self::N == 1 {
                    let r = *self;
                    self.set_mul(&r);
                    return;
                }

                // Compute the square over integers.
                let mut t = [0u64; Self::N << 1];

                // sum_{i<j} a_i*a_j*2^(64*(i+j)) < 2^(64*(2*N-1))
                // -> t[2*N-1] remains at zero
                let f = self.0[0];
                let (d, mut cc) = $crate::utils64::umull(f, self.0[1]);
                t[1] = d;
                for j in 2..Self::N {
                    let (d, hi) = $crate::utils64::umull_add(f, self.0[j], cc);
                    t[j] = d;
                    cc = hi;
                }
                t[Self::N] = cc;
                for i in 1..(Self::N - 1) {
                    let f = self.0[i];
                    let (d, mut cc) = $crate::utils64::umull_add(f, self.0[i + 1], t[(i << 1) + 1]);
                    t[(i << 1) + 1] = d;
                    for j in (i + 2)..Self::N {
                        let (d, hi) = $crate::utils64::umull_add2(f, self.0[j], t[i + j], cc);
                        t[i + j] = d;
                        cc = hi;
                    }
                    t[i + Self::N] = cc;
                }

                // Double the partial sum.
                // -> t contains sum_{i!=j} a_i*a_j*2^(64*(i+j))
                let mut cc = 0;
                for i in 1..((Self::N << 1) - 1) {
                    let w = t[i];
                    let ee = w >> 63;
                    t[i] = (w << 1) | cc;
                    cc = ee;
                }
                t[(Self::N << 1) - 1] = cc;

                // Add the squares a_i*a_i*w^(64*2*i).
                let mut cc = 0;
                for i in 0..Self::N {
                    let (lo, hi) = $crate::utils64::umull(self.0[i], self.0[i]);
                    let (d0, ee) = $crate::utils64::addcarry_u64(lo, t[i << 1], cc);
                    let (d1, ee) = $crate::utils64::addcarry_u64(hi, t[(i << 1) + 1], ee);
                    t[i << 1] = d0;
                    t[(i << 1) + 1] = d1;
                    cc = ee;
                }

                // Apply Montgomery reduction. We use the following facts:
                //  - upper half is necessarily less than p
                //  - set_montyred() accepts a full-limbs input and outputs a
                //    value of at most p
                //  - set_add() tolerates an input operand equal to p provided
                //    that the sum is less than 2*p
                self.0.copy_from_slice(&t[..Self::N]);
                self.set_montyred();
                let mut y = Self([0u64; Self::N]);
                y.0.copy_from_slice(&t[Self::N..]);
                self.set_add(&y);
            }

            /// Compute the square of this value.
            #[inline(always)]
            pub fn square(self) -> Self {
                let mut r = self;
                r.set_square();
                r
            }

            /// Square this value n times in place
            #[inline(always)]
            pub fn set_n_square(&mut self, n: u32) {
                for _ in 0..n {
                    self.set_square();
                }
            }

            /// Square this value n times
            #[inline(always)]
            pub fn n_square(self, n: u32) -> Self {
                let mut r = self;
                r.set_n_square(n);
                r
            }

            /// Halve this value.
            #[inline]
            pub fn set_half(&mut self) {
                let m = (self.0[0] & 1).wrapping_neg();
                let (mut dd, mut cc) =
                    $crate::utils64::addcarry_u64(self.0[0], m & Self::MODULUS[0], 0);
                dd >>= 1;
                for i in 1..Self::N {
                    let (x, ee) =
                        $crate::utils64::addcarry_u64(self.0[i], m & Self::MODULUS[i], cc);
                    cc = ee;
                    self.0[i - 1] = dd | (x << 63);
                    dd = x >> 1;
                }
                self.0[Self::N - 1] = dd | ((cc as u64) << 63);
            }

            /// Compute the half of this value.
            #[inline(always)]
            pub fn half(self) -> Self {
                let mut r = self;
                r.set_half();
                r
            }

            /// Double this value.
            #[inline]
            pub fn set_mul2(&mut self) {
                // Double (as an integer) and subtract the modulus.
                let mut cc = 0;
                let mut tb = 0;
                for i in 0..Self::N {
                    let w = self.0[i];
                    let t = (w << 1) | tb;
                    tb = w >> 63;
                    let (d, ee) = $crate::utils64::subborrow_u64(t, Self::MODULUS[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }

                // add back modulus if the result was negative
                let mm = tb.wrapping_sub(cc as u64);
                cc = 0;
                for i in 0..Self::N {
                    let (d, ee) =
                        $crate::utils64::addcarry_u64(self.0[i], mm & Self::MODULUS[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }
            }

            /// Compute the sum of this value with itself.
            #[inline(always)]
            pub fn mul2(self) -> Self {
                let mut r = self;
                r.set_mul2();
                r
            }

            /// Triple this value.
            #[inline]
            pub fn set_mul3(&mut self) {
                let r = self.mul2();
                *self += &r;
            }

            /// Compute the triple of this value.
            #[inline(always)]
            pub fn mul3(self) -> Self {
                let mut r = self;
                r.set_mul3();
                r
            }

            /// Quadruple this value.
            #[inline]
            pub fn set_mul4(&mut self) {
                self.set_mul2();
                self.set_mul2();
            }

            /// Compute the quadruple of this value.
            #[inline(always)]
            pub fn mul4(self) -> Self {
                let mut r = self;
                r.set_mul4();
                r
            }

            /// Multiply this value by 8
            #[inline]
            pub fn set_mul8(&mut self) {
                self.set_mul2();
                self.set_mul2();
                self.set_mul2();
            }

            /// Compute 8 times this value.
            #[inline(always)]
            pub fn mul8(self) -> Self {
                let mut r = self;
                r.set_mul8();
                r
            }

            /// Multiply this value by a small signed integer k.
            #[inline]
            pub fn set_mul_small(&mut self, k: i32) {
                // Special case: if the modulus fits in one word, just convert
                if Self::N == 1 {
                    *self *= Self::from(k as i64);
                    return;
                }

                // Get the absolute value of the multiplier (but remember the sign).
                let sk = (k >> 31) as u32;
                let ak = ((k as u32) ^ sk).wrapping_sub(sk);

                // Do the product over integers.
                let (d, mut hi) = $crate::utils64::umull(self.0[0], ak as u64);
                self.0[0] = d;
                for i in 1..Self::N {
                    let (d, ee) = $crate::utils64::umull_add(self.0[i], ak as u64, hi);
                    self.0[i] = d;
                    hi = ee;
                }

                // We write:
                //    p = p1*2^m + p0   (modulus)
                //    x = x1*2^m + x0   (unreduced product)
                // with:
                //    2^31 <= p1 < 2^32
                //    0 <= p0 < 2^m
                //    0 <= x0 < 2^m
                // Since the current value x is the product of the input (less
                // than p) by a multiplier of at most 2^31, we know that:
                //    0 <= x < p*2^31 < 2^(63+m)
                //    0 <= x1 < 2^63.
                // We compute:
                //    b = floor(x1/p1)
                // Analysis shows that floor(x/p) = b, b-1 or b+1.
                //
                // We thus obtain b, then increment it (unless b == p1), then
                // subtract b*p from x; we then add back p repeatedly until a
                // non-negative result is obtained. At most two conditional
                // additions are needed to achieve that result.
                //
                // Division by p1 can be done with the Granlund-Montgomery method:
                //    https://dl.acm.org/doi/10.1145/773473.178249
                // (LLVM usually applies that method, but may fail to do so if for
                // instance optimizing for code size on some platforms, thus it is
                // best to apply the method explicitly so that constant-time code
                // is more reliably achieved.)

                // Extract top word of x.
                let bl = Self::BIT_LENGTH & 63;
                let x1 = if bl == 0 {
                    (self.0[Self::N - 1] >> 32) | (hi << 32)
                } else if bl < 32 {
                    (self.0[Self::N - 1] << (32 - bl)) | (self.0[Self::N - 2] >> (32 + bl))
                } else if bl == 32 {
                    self.0[Self::N - 1]
                } else {
                    (hi << (96 - bl)) | (self.0[Self::N - 1] >> (bl - 32))
                };

                // Compute b = floor(x1/p1).
                let (_, t) = $crate::utils64::umull(x1, Self::P1DIV_M);
                let b = (x1.wrapping_sub(t) >> 1).wrapping_add(t) >> 31;

                // Add 1 to b, unless b == p1 (we cannot have b > p1).
                let b = b + (Self::P1.wrapping_sub(b) >> 63);

                // Subtract b*p from x.
                let mut cc1 = 0;
                let mut cc2 = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::umull_add(b, Self::MODULUS[i], cc1);
                    cc1 = ee;
                    let (d, ee) = $crate::utils64::subborrow_u64(self.0[i], d, cc2);
                    self.0[i] = d;
                    cc2 = ee;
                }
                let (mut hi, _) = $crate::utils64::subborrow_u64(hi, cc1, cc2);

                // Add p (at most twice) as long as the value is negative.
                for _ in 0..2 {
                    let m = $crate::utils64::sgnw(hi);
                    let mut cc = 0;
                    for i in 0..Self::N {
                        let (d, ee) =
                            $crate::utils64::addcarry_u64(self.0[i], m & Self::MODULUS[i], cc);
                        self.0[i] = d;
                        cc = ee;
                    }
                    hi = hi.wrapping_add(cc as u64);
                }

                // We computed self*|k|; we must adjust for the sign of k.
                self.set_cond_neg(sk);
            }

            /// Compute the product of this value by a small (unsigned) integer k.
            #[inline(always)]
            pub fn mul_small(self, k: i32) -> Self {
                let mut r = self;
                r.set_mul_small(k);
                r
            }

            /// Set this value to either a or b, depending on whether the control
            /// word ctl is 0x00000000 or 0xFFFFFFFF, respectively.
            /// The value of ctl MUST be either 0x00000000 or 0xFFFFFFFF.
            #[inline]
            pub fn set_select(&mut self, a: &Self, b: &Self, ctl: u32) {
                let c = (ctl as u64) | ((ctl as u64) << 32);
                for i in 0..Self::N {
                    let wa = a.0[i];
                    let wb = b.0[i];
                    self.0[i] = wa ^ (c & (wa ^ wb));
                }
            }

            /// Return a or b, if ctl is 0x00000000 or 0xFFFFFFFF, respectively.
            /// ctl MUST be either 0x00000000 or 0xFFFFFFFF.
            /// The value of ctl MUST be either 0x00000000 or 0xFFFFFFFF.
            #[inline]
            pub fn select(a: &Self, b: &Self, ctl: u32) -> Self {
                let mut r = Self::ZERO;
                r.set_select(a, b, ctl);
                r
            }

            /// Set this value to rhs if ctl is 0xFFFFFFFF; leave it unchanged if
            /// ctl is 0x00000000.
            /// The value of ctl MUST be either 0x00000000 or 0xFFFFFFFF.
            #[inline]
            pub fn set_cond(&mut self, rhs: &Self, ctl: u32) {
                let c = (ctl as u64) | ((ctl as u64) << 32);
                for i in 0..Self::N {
                    let wa = self.0[i];
                    let wb = rhs.0[i];
                    self.0[i] = wa ^ (c & (wa ^ wb));
                }
            }

            /// Negate this value if ctl is 0xFFFFFFFF; leave it unchanged if
            /// ctl is 0x00000000.
            /// The value of ctl MUST be either 0x00000000 or 0xFFFFFFFF.
            #[inline]
            pub fn set_cond_neg(&mut self, ctl: u32) {
                let v = -(self as &Self);
                self.set_cond(&v, ctl);
            }

            /// Exchange the values of a and b is ctl is 0xFFFFFFFF; leave both
            /// values unchanged if ctl is 0x00000000.
            /// The value of ctl MUST be either 0x00000000 or 0xFFFFFFFF.
            #[inline]
            pub fn cond_swap(a: &mut Self, b: &mut Self, ctl: u32) {
                let c = (ctl as u64) | ((ctl as u64) << 32);
                for i in 0..Self::N {
                    let wa = a.0[i];
                    let wb = b.0[i];
                    let wc = c & (wa ^ wb);
                    a.0[i] = wa ^ wc;
                    b.0[i] = wb ^ wc;
                }
            }

            // Set this value to (u*f+v*g)/2^64. Coefficients f
            // and g are provided as u64, but they are signed integers in the
            // [-2^62..+2^62] range.
            #[inline]
            fn set_montylin(&mut self, u: &Self, v: &Self, f: u64, g: u64) {
                // Make sure f and g are non-negative.
                let sf = $crate::utils64::sgnw(f);
                let f = (f ^ sf).wrapping_sub(sf);
                let tu = Self::select(u, &-u, sf as u32);
                let sg = $crate::utils64::sgnw(g);
                let g = (g ^ sg).wrapping_sub(sg);
                let tv = Self::select(v, &-v, sg as u32);

                let (d, mut cc) = $crate::utils64::umull_x2(tu.0[0], f, tv.0[0], g);
                self.0[0] = d;
                for i in 1..Self::N {
                    let (d, hi) = $crate::utils64::umull_x2_add(tu.0[i], f, tv.0[i], g, cc);
                    self.0[i] = d;
                    cc = hi;
                }
                let up = cc;

                // Montgomery reduction (one round)
                let k = self.0[0].wrapping_mul(Self::P0I);
                let (_, mut cc) = $crate::utils64::umull_add(k, Self::MODULUS[0], self.0[0]);
                for i in 1..Self::N {
                    let (d, hi) = $crate::utils64::umull_add2(k, Self::MODULUS[i], self.0[i], cc);
                    self.0[i - 1] = d;
                    cc = hi;
                }
                let (d, cc1) = $crate::utils64::addcarry_u64(up, cc, 0);
                self.0[Self::N - 1] = d;

                // |f| <= 2^62 and |g| <= 2^62, therefore |u*f + v*g| < p*2^63
                // We added less than p*2^64, and divided by 2^64, so the result
                // is less than 2*p and a single conditional subtraction is enough.
                let mut cc2 = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::subborrow_u64(self.0[i], Self::MODULUS[i], cc2);
                    self.0[i] = d;
                    cc2 = ee;
                }
                let mm = (cc1 as u64).wrapping_sub(cc2 as u64);
                let mut cc = 0;
                for i in 0..Self::N {
                    let (d, ee) =
                        $crate::utils64::addcarry_u64(self.0[i], mm & Self::MODULUS[i], cc);
                    self.0[i] = d;
                    cc = ee;
                }
            }

            #[inline(always)]
            fn montylin(a: &Self, b: &Self, f: u64, g: u64) -> Self {
                let mut r = Self::ZERO;
                r.set_montylin(a, b, f, g);
                r
            }

            // Set this value to abs((a*f + b*g)/2^31). Values a and b are
            // interpreted as plain integers (not modular). Coefficients f and
            // g are provided as u64 but they really are signed integers in the
            // [-2^31..+2^31] range (inclusive). The low 31 bits of a*f + b*g
            // are dropped (i.e. the division is assumed to be exact). The result
            // is assumed to fit in N limbs (extra high bits, if any, are
            // dropped). The absolute value of (a*f + b*g)/2^31 is computed.
            // Returned value is -1 (as a u64) if a*f + b*g was negative, 0
            // otherwise.
            #[inline]
            fn set_lindiv31abs(&mut self, a: &Self, b: &Self, f: u64, g: u64) -> u64 {
                // Replace f and g with abs(f) and abs(g), but remember the
                // original signs.
                let sf = $crate::utils64::sgnw(f);
                let f = (f ^ sf).wrapping_sub(sf);
                let sg = $crate::utils64::sgnw(g);
                let g = (g ^ sg).wrapping_sub(sg);

                // Compute a*f + b*g (upper word in 'up')
                let mut cc1 = 0;
                let mut cc2 = 0;
                let mut cc3 = 0;
                for i in 0..Self::N {
                    let (d1, ee1) = $crate::utils64::subborrow_u64(a.0[i] ^ sf, sf, cc1);
                    cc1 = ee1;
                    let (d2, ee2) = $crate::utils64::subborrow_u64(b.0[i] ^ sg, sg, cc2);
                    cc2 = ee2;
                    let (d3, hi3) = $crate::utils64::umull_x2_add(d1, f, d2, g, cc3);
                    self.0[i] = d3;
                    cc3 = hi3;
                }
                let up = cc3
                    .wrapping_sub((cc1 as u64).wrapping_neg() & f)
                    .wrapping_sub((cc2 as u64).wrapping_neg() & g);

                // Right-shift the result by 31 bits.
                for i in 0..(Self::N - 1) {
                    self.0[i] = (self.0[i] >> 31) | (self.0[i + 1] << 33);
                }
                self.0[Self::N - 1] = (self.0[Self::N - 1] >> 31) | (up << 33);

                // Negate the result if (a*f + b*g) was negative.
                let w = $crate::utils64::sgnw(up);
                let mut cc = 0;
                for i in 0..Self::N {
                    let (d, ee) = $crate::utils64::subborrow_u64(self.0[i] ^ w, w, cc);
                    self.0[i] = d;
                    cc = ee;
                }

                w
            }

            #[inline(always)]
            fn lindiv31abs(a: &Self, b: &Self, f: u64, g: u64) -> (Self, u64) {
                let mut r = Self::ZERO;
                let ng = r.set_lindiv31abs(a, b, f, g);
                (r, ng)
            }

            /// Divide this value by `y`. If `y` is zero, then this sets this value
            /// to zero.
            fn set_div(&mut self, y: &Self) {
                // a <- y
                // b <- p (modulus)
                // u <- x (self)
                // v <- 0
                //
                // Invariants:
                //    a*x = y*u mod p
                //    b*x = y*v mod p
                //    b is always odd
                //
                // At each step:
                //    if a is even, then:
                //        a <- a/2, u <- u/2 mod p
                //    else:
                //        if a < b:
                //            (a, u, b, v) <- (b, v, a, u)
                //        a <- (a - b)/2
                //        u <- (u - v)/2 mod p
                //
                // We optimize this algorithm following:
                //    https://eprint.iacr.org/2020/972

                let mut a = *y;
                let mut b = Self(Self::MODULUS);
                let mut u = *self;
                let mut v = Self::ZERO;

                // Generic loop; each iteration reduces the sum of the sizes
                // of a and b by at least 31, and that sum starts at 2*BITLEN
                // (at most). We need to run it until the sum of the two lengths
                // is at most 64.
                for _ in 0..Self::NUM1 {
                    // Get approximations of a and b over 64 bits:
                    //  - If len(a) <= 64 and len(b) <= 64, then we just
                    //    use their values (low limbs).
                    //  - Otherwise, with n = max(len(a), len(b)), we use:
                    //       (a mod 2^31) + 2^31*floor(a / 2^(n - 33))
                    //       (b mod 2^31) + 2^31*floor(b / 2^(n - 33))
                    let mut c_hi = 0xFFFFFFFFFFFFFFFFu64;
                    let mut c_lo = 0xFFFFFFFFFFFFFFFFu64;
                    let mut a_hi = 0u64;
                    let mut a_lo = 0u64;
                    let mut b_hi = 0u64;
                    let mut b_lo = 0u64;
                    for j in (0..Self::N).rev() {
                        let aw = a.0[j];
                        let bw = b.0[j];
                        a_hi ^= (a_hi ^ aw) & c_hi;
                        a_lo ^= (a_lo ^ aw) & c_lo;
                        b_hi ^= (b_hi ^ bw) & c_hi;
                        b_lo ^= (b_lo ^ bw) & c_lo;
                        c_lo = c_hi;
                        let mw = aw | bw;
                        c_hi &= ((mw | mw.wrapping_neg()) >> 63).wrapping_sub(1);
                    }

                    // If c_lo = 0, then we grabbed two words for a and b.
                    // If c_lo != 0, then c_hi = 0 (they cannot be both non-zero
                    // since that would mean that a = b = 0, but b is odd). In that
                    // case, we grabbed one word (in a_hi and b_hi) and both values
                    // fit in 64 bits.
                    let s = $crate::utils64::lzcnt(a_hi | b_hi);
                    let mut xa = (a_hi << s) | ((a_lo >> 1) >> (63 - s));
                    let mut xb = (b_hi << s) | ((b_lo >> 1) >> (63 - s));
                    xa = (xa & 0xFFFFFFFF80000000) | (a.0[0] & 0x000000007FFFFFFF);
                    xb = (xb & 0xFFFFFFFF80000000) | (b.0[0] & 0x000000007FFFFFFF);

                    // If c_lo != 0, then we should ignore the computed xa and xb,
                    // and instead use the low limbs directly.
                    xa ^= c_lo & (xa ^ a.0[0]);
                    xb ^= c_lo & (xb ^ b.0[0]);

                    // Compute the 31 inner iterations.
                    let mut fg0 = 1u64;
                    let mut fg1 = 1u64 << 32;
                    for _ in 0..31 {
                        let a_odd = (xa & 1).wrapping_neg();
                        let (_, cc) = $crate::utils64::subborrow_u64(xa, xb, 0);
                        let swap = a_odd & (cc as u64).wrapping_neg();
                        let t1 = swap & (xa ^ xb);
                        xa ^= t1;
                        xb ^= t1;
                        let t2 = swap & (fg0 ^ fg1);
                        fg0 ^= t2;
                        fg1 ^= t2;
                        xa = xa.wrapping_sub(a_odd & xb);
                        fg0 = fg0.wrapping_sub(a_odd & fg1);
                        xa >>= 1;
                        fg1 <<= 1;
                    }
                    fg0 = fg0.wrapping_add(0x7FFFFFFF7FFFFFFF);
                    fg1 = fg1.wrapping_add(0x7FFFFFFF7FFFFFFF);
                    let f0 = (fg0 & 0xFFFFFFFF).wrapping_sub(0x7FFFFFFF);
                    let g0 = (fg0 >> 32).wrapping_sub(0x7FFFFFFF);
                    let f1 = (fg1 & 0xFFFFFFFF).wrapping_sub(0x7FFFFFFF);
                    let g1 = (fg1 >> 32).wrapping_sub(0x7FFFFFFF);

                    // Propagate updates to a, b, u and v.
                    let (na, nega) = Self::lindiv31abs(&a, &b, f0, g0);
                    let (nb, negb) = Self::lindiv31abs(&a, &b, f1, g1);
                    let f0 = (f0 ^ nega).wrapping_sub(nega);
                    let g0 = (g0 ^ nega).wrapping_sub(nega);
                    let f1 = (f1 ^ negb).wrapping_sub(negb);
                    let g1 = (g1 ^ negb).wrapping_sub(negb);
                    let nu = Self::montylin(&u, &v, f0, g0);
                    let nv = Self::montylin(&u, &v, f1, g1);
                    a = na;
                    b = nb;
                    u = nu;
                    v = nv;
                }

                // If y is non-zero, then the final GCD is 1, and
                // len(a) + len(b) <= NUM2 + 2 at this point (initially,
                // len(a) + len(b) <= 2*BITLEN, and each outer iteration reduces
                // the total by at least 31). Thus, the two values fit in one word
                // and we can finish the computation that way. We only need NUM2
                // iterations to reach the point where b = 1.
                let mut xa = a.0[0];
                let mut xb = b.0[0];
                let mut f0 = 1u64;
                let mut g0 = 0u64;
                let mut f1 = 0u64;
                let mut g1 = 1u64;
                for _ in 0..Self::NUM2 {
                    let a_odd = (xa & 1).wrapping_neg();
                    let (_, cc) = $crate::utils64::subborrow_u64(xa, xb, 0);
                    let swap = a_odd & (cc as u64).wrapping_neg();
                    let t1 = swap & (xa ^ xb);
                    xa ^= t1;
                    xb ^= t1;
                    let t2 = swap & (f0 ^ f1);
                    f0 ^= t2;
                    f1 ^= t2;
                    let t3 = swap & (g0 ^ g1);
                    g0 ^= t3;
                    g1 ^= t3;
                    xa = xa.wrapping_sub(a_odd & xb);
                    f0 = f0.wrapping_sub(a_odd & f1);
                    g0 = g0.wrapping_sub(a_odd & g1);
                    xa >>= 1;
                    f1 <<= 1;
                    g1 <<= 1;
                }

                self.set_montylin(&u, &v, f1, g1);

                // If y != 0 then b = 1 at this point. If y == 0, then we
                // force the result to zero.
                let w = !y.is_zero();
                let w = ((w as u64) << 32) | (w as u64);
                for i in 0..Self::N {
                    self.0[i] &= w;
                }

                // At this point, each outer iteration injected 31 extra doublings,
                // plus NUM2 for the last loop, for a total of NUM1*31 + NUM2.
                // Each montylin() call divided by 2^64, so in total we really
                // divided the value by 2^(64*(NUM1+1) - 31*NUM1 - NUM2).
                //
                // Moreover, both divisor and dividend were in Montgomery
                // representation, so the result is not in Montgomery representation
                // (the two R factors canceled each other). We want the result
                // in Montgomery representation, i.e. multiplied by 2^(64*N).
                // Therefore, we must multiply by 2^(33*NUM1 + 64 - NUM2 + 64*N),
                // which we need in
                self.set_mul(&Self::TFIXDIV);
            }

            pub fn set_invert(&mut self) {
                let r = *self;
                *self = Self::ONE;
                self.set_div(&r);
            }

            pub fn invert(self) -> Self {
                let mut r = Self::ONE;
                r.set_div(&self);
                r
            }

            /// Set this value to its square root. Returned value is 0xFFFFFFFF if
            /// the operation succeeded (value was indeed a quadratic residue), or
            /// 0x00000000 otherwise. On success, the chosen root is the one whose
            /// least significant bit (as an integer in [0..p-1]) is zero. On
            /// failure, this value is set to 0.
            pub fn set_sqrt(&mut self) -> u32 {
                // Compute x^((p+1)/4)
                let x = *self;
                self.set_pow_pubexp(&Self::SQRT_EXP);

                // Check whether the square of the result equals the input and zeroize
                // on failure
                let r = self.square().equals(&x);
                let rw = (r as u64) | ((r as u64) << 32);
                for i in 0..Self::N {
                    self.0[i] &= rw;
                }

                // Normalise the output so that the LSB is zero
                let ctl = ((self.encode()[0] as u32) & 1).wrapping_neg();
                self.set_cond_neg(ctl);

                r
            }

            /// Compute the square root of this value. If this value is indeed a
            /// quadratic residue, then this returns (x, 0xFFFFFFFF), with x being
            /// the (unique) square root of this value whose least significant bit
            /// is zero (when normalized to an integer in [0..p-1]). If this value
            /// is not a quadratic residue, then this returns (zero, 0x00000000).
            pub fn sqrt(self) -> (Self, u32) {
                let mut x = self;
                let r = x.set_sqrt();
                (x, r)
            }

            /// Set this value to its square root. Returned value is 0xFFFFFFFF if
            /// the operation succeeded (value was indeed a quadratic residue), or
            /// 0x00000000 otherwise. On success, the chosen root is the one whose
            /// least significant bit (as an integer in [0..p-1]) is zero. On
            /// failure, this value is set to 0.
            ///
            /// When p = 7 mod 8 we can compute x^((p+1)/8), but for all other cases
            /// we fall back to the 2x slower method of computing x^((p+1)/4) twice.
            pub fn set_fourth_root(&mut self) -> u32 {
                let x = *self;

                if Self::MODULUS[0] & 7 == 7 {
                    // Compute x^((p+1)/8)
                    self.set_pow_pubexp(&Self::FOURTH_ROOT_EXP);
                } else {
                    // Fall back to the much slower, general case of two sqrt.
                    self.set_pow_pubexp(&Self::SQRT_EXP);
                    self.set_pow_pubexp(&Self::SQRT_EXP);
                }

                // Check whether the square of the result equals the input and zeroize
                // on failure
                let r = self.n_square(2).equals(&x);
                let rw = (r as u64) | ((r as u64) << 32);
                for i in 0..Self::N {
                    self.0[i] &= rw;
                }

                // Normalise the output so that the LSB is zero
                let ctl = ((self.encode()[0] as u32) & 1).wrapping_neg();
                self.set_cond_neg(ctl);

                r
            }

            /// Compute the fourth root of this value. If this value is indeed some
            /// element to the power of four, then this returns (x, 0xFFFFFFFF), with x being
            /// the (unique) fourth root of this value whose least significant bit
            /// is zero (when normalized to an integer in [0..p-1]). If this value
            /// is not some element to the power of four, then this returns (zero, 0x00000000).
            pub fn fourth_root(self) -> (Self, u32) {
                let mut x = self;
                let r = x.set_fourth_root();
                (x, r)
            }

            /// Raise this value to the provided exponent. The exponent is non-zero
            /// and is public. The exponent is encoded over N 64-bit limbs.
            pub fn set_pow_pubexp(&mut self, e: &[u64; Self::N]) {
                // Make a 4-bit window; win[i] contains x^(i+1)
                let mut win = [Self::ZERO; 15];
                win[0] = *self;
                for i in 1..8 {
                    let j = i * 2;
                    win[j - 1] = win[i - 1].square();
                    win[j] = win[j - 1] * win[0];
                }

                // Explore 4-bit chunks of the exponent, high to low. Skip leading
                // chunks of value 0.
                let mut z = false;
                for i in (0..Self::N).rev() {
                    let ew = e[i];
                    for j in (0..16).rev() {
                        if z {
                            self.set_n_square(4);
                        }
                        let c = ((ew >> (j << 2)) & 0x0F) as usize;
                        if c != 0 {
                            if z {
                                self.set_mul(&win[c - 1]);
                            } else {
                                z = true;
                                *self = win[c - 1];
                            }
                        }
                    }
                }
                if !z {
                    *self = Self::ONE;
                }
            }

            /// Return this value to the provided exponent. The exponent is non-zero
            /// and is public. The exponent is encoded over N 64-bit limbs.
            pub fn pow_pubexp(self, e: &[u64; Self::N]) -> Self {
                let mut r = self;
                r.set_pow_pubexp(e);
                r
            }

            /// Raise this value to the power e. Exponent e is encoded in
            /// unsigned little-endian convention over exactly ebitlen bits.
            pub fn set_pow(&mut self, e: &[u8], ebitlen: usize) {
                self.set_pow_ext(e, 0, ebitlen);
            }

            /// Raise this value to the power e. Exponent e is encoded in
            /// unsigned little-endian convention, over exactly ebitlen bits,
            /// and starting at the bit offset eoff.
            pub fn set_pow_ext(&mut self, e: &[u8], eoff: usize, ebitlen: usize) {
                // TODO: implement a window optimization to make fewer
                // multiplications.
                let x = *self;
                *self = Self::ONE;
                for i in (eoff..(eoff + ebitlen)).rev() {
                    let y = &*self * &x;
                    let ctl = (((e[i >> 3] >> (i & 7)) as u32) & 1).wrapping_neg();
                    self.set_cond(&y, ctl);
                    if i == eoff {
                        break;
                    }
                    self.set_square();
                }
            }

            /// Return this value to the power e (as a new element). Exponent e
            /// is encoded in unsigned little-endian convention over exactly
            /// ebitlen bits.
            pub fn pow(self, e: &[u8], ebitlen: usize) -> Self {
                let mut x = self;
                x.set_pow(e, ebitlen);
                x
            }

            /// Return this value to the power e (as a new element). Exponent e
            /// is encoded in unsigned little-endian convention over exactly
            /// ebitlen bits, and starting at the bit offset eoff.
            pub fn pow_ext(self, e: &[u8], eoff: usize, ebitlen: usize) -> Self {
                let mut x = self;
                x.set_pow_ext(e, eoff, ebitlen);
                x
            }

            /// Raise this value to the power e. The exponent length (in bits)
            /// MUST be at most ebitlen. This is constant-time for both the
            /// base value (self) and the exponent (e); the exponent maximum
            /// size (ebitlen) is considered non-secret.
            fn set_pow_u64(&mut self, e: u64, ebitlen: usize) {
                match ebitlen {
                    0 => {
                        *self = Self::ONE;
                    }
                    1 => {
                        self.set_cond(&Self::ONE, ((e as u32) & 1).wrapping_sub(1));
                    }
                    _ => {
                        let x = *self;
                        self.set_cond(
                            &Self::ONE,
                            (((e >> (ebitlen - 1)) as u32) & 1).wrapping_sub(1),
                        );
                        for i in (0..(ebitlen - 1)).rev() {
                            self.set_square();
                            let y = &*self * &x;
                            self.set_cond(&y, (((e >> i) as u32) & 1).wrapping_neg());
                        }
                    }
                }
            }

            /// Return this value to the power e. The exponent length (in bits)
            /// MUST be at most ebitlen. This is constant-time for both the
            /// base value (self) and the exponent (e); the exponent maximum
            /// size (ebitlen) is considered non-secret.
            fn pow_u64(self, e: u64, ebitlen: usize) -> Self {
                let mut x = self;
                x.set_pow_u64(e, ebitlen);
                x
            }

            /// Raise this value to the power e. The exponent is considered
            /// non-secret.
            fn set_pow_u64_vartime(&mut self, e: u64) {
                match e {
                    0 => {
                        *self = Self::ONE;
                    }
                    1 => {
                        return;
                    }
                    2 => {
                        self.set_square();
                    }
                    3 => {
                        *self *= self.square();
                    }
                    4 => {
                        self.set_square();
                        self.set_square();
                    }
                    _ => {
                        let xx = self.square();
                        let xw = [*self, xx, xx * &*self];
                        let mut j = 63 - e.leading_zeros();
                        j &= !1u32;
                        *self = xw[((e >> j) as usize) - 1];
                        while j > 0 {
                            j -= 2;
                            self.set_square();
                            self.set_square();
                            let k = ((e >> j) as usize) & 3;
                            if k > 0 {
                                self.set_mul(&xw[k - 1]);
                            }
                        }
                    }
                }
            }

            /// Return this value to the power e. The exponent is considered
            /// non-secret.
            fn pow_u64_vartime(self, e: u64) -> Self {
                let mut x = self;
                x.set_pow_u64_vartime(e);
                x
            }

            /// Legendre symbol on this value. Return value is:
            ///   0   if this value is zero
            ///  +1   if this value is a non-zero quadratic residue
            ///  -1   if this value is not a quadratic residue
            pub fn legendre(self) -> i32 {
                // This is the same optimized binary GCD as in division, except
                // that we do not need to keep track of u and v. We can also
                // work directly on the Montgomery representation because R = 2^1184
                // is a square.
                let mut a = self;
                let mut b = Self(Self::MODULUS);
                let mut ls = 0u64;

                // Generic loop; each iteration reduces the sum of the sizes
                // of a and b by at least 31, and that sum starts at 2*BITLEN
                // (at most). We need to run it until the sum of the two lengths
                // is at most 64.
                for _ in 0..Self::NUM1 {
                    // Get approximations of a and b over 64 bits:
                    //  - If len(a) <= 64 and len(b) <= 64, then we just
                    //    use their values (low limbs).
                    //  - Otherwise, with n = max(len(a), len(b)), we use:
                    //       (a mod 2^31) + 2^31*floor(a / 2^(n - 33))
                    //       (b mod 2^31) + 2^31*floor(b / 2^(n - 33))
                    let mut c_hi = 0xFFFFFFFFFFFFFFFFu64;
                    let mut c_lo = 0xFFFFFFFFFFFFFFFFu64;
                    let mut a_hi = 0u64;
                    let mut a_lo = 0u64;
                    let mut b_hi = 0u64;
                    let mut b_lo = 0u64;
                    for j in (0..Self::N).rev() {
                        let aw = a.0[j];
                        let bw = b.0[j];
                        a_hi ^= (a_hi ^ aw) & c_hi;
                        a_lo ^= (a_lo ^ aw) & c_lo;
                        b_hi ^= (b_hi ^ bw) & c_hi;
                        b_lo ^= (b_lo ^ bw) & c_lo;
                        c_lo = c_hi;
                        let mw = aw | bw;
                        c_hi &= ((mw | mw.wrapping_neg()) >> 63).wrapping_sub(1);
                    }

                    // If c_lo = 0, then we grabbed two words for a and b.
                    // If c_lo != 0, then c_hi = 0 (they cannot be both non-zero
                    // since that would mean that a = b = 0, but b is odd). In that
                    // case, we grabbed one word (in a_hi and b_hi) and both values
                    // fit in 64 bits.
                    let s = $crate::utils64::lzcnt(a_hi | b_hi);
                    let mut xa = (a_hi << s) | ((a_lo >> 1) >> (63 - s));
                    let mut xb = (b_hi << s) | ((b_lo >> 1) >> (63 - s));
                    xa = (xa & 0xFFFFFFFF80000000) | (a.0[0] & 0x000000007FFFFFFF);
                    xb = (xb & 0xFFFFFFFF80000000) | (b.0[0] & 0x000000007FFFFFFF);

                    // If c_lo != 0, then we should ignore the computed xa and xb,
                    // and instead use the low limbs directly.
                    xa ^= c_lo & (xa ^ a.0[0]);
                    xb ^= c_lo & (xb ^ b.0[0]);

                    // First 29 inner iterations.
                    let mut fg0 = 1u64;
                    let mut fg1 = 1u64 << 32;
                    for _ in 0..29 {
                        let a_odd = (xa & 1).wrapping_neg();
                        let (_, cc) = $crate::utils64::subborrow_u64(xa, xb, 0);
                        let swap = a_odd & (cc as u64).wrapping_neg();
                        ls ^= swap & ((xa & xb) >> 1);
                        let t1 = swap & (xa ^ xb);
                        xa ^= t1;
                        xb ^= t1;
                        let t2 = swap & (fg0 ^ fg1);
                        fg0 ^= t2;
                        fg1 ^= t2;
                        xa = xa.wrapping_sub(a_odd & xb);
                        fg0 = fg0.wrapping_sub(a_odd & fg1);
                        xa >>= 1;
                        fg1 <<= 1;
                        ls ^= xb.wrapping_add(2) >> 2;
                    }

                    // Compute the updated a and b (low words only) to get enough
                    // bits for the next two iterations.
                    let fg0z = fg0.wrapping_add(0x7FFFFFFF7FFFFFFF);
                    let fg1z = fg1.wrapping_add(0x7FFFFFFF7FFFFFFF);
                    let f0 = (fg0z & 0xFFFFFFFF).wrapping_sub(0x7FFFFFFF);
                    let g0 = (fg0z >> 32).wrapping_sub(0x7FFFFFFF);
                    let f1 = (fg1z & 0xFFFFFFFF).wrapping_sub(0x7FFFFFFF);
                    let g1 = (fg1z >> 32).wrapping_sub(0x7FFFFFFF);
                    let mut a0 = a.0[0]
                        .wrapping_mul(f0)
                        .wrapping_add(b.0[0].wrapping_mul(g0))
                        >> 29;
                    let mut b0 = a.0[0]
                        .wrapping_mul(f1)
                        .wrapping_add(b.0[0].wrapping_mul(g1))
                        >> 29;
                    for _ in 0..2 {
                        let a_odd = (xa & 1).wrapping_neg();
                        let (_, cc) = $crate::utils64::subborrow_u64(xa, xb, 0);
                        let swap = a_odd & (cc as u64).wrapping_neg();
                        ls ^= swap & ((a0 & b0) >> 1);
                        let t1 = swap & (xa ^ xb);
                        xa ^= t1;
                        xb ^= t1;
                        let t2 = swap & (fg0 ^ fg1);
                        fg0 ^= t2;
                        fg1 ^= t2;
                        let t3 = swap & (a0 ^ b0);
                        a0 ^= t3;
                        b0 ^= t3;
                        xa = xa.wrapping_sub(a_odd & xb);
                        fg0 = fg0.wrapping_sub(a_odd & fg1);
                        a0 = a0.wrapping_sub(a_odd & b0);
                        xa >>= 1;
                        fg1 <<= 1;
                        a0 >>= 1;
                        ls ^= b0.wrapping_add(2) >> 2;
                    }

                    // Propagate updates to a and b.
                    fg0 = fg0.wrapping_add(0x7FFFFFFF7FFFFFFF);
                    fg1 = fg1.wrapping_add(0x7FFFFFFF7FFFFFFF);
                    let f0 = (fg0 & 0xFFFFFFFF).wrapping_sub(0x7FFFFFFF);
                    let g0 = (fg0 >> 32).wrapping_sub(0x7FFFFFFF);
                    let f1 = (fg1 & 0xFFFFFFFF).wrapping_sub(0x7FFFFFFF);
                    let g1 = (fg1 >> 32).wrapping_sub(0x7FFFFFFF);

                    // Propagate updates to a, b, u and v.
                    let (na, nega) = Self::lindiv31abs(&a, &b, f0, g0);
                    let (nb, _) = Self::lindiv31abs(&a, &b, f1, g1);
                    ls ^= nega & (nb.0[0] >> 1);
                    a = na;
                    b = nb;
                }

                // If y is non-zero, then the final GCD is 1, and
                // len(a) + len(b) <= NUM2 + 2 at this point (initially,
                // len(a) + len(b) <= 2*BITLEN, and each outer iteration reduces
                // the total by at least 31). Thus, the two values fit in one word
                // and we can finish the computation that way. We only need NUM2
                // iterations to reach the point where b = 1.
                let mut xa = a.0[0];
                let mut xb = b.0[0];
                for _ in 0..Self::NUM2 {
                    let a_odd = (xa & 1).wrapping_neg();
                    let (_, cc) = $crate::utils64::subborrow_u64(xa, xb, 0);
                    let swap = a_odd & (cc as u64).wrapping_neg();
                    ls ^= swap & ((xa & xb) >> 1);
                    let t1 = swap & (xa ^ xb);
                    xa ^= t1;
                    xb ^= t1;
                    xa = xa.wrapping_sub(a_odd & xb);
                    xa >>= 1;
                    ls ^= xb.wrapping_add(2) >> 2;
                }

                // At this point, if the source value was not zero, then the low
                // bit of ls contains the QR status (0 = square, 1 = non-square),
                // which we need to convert to the expected value (+1 or -1).
                // If y == 0, then we return 0, per the API.
                let r = 1u32.wrapping_sub(((ls as u32) & 1) << 1);
                (r & !(self.is_zero() as u32)) as i32
            }

            /// Return `0xFFFFFFFF` when this value is a square in GF(p^2) and
            /// `0x00000000` otherwise.
            #[inline]
            fn is_square(self) -> u32 {
                !((self.legendre() >> 1) as u32)
            }

            /// Encode this value into bytes. Encoding uses little-endian, has
            /// a fixed size (for a given field), and is canonical.
            #[inline(always)]
            pub fn encode(self) -> [u8; Self::ENCODED_LENGTH] {
                let mut r = self;
                r.set_montyred();
                let mut d = [0u8; Self::ENCODED_LENGTH];
                for i in 0..(Self::N - 1) {
                    d[(i * 8)..(i * 8 + 8)].copy_from_slice(&r.0[i].to_le_bytes());
                }
                d[((Self::N - 1) * 8)..].copy_from_slice(
                    &(r.0[Self::N - 1].to_le_bytes()[..Self::ENCODED_LENGTH - (Self::N - 1) * 8]),
                );
                d
            }

            #[inline]
            fn set_decode_nocheck(&mut self, buf: &[u8]) {
                for i in 0..(Self::N - 1) {
                    self.0[i] = u64::from_le_bytes(
                        *<&[u8; 8]>::try_from(&buf[(8 * i)..(8 * i + 8)]).unwrap(),
                    );
                }
                let mut w = 0u64;
                for j in 0..(Self::ENCODED_LENGTH - (Self::N - 1) * 8) {
                    w |= (buf[(Self::N - 1) * 8 + j] as u64) << (8 * j);
                }
                self.0[Self::N - 1] = w;
            }

            /// Decode the provided bytes into a field element. Returned values
            /// are the element and 0xFFFFFFFF on success, or the zero element and
            /// 0x00000000 on failure. A failure is reported if the source slice
            /// does not have exactly the canonical encoding length of a field
            /// element (Self::ENCODED_LENGTH), or if the source encodes
            /// an integer which is not in the [0..(p-1)] range.
            #[inline(always)]
            pub fn decode(buf: &[u8]) -> (Self, u32) {
                if buf.len() != Self::ENCODED_LENGTH {
                    return (Self::ZERO, 0);
                }

                // decode raw value
                let mut r = Self::ZERO;
                r.set_decode_nocheck(buf);

                // check that the source is canonical; clear if invalid
                let (_, mut cc) = $crate::utils64::subborrow_u64(r.0[0], Self::MODULUS[0], 0);
                for i in 1..Self::N {
                    let (_, ee) = $crate::utils64::subborrow_u64(r.0[i], Self::MODULUS[i], cc);
                    cc = ee;
                }
                let m = (cc as u64).wrapping_neg();
                for i in 0..Self::N {
                    r.0[i] &= m;
                }

                // convert to Montgomery representation
                r.set_mul(&Self::R2);
                (r, m as u32)
            }

            /// Set this element by decoding the provided bytes. The source slice
            /// can have arbitrary length; the bytes are interpreted with the
            /// unsigned little-endian convention (no sign bit), and the resulting
            /// integer is reduced modulo the field modulus p. By definition, this
            /// function does not enforce canonicality of the source value.
            #[inline]
            pub fn set_decode_reduce(&mut self, buf: &[u8]) {
                let mut n = buf.len();
                if n == 0 {
                    *self = Self::ZERO;
                    return;
                }

                let mut tmp = [0u8; Self::ENCODED_LENGTH];
                let mut nn = n % Self::CLEN;
                if nn == 0 {
                    nn = Self::CLEN;
                }
                n -= nn;
                tmp[..nn].copy_from_slice(&buf[n..]);
                self.set_decode_nocheck(&tmp);

                while n > 0 {
                    n -= Self::CLEN;
                    tmp[..Self::CLEN].copy_from_slice(&buf[n..(n + Self::CLEN)]);
                    let mut d = Self::ZERO;
                    d.set_decode_nocheck(&tmp);
                    self.set_mul(&Self::TDEC);
                    self.set_add(&d);
                }

                self.set_mul(&Self::R2);
            }

            /// Decode the provided bytes into a field element. The source slice
            /// can have arbitrary length; the bytes are interpreted with the
            /// unsigned little-endian convention (no sign bit), and the resulting
            /// integer is reduced modulo the field modulus p. By definition, this
            /// function does not enforce canonicality of the source value.
            #[inline(always)]
            pub fn decode_reduce(buf: &[u8]) -> Self {
                let mut x = Self::ZERO;
                x.set_decode_reduce(buf);
                x
            }

            /// Set this structure to a random field element (indistinguishable
            /// from uniform generation).
            pub fn set_rand<T: ::rand_core::CryptoRng + ::rand_core::RngCore>(
                &mut self,
                rng: &mut T,
            ) {
                let mut tmp = [0u8; Self::ENCODED_LENGTH + 16];
                rng.fill_bytes(&mut tmp);
                self.set_decode_reduce(&tmp);
            }

            /// Return a new random field element (indistinguishable from
            /// uniform generation).
            pub fn rand<T: ::rand_core::CryptoRng + ::rand_core::RngCore>(rng: &mut T) -> Self {
                let mut x = Self::ZERO;
                x.set_rand(rng);
                x
            }

            /// Get the "hash" of the value (low 64 bits of the Montgomery
            /// representation).
            pub fn hashcode(self) -> u64 {
                self.0[0]
            }

            /// Implements Algorithm 2 from Patrick Longa's
            /// [ePrint 2022-367](https://eprint.iacr.org/2022/367) §3.
            /// Computes a1 * b1 + a2 * b2 using an optimised method intended
            /// for use in Fp2 multiplications
            #[inline(always)]
            pub fn sum_of_products(a1: &Self, b1: &Self, a2: &Self, b2: &Self) -> Self {
                // Line 1: u <- 0
                let mut u = Self::ZERO;

                let mut cc1: u64;
                let mut cc2: u64;
                let mut cc3: u64;

                let mut cch: u8 = 0;
                let mut cch1: u8 = 0;
                let mut cch2: u8 = 0;

                // Line 2: for j = 0 to N - 1
                for j in 0..Self::N {
                    // Line 3: u <- u + a0,j * b0 + a1,j * b1

                    // We can always compute a * b + c + d = lo + 2^64 * hi
                    // This is what we do in the inner loop to compute
                    //     u <- u + a0,j * b0
                    (u.0[0], cc1) = $crate::utils64::umull_add(a1.0[j], b1.0[0], u.0[0]);
                    for k in 1..Self::N {
                        (u.0[k], cc1) = $crate::utils64::umull_add2(a1.0[j], b1.0[k], u.0[k], cc1);
                    }
                    //     u <- u + a1,j * b1
                    (u.0[0], cc2) = $crate::utils64::umull_add(a2.0[j], b2.0[0], u.0[0]);
                    for k in 1..Self::N {
                        (u.0[k], cc2) = $crate::utils64::umull_add2(a2.0[j], b2.0[k], u.0[k], cc2);
                    }

                    // Line 4: q <- u * p' mod 2^64
                    let q = u.0[0].wrapping_mul(Self::P0I);

                    // Line 5: u <- (u + q * p') / 2^64
                    (_, cc3) = $crate::utils64::umull_add(q, Self::MODULUS[0], u.0[0]);
                    for k in 1..Self::N {
                        (u.0[k - 1], cc3) =
                            $crate::utils64::umull_add2(q, Self::MODULUS[k], u.0[k], cc3);
                    }

                    // We now have to handle all the carries, which means adding cc1, cc2 and cc3 as well
                    // as the carry cch from the last iteration
                    (u.0[Self::N - 1], cch1) = $crate::utils64::addcarry_u64(cc1, cc2, cch);
                    (u.0[Self::N - 1], cch2) =
                        $crate::utils64::addcarry_u64(u.0[Self::N - 1], cc3, 0);
                    cch = cch1 + cch2;
                }

                // From the paper we have something in the range [0, 2p) if 2*(p - 1)^2 < p*R
                // Which means we need to subtract p and then conditionally add p back
                let mut borrow: u8 = 0;
                for i in 0..Self::N {
                    (u.0[i], borrow) =
                        $crate::utils64::subborrow_u64(u.0[i], Self::MODULUS[i], borrow);
                }
                let mask = (cch as u64).wrapping_sub(borrow as u64);
                let mut cc = 0;
                for i in 0..Self::N {
                    (u.0[i], cc) =
                        $crate::utils64::addcarry_u64(u.0[i], mask & Self::MODULUS[i], cc);
                }

                // We are only guarenteed that the above result is in [0, 2p) when we have the
                // condition: 2*(p - 1)^2 < p*R which is true most of the time, but not when
                // p is close to R. When the user's modulus is too close to 2^(64 * N) then we
                // may need an additional conditional subtraction to get a result within the canonical
                // range.
                if Self::SUM_OF_PRODUCTS_ADDITIONAL_SUB {
                    borrow = 0;
                    for i in 0..Self::N {
                        (u.0[i], borrow) =
                            $crate::utils64::subborrow_u64(u.0[i], Self::MODULUS[i], borrow);
                    }
                    let mask = (borrow as u64).wrapping_neg();
                    cc = 0;
                    for i in 0..Self::N {
                        (u.0[i], cc) =
                            $crate::utils64::addcarry_u64(u.0[i], mask & Self::MODULUS[i], cc);
                    }
                }

                u
            }

            /// Implements Algorithm 2 from Patrick Longa's
            /// [ePrint 2022-367](https://eprint.iacr.org/2022/367) §3.
            /// Computes a1 * b1 - a2 * b2 using an optimised method intended
            /// for use in Fp2 multiplications
            #[inline(always)]
            pub fn difference_of_products(a1: &Self, b1: &Self, a2: &Self, b2: &Self) -> Self {
                // Compute p - b2
                let mut borrow = 0;
                let mut b2_minus = *b2;
                for i in 0..Self::N {
                    (b2_minus.0[i], borrow) =
                        $crate::utils64::subborrow_u64(Self::MODULUS[i], b2_minus.0[i], borrow);
                }
                // Regular sum of products (could the above be optimised into the main loop?)
                Self::sum_of_products(a1, b1, a2, &b2_minus)
            }

            pub fn batch_invert(xx: &mut [Self]) {
                // We use Montgomery's trick:
                //   1/u = v*(1/(u*v))
                //   1/v = u*(1/(u*v))
                // Applied recursively on n elements, this computes an inversion
                // with a single inversion in the field, and 3*(n-1) multiplications.
                // We use batches of 200 elements; larger batches only yield
                // moderate improvements, while sticking to a fixed moderate batch
                // size allows stack-based allocation.
                let n = xx.len();
                let mut i = 0;
                while i < n {
                    let blen = if (n - i) > 200 { 200 } else { n - i };
                    let mut tt = [Self::ZERO; 200];
                    tt[0] = xx[i];
                    let zz0 = tt[0].is_zero();
                    tt[0].set_cond(&Self::ONE, zz0);
                    for j in 1..blen {
                        tt[j] = xx[i + j];
                        tt[j].set_cond(&Self::ONE, tt[j].is_zero());
                        tt[j] *= tt[j - 1];
                    }
                    let mut k = Self::ONE / tt[blen - 1];
                    for j in (1..blen).rev() {
                        let mut x = xx[i + j];
                        let zz = x.is_zero();
                        x.set_cond(&Self::ONE, zz);
                        xx[i + j].set_cond(&(k * tt[j - 1]), !zz);
                        k *= x;
                    }
                    xx[i].set_cond(&k, !zz0);
                    i += blen;
                }
            }

            /*
             * Support functions which compute constants at compile time to
             * generate macro constants, this allows a user to create the field
             * with only the modulus as input.
             *
             * These have been adapted from https://github.com/pornin/crrl/src/backend/w64/gfgen.rs
             */
            // Return -1/x mod 2^64. It is assumed that x is odd.
            const fn ninv64(x: u64) -> u64 {
                let y = 2u64.wrapping_sub(x);
                let y = y.wrapping_mul(2u64.wrapping_sub(y.wrapping_mul(x)));
                let y = y.wrapping_mul(2u64.wrapping_sub(y.wrapping_mul(x)));
                let y = y.wrapping_mul(2u64.wrapping_sub(y.wrapping_mul(x)));
                let y = y.wrapping_mul(2u64.wrapping_sub(y.wrapping_mul(x)));
                let y = y.wrapping_mul(2u64.wrapping_sub(y.wrapping_mul(x)));
                y.wrapping_neg()
            }

            // Custom add-with-carry, for use in const (compile-time) contexts.
            const fn adc(x: u64, y: u64, cc: u64) -> (u64, u64) {
                let z = (x as u128) + (y as u128) + (cc as u128);
                (z as u64, (z >> 64) as u64)
            }

            // Custom sub-with-borrow, for use in const (compile-time) contexts.
            const fn sbb(x: u64, y: u64, cc: u64) -> (u64, u64) {
                let z = (x as u128).wrapping_sub(y as u128).wrapping_sub(cc as u128);
                (z as u64, ((z >> 64) as u64).wrapping_neg())
            }

            // Add the modulus, return borrow (compile-time).
            const fn addm(a: Self) -> (Self, u64) {
                let mut r = a;
                let mut cc = 0u64;
                let mut i = 0;
                while i < Self::N {
                    let (d, c) = Self::adc(r.0[i], Self::MODULUS[i], cc);
                    r.0[i] = d;
                    cc = c;
                    i += 1;
                }
                (r, cc)
            }

            // Subtract the modulus, return borrow (compile-time).
            const fn subm(a: Self) -> (Self, u64) {
                let mut r = a;
                let mut cc = 0u64;
                let mut i = 0;
                while i < Self::N {
                    let (d, c) = Self::sbb(r.0[i], Self::MODULUS[i], cc);
                    r.0[i] = d;
                    cc = c;
                    i += 1;
                }
                (r, cc)
            }

            // For the result of sum of products to work, we need 2*p - 4 to fit into
            // N words, otherwise we need to do an additional subtraction. This function
            // computes 2*p - 4 and if an overflow is detected, then a boolean is set
            // at compile time to ensure an addtional conditional subtraction is performed.
            const fn sum_of_products_check() -> bool {
                const fn sub4_inner(a: $typename, cc: u64, i: usize) -> ($typename, u64) {
                    if i == a.0.len() {
                        (a, cc)
                    } else {
                        let n = if i == 0 { 4 } else { 0 };
                        let (d, cc) = $typename::sbb(a.0[i], n, cc);
                        let mut aa = a;
                        aa.0[i] = d;
                        sub4_inner(aa, cc, i + 1)
                    }
                }

                // Subtract 4 from the modulus, we ignore the borrow here
                // as it would only be an issue for p = 2, 3 and we assume
                // p is large.
                let (r1, _) = sub4_inner(Self::new(Self::MODULUS), 0, 0);

                // Add the modulus to compute 2*p - 4
                let (r2, b1) = Self::addm(r1);

                // If the addition of the modulus has a carry, then return true
                b1 == 1
            }

            // Add the modulus if mm == -1; return a unchanged with mm == 0
            // (compile-time).
            const fn addm_cond(a: Self, mm: u64) -> Self {
                let mut r = a;
                let mut cc = 0u64;
                let mut i = 0;
                while i < Self::N {
                    let (d, c) = Self::adc(r.0[i], Self::MODULUS[i] & mm, cc);
                    r.0[i] = d;
                    cc = c;
                    i += 1;
                }
                r
            }

            // Get index of the top non-zero word of the modulus
            // (from the parameters).
            const fn top_word_index() -> usize {
                const fn top_word_index_inner(j: usize) -> usize {
                    if $modulus[j] != 0 {
                        j
                    } else {
                        top_word_index_inner(j - 1)
                    }
                }
                top_word_index_inner($modulus.len() - 1)
            }

            // Compute the modulus exact bit length (compile-time).
            const fn mod_bitlen() -> usize {
                const fn bitlen(x: u64, max: usize) -> usize {
                    if max == 1 {
                        x as usize
                    } else {
                        let hm = max >> 1;
                        let y = x >> hm;
                        if y == 0 {
                            bitlen(x, hm)
                        } else {
                            bitlen(y, max - hm) + hm
                        }
                    }
                }
                (Self::N - 1) * 64 + bitlen($typename::MODULUS[Self::N - 1], 64)
            }

            // Get the top 32 bits of the actual modulus value (if the modulus
            // is less than 32 bits in length, then this returns the modulus).
            const fn top_u32() -> u64 {
                if Self::BIT_LENGTH < 32 {
                    Self::MODULUS[0]
                } else {
                    let hi = Self::MODULUS[Self::N - 1];
                    let bl = Self::BIT_LENGTH & 63;
                    if bl == 0 {
                        hi >> 32
                    } else if bl < 32 {
                        let lo = Self::MODULUS[Self::N - 2];
                        (hi << (32 - bl)) | (lo >> (bl + 32))
                    } else {
                        hi >> (bl - 32)
                    }
                }
            }

            // Compute 2^n mod p as a plain integer (not in Montgomery form).
            const fn pow2mod(n: usize) -> Self {
                // One modular doubling: (2*a) mod p, returned as a plain integer.
                const fn double1(a: $typename) -> $typename {
                    // Left-shift the limb array by one bit.
                    let mut r = a;
                    let mut carry = 0u64;
                    let mut i = 0;
                    while i < $typename::N {
                        let w = r.0[i];
                        r.0[i] = (w << 1) | carry;
                        carry = w >> 63;
                        i += 1;
                    }
                    let dh = carry; // top bit shifted out of the most-significant limb

                    // Subtract the modulus (may underflow).
                    let (r2, cc) = $typename::subm(r);

                    // Add the modulus back if the subtraction underflowed and dh == 0.
                    $typename::addm_cond(r2, (cc & !dh).wrapping_neg())
                }

                let bl = Self::mod_bitlen();

                // Fast path: 2^n < p, so it is already fully reduced.
                if n < bl {
                    let mut d = Self([0u64; Self::N]);
                    d.0[n >> 6] = 1u64 << (n & 63);
                    return d;
                }

                // Start from 2^(bl-1) — fits in one limb, no reduction needed.
                let mut d = Self([0u64; Self::N]);
                d.0[(bl - 1) >> 6] = 1u64 << ((bl - 1) & 63);

                // Apply the remaining doublings one at a time.
                // Maximum iterations: for R2 of a 4096-bit field ≈ 4097.
                let mut remaining = n - (bl - 1);
                while remaining > 0 {
                    d = double1(d);
                    remaining -= 1;
                }

                d
            }

            // Const implementation of modular negation. This MUST NOT be
            // applied on zero.
            const fn const_neg(a: Self) -> Self {
                let mut r = Self([0u64; Self::N]);
                let mut cc = 0u64;
                let mut j = 0;
                while j < Self::N {
                    let (x, c) = Self::sbb(Self::MODULUS[j], a.0[j], cc);
                    r.0[j] = x;
                    cc = c;
                    j += 1;
                }
                r
            }

            // Const implementation of Montgomery multiplication. It uses
            // recursion in order to be compatible with the constraints of
            // const code; at runtime, it would be slower than the normal
            // implementation, but still constant-time (in case it gets
            // mistakenly used).
            const fn const_mmul(a: Self, b: Self) -> Self {
                // 128-bit fused multiply-accumulate: x*y + acc + carry → (lo, hi).
                const fn umaal(x: u64, y: u64, acc: u64, carry: u64) -> (u64, u64) {
                    let z = (x as u128) * (y as u128) + (acc as u128) + (carry as u128);
                    (z as u64, (z >> 64) as u64)
                }

                let mut d = Self([0u64; Self::N]);
                let mut dh: u64 = 0;

                // Outer loop — one Montgomery reduction step per limb of b.
                let mut j = 0;
                while j < Self::N {
                    let bj = b.0[j];
                    // Montgomery factor for this step.
                    let fm = a.0[0]
                        .wrapping_mul(bj)
                        .wrapping_add(d.0[0])
                        .wrapping_mul(Self::P0I);

                    let mut cc1 = 0u64; // carry from a[i]*bj column
                    let mut cc2 = 0u64; // carry from MODULUS[i]*fm column

                    // Inner loop — accumulate a[i]*bj and reduce with fm*MODULUS[i].
                    // The word at i=0 is consumed by the reduction and not stored;
                    // for i>0 each result shifts into d[i-1].
                    let mut i = 0;
                    while i < Self::N {
                        let (z, c1) = umaal(a.0[i], bj, d.0[i], cc1);
                        let (z, c2) = umaal(Self::MODULUS[i], fm, z, cc2);
                        if i > 0 {
                            d.0[i - 1] = z;
                        }
                        cc1 = c1;
                        cc2 = c2;
                        i += 1;
                    }

                    // Fold the two high-word carries into d[N-1] and dh.
                    let (z, zh1) = Self::adc(dh, cc1, 0);
                    let (z, zh2) = Self::adc(z, cc2, 0);
                    d.0[Self::N - 1] = z;
                    dh = zh1 + zh2;

                    j += 1;
                }

                // Final conditional subtraction (identical logic to original).
                let (d2, cc) = Self::subm(d);
                Self::addm_cond(d2, (cc & !dh).wrapping_neg())
            }

            const fn const_small(x: u64) -> Self {
                let mut d = [0u64; Self::N];
                d[0] = x;
                Self::const_mmul(Self(d), Self::R2)
            }

            const fn const_sqrt_exp() -> [u64; Self::N] {
                let mut d = [0u64; Self::N];
                let mut dd = 0u64;
                let mut cc = 1u64; // start with carry=1 to compute (p+1)/4
                let mut i = 0;
                while i < Self::N {
                    let (x, c) = Self::adc(Self::MODULUS[i], 0, cc);
                    cc = c;
                    if i > 0 {
                        d[i - 1] = dd | (x << 62);
                    }
                    dd = x >> 2;
                    i += 1;
                }
                d[Self::N - 1] = dd;
                d
            }

            const fn const_fourth_root_exp() -> [u64; Self::N] {
                let mut d = [0u64; Self::N];
                let mut dd = 0u64;
                let mut cc = 1u64;
                let mut i = 0;
                while i < Self::N {
                    let (x, c) = Self::adc(Self::MODULUS[i], 0, cc);
                    cc = c;
                    if i > 0 {
                        d[i - 1] = dd | (x << 61);
                    }
                    dd = x >> 3;
                    i += 1;
                }
                d[Self::N - 1] = dd;
                d
            }

            /// Decode an element from bytes, no check is made that the input
            /// value is reduced except that the buffer is of the excpected
            /// length of `Self::ENCODED_LENGTH`.
            pub const fn const_decode_no_check(buf: &[u8]) -> Self {
                let mut r = Self::ZERO;
                if buf.len() != Self::ENCODED_LENGTH {
                    return r;
                }

                // Fill the first N-1 elements
                let mut i = 0;
                while i < Self::N - 1 {
                    r.0[i] = u64::from_le_bytes([
                        buf[i * 8],
                        buf[i * 8 + 1],
                        buf[i * 8 + 2],
                        buf[i * 8 + 3],
                        buf[i * 8 + 4],
                        buf[i * 8 + 5],
                        buf[i * 8 + 6],
                        buf[i * 8 + 7],
                    ]);
                    i += 1;
                }

                // Fill the last element
                let mut w = 0u64;
                let mut j = 0;
                while j < Self::ENCODED_LENGTH - (Self::N - 1) * 8 {
                    w |= (buf[(Self::N - 1) * 8 + j] as u64) << (8 * j);
                    j += 1;
                }
                r.0[Self::N - 1] = w;

                Self::const_mmul(r, Self::R2)
            }
        }

        /*
         * Implementations of from methods from simple integer types
         */

        impl From<u64> for $typename {
            fn from(x: u64) -> $typename {
                let mut r = Self::ZERO;
                if Self::N == 1 {
                    r.set_decode_reduce(&x.to_le_bytes());
                } else {
                    r.0[0] = x;
                    r.set_mul(&Self::R2);
                }
                r
            }
        }

        impl From<i64> for $typename {
            fn from(x: i64) -> $typename {
                let sx = (x >> 63) as u64;
                let ax = ((x as u64) ^ sx).wrapping_sub(sx);
                let mut r = Self::from(ax);
                r.set_cond_neg(sx as u32);
                r
            }
        }

        impl From<u32> for $typename {
            fn from(x: u32) -> $typename {
                Self::from(x as u64)
            }
        }

        impl From<i32> for $typename {
            fn from(x: i32) -> $typename {
                let mut r = Self::ONE;
                r.set_mul_small(x);
                r
            }
        }

        /*
         * Implementations of all the traits needed to use the simple operators
         * (+, *, /...) on field element instances, with or without references.
         * as well as a display method.
         */

        impl ::std::fmt::Display for $typename {
            fn fmt(&self, f: &mut ::std::fmt::Formatter) -> ::std::fmt::Result {
                // Encode the value to get a canonical value
                let v_bytes = self.encode();

                // Convert the bytes into a base-16 string
                let mut v_hex_string: String = v_bytes
                    .iter()
                    .rev()
                    .map(|byte| format!("{:02x}", byte)) // Format each byte as a two-digit hex
                    .collect();

                // Remove leading zeros
                v_hex_string = v_hex_string.trim_start_matches("0").to_string();

                // If the value was zero, we need to add a zero back
                if v_hex_string.is_empty() {
                    v_hex_string = String::from("0");
                }

                write!(f, "0x{}", v_hex_string)
            }
        }

        impl ::core::ops::Add<$typename> for $typename {
            type Output = $typename;

            #[inline(always)]
            fn add(self, other: $typename) -> $typename {
                let mut r = self;
                r.set_add(&other);
                r
            }
        }

        impl ::core::ops::Add<&$typename> for $typename {
            type Output = $typename;

            #[inline(always)]
            fn add(self, other: &$typename) -> $typename {
                let mut r = self;
                r.set_add(other);
                r
            }
        }

        impl ::core::ops::Add<$typename> for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn add(self, other: $typename) -> $typename {
                let mut r = *self;
                r.set_add(&other);
                r
            }
        }

        impl ::core::ops::Add<&$typename> for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn add(self, other: &$typename) -> $typename {
                let mut r = *self;
                r.set_add(other);
                r
            }
        }

        impl ::core::ops::AddAssign<$typename> for $typename {
            #[inline(always)]
            fn add_assign(&mut self, other: $typename) {
                self.set_add(&other);
            }
        }

        impl ::core::ops::AddAssign<&$typename> for $typename {
            #[inline(always)]
            fn add_assign(&mut self, other: &$typename) {
                self.set_add(other);
            }
        }

        impl ::core::ops::Div<$typename> for $typename {
            type Output = $typename;

            #[inline(always)]
            fn div(self, other: $typename) -> $typename {
                let mut r = self;
                r.set_div(&other);
                r
            }
        }

        impl ::core::ops::Div<&$typename> for $typename {
            type Output = $typename;

            #[inline(always)]
            fn div(self, other: &$typename) -> $typename {
                let mut r = self;
                r.set_div(other);
                r
            }
        }

        impl ::core::ops::Div<$typename> for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn div(self, other: $typename) -> $typename {
                let mut r = *self;
                r.set_div(&other);
                r
            }
        }

        impl ::core::ops::Div<&$typename> for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn div(self, other: &$typename) -> $typename {
                let mut r = *self;
                r.set_div(other);
                r
            }
        }

        impl ::core::ops::DivAssign<$typename> for $typename {
            #[inline(always)]
            fn div_assign(&mut self, other: $typename) {
                self.set_div(&other);
            }
        }

        impl ::core::ops::DivAssign<&$typename> for $typename {
            #[inline(always)]
            fn div_assign(&mut self, other: &$typename) {
                self.set_div(other);
            }
        }

        impl ::core::ops::Mul<$typename> for $typename {
            type Output = $typename;

            #[inline(always)]
            fn mul(self, other: $typename) -> $typename {
                let mut r = self;
                r.set_mul(&other);
                r
            }
        }

        impl ::core::ops::Mul<&$typename> for $typename {
            type Output = $typename;

            #[inline(always)]
            fn mul(self, other: &$typename) -> $typename {
                let mut r = self;
                r.set_mul(other);
                r
            }
        }

        impl ::core::ops::Mul<$typename> for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn mul(self, other: $typename) -> $typename {
                let mut r = *self;
                r.set_mul(&other);
                r
            }
        }

        impl ::core::ops::Mul<&$typename> for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn mul(self, other: &$typename) -> $typename {
                let mut r = *self;
                r.set_mul(other);
                r
            }
        }

        impl ::core::ops::MulAssign<$typename> for $typename {
            #[inline(always)]
            fn mul_assign(&mut self, other: $typename) {
                self.set_mul(&other);
            }
        }

        impl ::core::ops::MulAssign<&$typename> for $typename {
            #[inline(always)]
            fn mul_assign(&mut self, other: &$typename) {
                self.set_mul(other);
            }
        }

        impl ::core::ops::Neg for $typename {
            type Output = $typename;

            #[inline(always)]
            fn neg(self) -> $typename {
                let mut r = self;
                r.set_neg();
                r
            }
        }

        impl ::core::ops::Neg for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn neg(self) -> $typename {
                let mut r = *self;
                r.set_neg();
                r
            }
        }

        impl ::core::ops::Sub<$typename> for $typename {
            type Output = $typename;

            #[inline(always)]
            fn sub(self, other: $typename) -> $typename {
                let mut r = self;
                r.set_sub(&other);
                r
            }
        }

        impl ::core::ops::Sub<&$typename> for $typename {
            type Output = $typename;

            #[inline(always)]
            fn sub(self, other: &$typename) -> $typename {
                let mut r = self;
                r.set_sub(other);
                r
            }
        }

        impl ::core::ops::Sub<$typename> for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn sub(self, other: $typename) -> $typename {
                let mut r = *self;
                r.set_sub(&other);
                r
            }
        }

        impl ::core::ops::Sub<&$typename> for &$typename {
            type Output = $typename;

            #[inline(always)]
            fn sub(self, other: &$typename) -> $typename {
                let mut r = *self;
                r.set_sub(other);
                r
            }
        }

        impl ::core::ops::SubAssign<$typename> for $typename {
            #[inline(always)]
            fn sub_assign(&mut self, other: $typename) {
                self.set_sub(&other);
            }
        }

        impl ::core::ops::SubAssign<&$typename> for $typename {
            #[inline(always)]
            fn sub_assign(&mut self, other: &$typename) {
                self.set_sub(other);
            }
        }

        impl $crate::traits::Fq for $typename {
            // Reexport constants for base field Trait
            const N: usize = Self::N;
            const ENCODED_LENGTH: usize = Self::ENCODED_LENGTH;
            const ZERO: Self = Self::ZERO;
            const ONE: Self = Self::ONE;
            const TWO: Self = Self::TWO;
            const THREE: Self = Self::THREE;
            const FOUR: Self = Self::FOUR;
            const MINUS_ONE: Self = Self::MINUS_ONE;

            fn is_zero(self) -> u32 {
                self.is_zero()
            }

            fn equals(self, rhs: &Self) -> u32 {
                self.equals(rhs)
            }

            fn set_neg(&mut self) {
                self.set_neg()
            }

            fn set_half(&mut self) {
                self.set_half()
            }

            fn set_mul2(&mut self) {
                self.set_mul2()
            }

            fn set_mul3(&mut self) {
                self.set_mul3()
            }

            fn set_mul4(&mut self) {
                self.set_mul4()
            }

            fn set_mul8(&mut self) {
                self.set_mul8()
            }

            fn half(self) -> Self {
                self.half()
            }
            fn mul2(self) -> Self {
                self.mul2()
            }
            fn mul3(self) -> Self {
                self.mul3()
            }
            fn mul4(self) -> Self {
                self.mul4()
            }
            fn mul8(self) -> Self {
                self.mul8()
            }

            fn set_mul_small(&mut self, k: i32) {
                self.set_mul_small(k)
            }
            fn set_square(&mut self) {
                self.set_square()
            }
            fn set_n_square(&mut self, n: u32) {
                self.set_n_square(n)
            }
            fn set_invert(&mut self) {
                self.set_invert()
            }
            fn mul_small(self, k: i32) -> Self {
                self.mul_small(k)
            }
            fn square(self) -> Self {
                self.square()
            }
            fn n_square(self, n: u32) -> Self {
                self.n_square(n)
            }
            fn invert(self) -> Self {
                self.invert()
            }
            fn set_sqrt(&mut self) -> u32 {
                self.set_sqrt()
            }
            fn sqrt(self) -> (Self, u32) {
                self.sqrt()
            }
            fn legendre(self) -> i32 {
                self.legendre()
            }
            fn is_square(self) -> u32 {
                self.is_square()
            }
            fn batch_invert(xx: &mut [Self]) {
                <$typename>::batch_invert(xx)
            }

            fn set_select(&mut self, a: &Self, b: &Self, ctl: u32) {
                self.set_select(a, b, ctl)
            }
            fn set_cond(&mut self, rhs: &Self, ctl: u32) {
                self.set_cond(rhs, ctl)
            }
            fn set_cond_neg(&mut self, ctl: u32) {
                self.set_cond_neg(ctl)
            }
            fn select(a: &Self, b: &Self, ctl: u32) -> Self {
                <$typename>::select(a, b, ctl)
            }
            fn cond_swap(a: &mut Self, b: &mut Self, ctl: u32) {
                <$typename>::cond_swap(a, b, ctl)
            }

            fn encode(self) -> [u8; Self::ENCODED_LENGTH] {
                self.encode()
            }
            fn decode(buf: &[u8]) -> (Self, u32) {
                <$typename>::decode(buf)
            }
            fn decode_reduce(buf: &[u8]) -> Self {
                <$typename>::decode_reduce(buf)
            }

            fn hashcode(self) -> u64 {
                self.hashcode()
            }
        }

        impl $crate::traits::FqRoots for $typename {
            fn set_fourth_root(&mut self) -> u32 {
                self.set_fourth_root()
            }
            fn fourth_root(self) -> (Self, u32) {
                self.fourth_root()
            }
        }

        impl $crate::traits::FqExp for $typename {
            fn set_pow(&mut self, e: &[u8], ebitlen: usize) {
                self.set_pow(e, ebitlen)
            }
            fn set_pow_ext(&mut self, e: &[u8], eoff: usize, ebitlen: usize) {
                self.set_pow_ext(e, eoff, ebitlen)
            }
            fn set_pow_u64(&mut self, e: u64, ebitlen: usize) {
                self.set_pow_u64(e, ebitlen)
            }
            fn set_pow_u64_vartime(&mut self, e: u64) {
                self.set_pow_u64_vartime(e)
            }
            fn set_pow_pubexp(&mut self, e: &[u64; Self::N]) {
                self.set_pow_pubexp(e)
            }
            fn pow(self, e: &[u8], ebitlen: usize) -> Self {
                self.pow(e, ebitlen)
            }
            fn pow_ext(self, e: &[u8], eoff: usize, ebitlen: usize) -> Self {
                self.pow_ext(e, eoff, ebitlen)
            }
            fn pow_u64(self, e: u64, ebitlen: usize) -> Self {
                self.pow_u64(e, ebitlen)
            }
            fn pow_u64_vartime(self, e: u64) -> Self {
                self.pow_u64_vartime(e)
            }
            fn pow_pubexp(self, e: &[u64; Self::N]) -> Self {
                self.pow_pubexp(e)
            }
        }

        impl $crate::traits::FqRnd for $typename {
            fn set_rand<R: ::rand_core::CryptoRng + ::rand_core::RngCore>(&mut self, rng: &mut R) {
                self.set_rand(rng)
            }
            fn rand<R: ::rand_core::CryptoRng + ::rand_core::RngCore>(rng: &mut R) -> Self {
                <$typename>::rand(rng)
            }
        }
    };
} // End of macro: define_fp_core
