diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs index d688fca..7c5ff70 100644 --- a/src/backend/avx2/edwards.rs +++ b/src/backend/avx2/edwards.rs @@ -86,7 +86,7 @@ impl ExtendedPoint { // Set tmp1 = tmp0^2, negating the D values tmp1 = tmp0.square_and_negate_D(); - // Now tmp1 = (S1 S2 S3 -S4) + // Now tmp1 = (S1 S2 S3 -S4) with b < 0.007 // See discussion of bounds in the module-level documentation. // We want to compute @@ -105,23 +105,27 @@ impl ExtendedPoint { let S_1 = tmp1.shuffle(Shuffle::AAAA); let S_2 = tmp1.shuffle(Shuffle::BBBB); - // tmp0 = (0, 0, 2S_3, 0) tmp0 = zero.blend(tmp1 + tmp1, Lanes::C); - // tmp0 = (0, 0, 2S_3, -S_4) + // tmp0 = (0, 0, 2S_3, 0) tmp0 = tmp0.blend(tmp1, Lanes::D); - // tmp0 = (S_1, S_1, S_1 + 2S_3, S_1 - S_4) + // tmp0 = (0, 0, 2S_3, -S_4) tmp0 = tmp0 + S_1; - // tmp0 = (S_1 + S_2, S_1, S_1 + 2S_3, S_1 + S_2 - S_4) + // tmp0 = ( S_1, S_1, S_1 + 2S_3, S_1 - S_4) tmp0 = tmp0 + zero.blend(S_2, Lanes::AD); - // tmp0 = (S_1 + S_2, S_1 - S_2, S_1 - S_2 + 2S_3, S_1 + S_2 - S_4) + // tmp0 = (S_1 + S_2, S_1, S_1 + 2S_3, S_1 + S_2 - S_4) tmp0 = tmp0 + zero.blend(S_2.negate_lazy(), Lanes::BC); + // tmp0 = (S_1 + S_2, S_1 - S_2, S_1 - S_2 + 2S_3, S_1 + S_2 - S_4) + // b < ( 1.01, 1.6, 2.33, 1.6) // Now tmp0 = (S_5, S_6, S_8, S_9) - // Set tmp1 = (S_9, S_6, S_6, S_9) + // Set tmp1 = ( S_9, S_6, S_6, S_9) + // b < ( 1.6, 1.6, 1.6, 1.6) tmp1 = tmp0.shuffle(Shuffle::DBBD); - // Set tmp1 = (S_8, S_5, S_8, S_5) + // Set tmp1 = ( S_8, S_5, S_8, S_5) + // b < (2.33, 1.01, 2.33, 1.01) tmp0 = tmp0.shuffle(Shuffle::CACA); + // Bounds on (tmp0, tmp1) are (2.33, 1.6) < (2.5, 1.75). ExtendedPoint(&tmp0 * &tmp1) } @@ -142,15 +146,16 @@ impl From for CachedPoint { fn from(P: ExtendedPoint) -> CachedPoint { let mut x = P.0; - // x = (S2 S3 Z2 T2) x = x.blend(x.diff_sum(), Lanes::AB); + // x = (X1 - Y1, X2 + Y2, Z2, T2) = (S2 S3 Z2 T2) - // x = (121666*S2 121666*S3 2*121666*Z2 2*121665*T2) x = x * (121666, 121666, 2*121666, 2*121665); + // x = (121666*S2 121666*S3 2*121666*Z2 2*121665*T2) - // x = (121666*S2 121666*S3 2*121666*Z2 -2*121665*T2) x = x.blend(-x, Lanes::D); + // x = (121666*S2 121666*S3 2*121666*Z2 -2*121665*T2) + // The coefficients of the output are bounded with b < 0.007. CachedPoint(x) } } @@ -175,7 +180,12 @@ impl ConditionallyAssignable for CachedPoint { impl<'a> Neg for &'a CachedPoint { type Output = CachedPoint; - + /// Lazily negate the point. + /// + /// # Warning + /// + /// Because this method does not perform a reduction, it is not + /// safe to repeatedly negate a point. fn neg(self) -> CachedPoint { let swapped = self.0.shuffle(Shuffle::BACD); CachedPoint(swapped.blend(swapped.negate_lazy(), Lanes::D)) @@ -185,28 +195,30 @@ impl<'a> Neg for &'a CachedPoint { impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { type Output = ExtendedPoint; - /// Uses a slight tweak of the parallel unified formulas of HWCD'08 + /// Add an `ExtendedPoint` and a `CachedPoint`. fn add(self, other: &'b CachedPoint) -> ExtendedPoint { let mut tmp = self.0; - // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) tmp = tmp.blend(tmp.diff_sum(), Lanes::AB); + // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) with b < 1.6 - // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) + // (tmp, other) bounded with b < (1.6, 1.0) < (2.5, 1.75). tmp = &tmp * &other.0; + // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) - // tmp = (S8 S9 S11 S10) tmp = tmp.shuffle(Shuffle::ABDC); + // tmp = (S8 S9 S11 S10) - // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) tmp = tmp.diff_sum(); + // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) - // set t0 = (S12 S15 S15 S12) let t0 = tmp.shuffle(Shuffle::ADDA); - // set t1 = (S14 S13 S14 S13) + // t0 = (S12 S15 S15 S12) let t1 = tmp.shuffle(Shuffle::CBCB); + // t1 = (S14 S13 S14 S13) - // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) + // All coefficients of t0, t1 are bounded with b < 1.6. + // Return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) ExtendedPoint(&t0 * &t1) } } @@ -216,8 +228,9 @@ impl<'a, 'b> Sub<&'b CachedPoint> for &'a ExtendedPoint { /// Implement subtraction by negating the point and adding. /// - /// Empirically, this seems about the same cost as a custom subtraction impl (maybe because the - /// benefit is cancelled by increased code size?) + /// Empirically, this seems about the same cost as a custom + /// subtraction impl (maybe because the benefit is cancelled by + /// increased code size?) fn sub(self, other: &'b CachedPoint) -> ExtendedPoint { self + &(-other) } diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs index 297ba34..1484d75 100644 --- a/src/backend/avx2/field.rs +++ b/src/backend/avx2/field.rs @@ -112,7 +112,7 @@ pub enum Shuffle { ABDC, } -/// A vector of four `FieldElements`, implemented using AVX2. +/// A vector of four field elements, in an AVX2-friendly format. #[derive(Clone, Copy, Debug)] pub struct FieldElement32x4(pub(crate) [u32x8; 5]); @@ -258,14 +258,21 @@ impl FieldElement32x4 { ]) } + /// Construct a vector of zeros. pub fn zero() -> FieldElement32x4 { FieldElement32x4([u32x8::splat(0); 5]) } + /// Convenience wrapper around `new(x,x,x,x)`. pub fn splat(x: &FieldElement64) -> FieldElement32x4 { FieldElement32x4::new(x, x, x, x) } + /// Create a `FieldElement32x4` from four `FieldElement64`s. + /// + /// # Postconditions + /// + /// The resulting `FieldElement32x4` is bounded with \\( b < 0.0002 \\). pub fn new( x0: &FieldElement64, x1: &FieldElement64, @@ -296,11 +303,18 @@ impl FieldElement32x4 { /// Given \\((A,B,C,D)\\), compute \\((-A,-B,-C,-D)\\), without /// performing a reduction. /// - /// Input limbs must be less than the limbs of \\(2p\\), i.e., freshly reduced. + /// # Preconditions /// - /// The output limbs are bounded by \\(2p\\). + /// The coefficients of `self` must be bounded with \\( b < 0.999 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 1 \\). #[inline] pub fn negate_lazy(&self) -> FieldElement32x4 { + // The limbs of self are bounded with b < 0.999, while the + // smallest limb of 2*p is 67108845 > 2^{26+0.9999}, so + // underflows are not possible. FieldElement32x4([ P_TIMES_2_LO - self.0[0], P_TIMES_2_HI - self.0[1], @@ -311,17 +325,29 @@ impl FieldElement32x4 { } /// Given `self = (A,B,C,D)`, compute `(B - A, B + A, D - C, D + C)`. + /// + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 0.01 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 1.6 \\). #[inline] pub fn diff_sum(&self) -> FieldElement32x4 { // tmp1 = (B, A, D, C) let tmp1 = self.shuffle(Shuffle::BADC); // tmp2 = (-A, B, -C, D) let tmp2 = self.blend(self.negate_lazy(), Lanes::AC); - // (B - A, B + A, D - C, D + C) + // (B - A, B + A, D - C, D + C) bounded with b < 1.6 tmp1 + tmp2 } - /// Compute the reduced form of the field elements. + /// Reduce this vector of field elements \\(\mathrm{mod} p\\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.0002 \\). #[inline] pub fn reduce(&self) -> FieldElement32x4 { let shifts = i32x8::new(26, 26, 25, 25, 26, 26, 25, 25); @@ -336,22 +362,36 @@ impl FieldElement32x4 { (1 << 25) - 1, ); - // Compute the carryout of v. - let carry = |v: u32x8| -> u32x8 { + // Let c(x) denote the carryout of the coefficient x. + // + // Given ( x0, y0, x1, y1, z0, w0, z1, w1), + // compute (c(x1), c(y1), c(x0), c(y0), c(z1), c(w1), c(z0), c(w0)). + // + // The carryouts are bounded by 2^(32 - 25) = 2^7. + let rotated_carryout = |v: u32x8| -> u32x8 { unsafe { use core::arch::x86_64::_mm256_srlv_epi32; - _mm256_srlv_epi32(v.into_bits(), shifts.into_bits()).into_bits() - } - }; - - // Swap adjacent 32-bit lanes. - let swap_lanes = |v: u32x8| -> u32x8 { - unsafe { use core::arch::x86_64::_mm256_shuffle_epi32; - _mm256_shuffle_epi32(v.into_bits(), 0b01_00_11_10).into_bits() + + let c = _mm256_srlv_epi32(v.into_bits(), shifts.into_bits()); + _mm256_shuffle_epi32(c, 0b01_00_11_10).into_bits() } }; + // Combine (lo, lo, lo, lo, lo, lo, lo, lo) + // with (hi, hi, hi, hi, hi, hi, hi, hi) + // to (lo, lo, hi, hi, lo, lo, hi, hi) + // + // This allows combining carryouts, e.g., + // + // lo (c(x1), c(y1), c(x0), c(y0), c(z1), c(w1), c(z0), c(w0)) + // hi (c(x3), c(y3), c(x2), c(y2), c(z3), c(w3), c(z2), c(w2)) + // -> (c(x1), c(y1), c(x2), c(y2), c(z1), c(w1), c(z2), c(w2)) + // + // which is exactly the vector of carryins for + // + // ( x2, y2, x3, y3, z2, w2, z3, w3). + // let combine = |v_lo: u32x8, v_hi: u32x8| -> u32x8 { unsafe { use core::arch::x86_64::_mm256_blend_epi32; @@ -361,34 +401,62 @@ impl FieldElement32x4 { let mut v = self.0; - let c10 = swap_lanes(carry(v[0])); + let c10 = rotated_carryout(v[0]); v[0] = (v[0] & masks) + combine(u32x8::splat(0), c10); - let c32 = swap_lanes(carry(v[1])); + + let c32 = rotated_carryout(v[1]); v[1] = (v[1] & masks) + combine(c10, c32); - let c54 = swap_lanes(carry(v[2])); + + let c54 = rotated_carryout(v[2]); v[2] = (v[2] & masks) + combine(c32, c54); - let c76 = swap_lanes(carry(v[3])); + + let c76 = rotated_carryout(v[3]); v[3] = (v[3] & masks) + combine(c54, c76); - let c98 = swap_lanes(carry(v[4])); + + let c98 = rotated_carryout(v[4]); v[4] = (v[4] & masks) + combine(c76, c98); - // Still need to account for c9 - // c98 = (c9, c9, c8, c8, c9, c9, c8, c8) - // - let c9_19: u32x8; - unsafe { + let c9_19: u32x8 = unsafe { use core::arch::x86_64::_mm256_mul_epu32; use core::arch::x86_64::_mm256_shuffle_epi32; - let c9_spread = _mm256_shuffle_epi32(c98.into_bits(), 0b11_01_10_00); - let c9_19_spread = _mm256_mul_epu32(c9_spread, u64x4::splat(19).into_bits()); - c9_19 = _mm256_shuffle_epi32(c9_19_spread, 0b11_01_10_00).into_bits(); - } + // Need to rearrange c98, since vpmuludq uses the low + // 32-bits of each 64-bit lane to compute the product: + // + // c98 = (c(x9), c(y9), c(x8), c(y8), c(z9), c(w9), c(z8), c(w8)); + // c9_spread = (c(x9), c(x8), c(y9), c(y8), c(z9), c(z8), c(w9), c(w8)). + let c9_spread = _mm256_shuffle_epi32(c98.into_bits(), 0b11_01_10_00); + + // Since the carryouts are bounded by 2^7, their products with 19 + // are bounded by 2^11.25. This means that + // + // c9_19_spread = (19*c(x9), 0, 19*c(y9), 0, 19*c(z9), 0, 19*c(w9), 0). + let c9_19_spread = _mm256_mul_epu32(c9_spread, u64x4::splat(19).into_bits()); + + // Unshuffle: + // c9_19 = (19*c(x9), 19*c(y9), 0, 0, 19*c(z9), 19*c(w9), 0, 0). + _mm256_shuffle_epi32(c9_19_spread, 0b11_01_10_00).into_bits() + }; + + // Add the final carryin. v[0] = v[0] + c9_19; + // Each output coefficient has exactly one carryin, which is + // bounded by 2^11.25, so they are bounded as + // + // c_even < 2^26 + 2^11.25 < 26.00006 < 2^{26+b} + // c_odd < 2^25 + 2^11.25 < 25.0001 < 2^{25+b} + // + // where b = 0.0002. FieldElement32x4(v) } + /// Given an array of wide coefficients, reduce them to a `FieldElement32x4`. + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). + #[inline] fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { // These aren't const because splat isn't a const fn let LOW_25_BITS: u64x4 = u64x4::splat((1 << 25) - 1); @@ -442,8 +510,13 @@ impl FieldElement32x4 { z[1] = z[1] + c1; // z1 < 2^25 + 2^17.25 < 2^25.0067 carry(&mut z, 0); // z0 < 2^26, z1 < 2^25.0067 + 2^4.33 = 2^25.007 - // Now repack the [u64x4; 10] into a FieldElement32x4 - + // The output coefficients are bounded with + // + // b = 0.007 for z[1] + // b = 0.0004 for z[5] + // b = 0 for other z[i]. + // + // So the packed result is bounded with b = 0.007. FieldElement32x4([ repack_pair(z[0].into_bits(), z[1].into_bits()), repack_pair(z[2].into_bits(), z[3].into_bits()), @@ -453,13 +526,15 @@ impl FieldElement32x4 { ]) } - /// Square this field element, then conditionally negate according - /// to `neg_mask`. This parameter is hardcoded as `neg_mask = - /// D_LANES64` to negate the \\( D \\) value. + /// Square this field element, and negate the result's \\(D\\) value. /// - /// # Precondition + /// # Preconditions /// - /// Limbs must be bounded by bit-excess \\( b < 2.0 \\). + /// The coefficients of `self` must be bounded with \\( b < 1.5 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). pub fn square_and_negate_D(&self) -> FieldElement32x4 { #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { @@ -552,12 +627,18 @@ impl FieldElement32x4 { impl Neg for FieldElement32x4 { type Output = FieldElement32x4; - /// Given \\((A,B,C,D)\\), compute \\((-A,-B,-C,-D)\\), and - /// perform a reduction. + /// Negate this field element, performing a reduction. /// - /// The input limbs can be any size. + /// If the coefficients are known to be small, use `negate_lazy` + /// to avoid performing a reduction. /// - /// The output limbs are freshly reduced. + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 4.0 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.0002 \\). #[inline] fn neg(self) -> FieldElement32x4 { FieldElement32x4([ @@ -572,6 +653,7 @@ impl Neg for FieldElement32x4 { impl Add for FieldElement32x4 { type Output = FieldElement32x4; + /// Add two `FieldElement32x4`s, without performing a reduction. #[inline] fn add(self, rhs: FieldElement32x4) -> FieldElement32x4 { FieldElement32x4([ @@ -586,6 +668,11 @@ impl Add for FieldElement32x4 { impl Mul<(u32, u32, u32, u32)> for FieldElement32x4 { type Output = FieldElement32x4; + /// Perform a multiplication by a vector of small constants. + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). #[inline] fn mul(self, scalars: (u32, u32, u32, u32)) -> FieldElement32x4 { unsafe { @@ -617,7 +704,19 @@ impl Mul<(u32, u32, u32, u32)> for FieldElement32x4 { impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { type Output = FieldElement32x4; - fn mul(self, _rhs: &'b FieldElement32x4) -> FieldElement32x4 { + /// Multiply `self` by `rhs`. + /// + /// # Preconditions + /// + /// The coefficients of `self` must be bounded with \\( b < 2.5 \\). + /// + /// The coefficients of `rhs` must be bounded with \\( b < 1.75 \\). + /// + /// # Postconditions + /// + /// The coefficients of the result are bounded with \\( b < 0.007 \\). + /// + fn mul(self, rhs: &'b FieldElement32x4) -> FieldElement32x4 { #[inline(always)] fn m(x: u32x8, y: u32x8) -> u64x4 { use core::arch::x86_64::_mm256_mul_epu32; @@ -636,11 +735,11 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let (x6, x7) = unpack_pair(self.0[3]); let (x8, x9) = unpack_pair(self.0[4]); - let (y0, y1) = unpack_pair(_rhs.0[0]); - let (y2, y3) = unpack_pair(_rhs.0[1]); - let (y4, y5) = unpack_pair(_rhs.0[2]); - let (y6, y7) = unpack_pair(_rhs.0[3]); - let (y8, y9) = unpack_pair(_rhs.0[4]); + let (y0, y1) = unpack_pair(rhs.0[0]); + let (y2, y3) = unpack_pair(rhs.0[1]); + let (y4, y5) = unpack_pair(rhs.0[2]); + let (y6, y7) = unpack_pair(rhs.0[3]); + let (y8, y9) = unpack_pair(rhs.0[4]); let v19 = u32x8::new(19, 0, 19, 0, 19, 0, 19, 0); @@ -648,9 +747,9 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let y2_19 = m_lo(v19, y2); // iff 26 + b + lg(19) < 32 let y3_19 = m_lo(v19, y3); // if b < 32 - 26 - 4.248 = 1.752 let y4_19 = m_lo(v19, y4); - let y5_19 = m_lo(v19, y5); // below, b<2.5: this is a bottleneck, - let y6_19 = m_lo(v19, y6); // could be avoided by promoting to - let y7_19 = m_lo(v19, y7); // u64 here instead of in m() + let y5_19 = m_lo(v19, y5); + let y6_19 = m_lo(v19, y6); + let y7_19 = m_lo(v19, y7); let y8_19 = m_lo(v19, y8); let y9_19 = m_lo(v19, y9); @@ -671,6 +770,44 @@ impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { let z8 = m(x0,y8) + m(x1_2,y7) + m(x2,y6) + m(x3_2,y5) + m(x4,y4) + m(x5_2,y3) + m(x6,y2) + m(x7_2,y1) + m(x8,y0) + m(x9_2,y9_19); let z9 = m(x0,y9) + m(x1,y8) + m(x2,y7) + m(x3,y6) + m(x4,y5) + m(x5,y4) + m(x6,y3) + m(x7,y2) + m(x8,y1) + m(x9,y0); + // The bounds on z[i] are the same as in the serial 32-bit code + // and the comment below is copied from there: + + // How big is the contribution to z[i+j] from x[i], y[j]? + // + // Using the bounds above, we get: + // + // i even, j even: x[i]*y[j] < 2^(26+b)*2^(26+b) = 2*2^(51+2*b) + // i odd, j even: x[i]*y[j] < 2^(25+b)*2^(26+b) = 1*2^(51+2*b) + // i even, j odd: x[i]*y[j] < 2^(26+b)*2^(25+b) = 1*2^(51+2*b) + // i odd, j odd: 2*x[i]*y[j] < 2*2^(25+b)*2^(25+b) = 1*2^(51+2*b) + // + // We perform inline reduction mod p by replacing 2^255 by 19 + // (since 2^255 - 19 = 0 mod p). This adds a factor of 19, so + // we get the bounds (z0 is the biggest one, but calculated for + // posterity here in case finer estimation is needed later): + // + // z0 < ( 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 249*2^(51 + 2*b) + // z1 < ( 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 154*2^(51 + 2*b) + // z2 < ( 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 195*2^(51 + 2*b) + // z3 < ( 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 118*2^(51 + 2*b) + // z4 < ( 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 141*2^(51 + 2*b) + // z5 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 82*2^(51 + 2*b) + // z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 87*2^(51 + 2*b) + // z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 )*2^(51 + 2b) = 46*2^(51 + 2*b) + // z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 )*2^(51 + 2b) = 33*2^(51 + 2*b) + // z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 )*2^(51 + 2b) = 10*2^(51 + 2*b) + // + // So z[0] fits into a u64 if 51 + 2*b + lg(249) < 64 + // if b < 2.5. + + // In fact this bound is slightly sloppy, since it treats both + // inputs x and y as being bounded by the same parameter b, + // while they are in fact bounded by b_x and b_y, and we + // already require that b_y < 1.75 in order to fit the + // multiplications by 19 into a u32. The tighter bound on b_y + // means we could get a tighter bound on the outputs, or a + // looser bound on b_x. FieldElement32x4::reduce64([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9]) } }