From d8d235fb485710886b8cbda5e6d85375e9933d29 Mon Sep 17 00:00:00 2001 From: Henry de Valence Date: Fri, 19 Jan 2018 13:49:10 -0800 Subject: [PATCH] Move pow2k into the backends and use it to implement square() --- src/backend/u32/field.rs | 10 +++ src/backend/u64/field.rs | 153 ++++++++++++++++++++++++++------------- src/field.rs | 9 --- 3 files changed, 112 insertions(+), 60 deletions(-) diff --git a/src/backend/u32/field.rs b/src/backend/u32/field.rs index a7a6419..fef4615 100644 --- a/src/backend/u32/field.rs +++ b/src/backend/u32/field.rs @@ -264,6 +264,16 @@ impl FieldElement32 { ]) } + /// Given `k > 0`, return `self^(2^k)`. + pub fn pow2k(&self, k: u32) -> FieldElement32 { + debug_assert!( k > 0 ); + let mut z = self.square(); + for _ in 1..k { + z = z.square(); + } + z + } + /// Given unreduced coefficients `z[0], ..., z[9]` of any size, /// carry and reduce them mod p to obtain a `FieldElement32` /// whose coefficients have excess `b < 0.007`. diff --git a/src/backend/u64/field.rs b/src/backend/u64/field.rs index 70d7fa8..1239491 100644 --- a/src/backend/u64/field.rs +++ b/src/backend/u64/field.rs @@ -337,73 +337,124 @@ impl FieldElement64 { s } - #[inline(always)] - fn square_inner(&self) -> [u64; 5] { + /// Given `k > 0`, return `self^(2^k)`. + pub fn pow2k(&self, mut k: u32) -> FieldElement64 { + + debug_assert!( k > 0 ); + /// Multiply two 64-bit integers with 128 bits of output. #[inline(always)] fn m(x: u64, y: u64) -> u128 { (x as u128) * (y as u128) } - // Alias self, _rhs for more readable formulas - let a: &[u64; 5] = &self.0; + let mut a: [u64; 5] = self.0; - // Precomputation: 64-bit multiply by 19 - let a3_19 = 19 * a[3]; - let a4_19 = 19 * a[4]; + loop { + // Precondition: assume input limbs a[i] are bounded as + // + // a[i] < 2^(51 + b) + // + // where b is a real parameter measuring the "bit excess" of the limbs. - // Multiply to get 128-bit coefficients of output - let c0: u128 = m(a[0], a[0]) + 2*( m(a[1], a4_19) + m(a[2], a3_19) ); - let mut c1: u128 = m(a[3], a3_19) + 2*( m(a[0], a[1]) + m(a[2], a4_19) ); - let mut c2: u128 = m(a[1], a[1]) + 2*( m(a[0], a[2]) + m(a[4], a3_19) ); - let mut c3: u128 = m(a[4], a4_19) + 2*( m(a[0], a[3]) + m(a[1], a[2]) ); - let mut c4: u128 = m(a[2], a[2]) + 2*( m(a[0], a[4]) + m(a[1], a[3]) ); + // Precomputation: 64-bit multiply by 19. + // + // This fits into a u64 whenever 51 + b + lg(19) < 64. + // + // Since 51 + b + lg(19) < 51 + 4.25 + b + // = 55.25 + b, + // this fits if b < 8.75. + let a3_19 = 19 * a[3]; + let a4_19 = 19 * a[4]; - // Same bound as in multiply: - // c[i] < 2^2b * (1+i + (4-i)*19) < 2^(2b + lg(1+4*19)) < 2^(2b + 6.27) - // where b is the bitlength of the input limbs. - // - // The carry (c[i] >> 51) fits into a u64 iff 2b+6.27 < 64+51 iff b <= 54. - // After the first carry pass, all c[i] fit into u64. - debug_assert!(a[0] < (1 << 54)); - debug_assert!(a[1] < (1 << 54)); - debug_assert!(a[2] < (1 << 54)); - debug_assert!(a[3] < (1 << 54)); - debug_assert!(a[4] < (1 << 54)); + // Multiply to get 128-bit coefficients of output. + // + // The 128-bit multiplications by 2 turn into 1 slr + 1 slrd each, + // which doesn't seem any better or worse than doing them as precomputations + // on the 64-bit inputs. + let c0: u128 = m(a[0], a[0]) + 2*( m(a[1], a4_19) + m(a[2], a3_19) ); + let mut c1: u128 = m(a[3], a3_19) + 2*( m(a[0], a[1]) + m(a[2], a4_19) ); + let mut c2: u128 = m(a[1], a[1]) + 2*( m(a[0], a[2]) + m(a[4], a3_19) ); + let mut c3: u128 = m(a[4], a4_19) + 2*( m(a[0], a[3]) + m(a[1], a[2]) ); + let mut c4: u128 = m(a[2], a[2]) + 2*( m(a[0], a[4]) + m(a[1], a[3]) ); - // The 128-bit output limbs are stored in two 64-bit registers (low/high part). - // By rebinding the names after carrying, we free the upper registers for reuse. - let low_51_bit_mask = (1u64 << 51) - 1; - c1 += (c0 >> 51) as u128; - let mut c0: u64 = (c0 as u64) & low_51_bit_mask; - c2 += (c1 >> 51) as u128; - let c1: u64 = (c1 as u64) & low_51_bit_mask; - c3 += (c2 >> 51) as u128; - let c2: u64 = (c2 as u64) & low_51_bit_mask; - c4 += (c3 >> 51) as u128; - let c3: u64 = (c3 as u64) & low_51_bit_mask; - c0 += ((c4 >> 51) as u64) * 19; - let c4: u64 = (c4 as u64) & low_51_bit_mask; + // Same bound as in multiply: + // c[i] < 2^(102 + 2*b) * (1+i + (4-i)*19) + // < 2^(102 + lg(1 + 4*19) + 2*b) + // < 2^(108.27 + 2*b) + // + // The carry (c[i] >> 51) fits into a u64 when + // 108.27 + 2*b - 51 < 64 + // 2*b < 6.73 + // b < 3.365. + // + // So we require b < 3 to ensure this fits. + debug_assert!(a[0] < (1 << 54)); + debug_assert!(a[1] < (1 << 54)); + debug_assert!(a[2] < (1 << 54)); + debug_assert!(a[3] < (1 << 54)); + debug_assert!(a[4] < (1 << 54)); - // Now c_i all fit into u64, but are not yet bounded by 2^51. - [c0,c1,c2,c3,c4] + const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1; + + // Casting to u64 and back tells the compiler that the carry is bounded by 2^64, so + // that the addition is a u128 + u64 rather than u128 + u128. + c1 += ((c0 >> 51) as u64) as u128; + a[0] = (c0 as u64) & LOW_51_BIT_MASK; + + c2 += ((c1 >> 51) as u64) as u128; + a[1] = (c1 as u64) & LOW_51_BIT_MASK; + + c3 += ((c2 >> 51) as u64) as u128; + a[2] = (c2 as u64) & LOW_51_BIT_MASK; + + c4 += ((c3 >> 51) as u64) as u128; + a[3] = (c3 as u64) & LOW_51_BIT_MASK; + + let carry: u64 = (c4 >> 51) as u64; + a[4] = (c4 as u64) & LOW_51_BIT_MASK; + + // To see that this does not overflow, we need a[0] + carry * 19 < 2^64. + // + // c4 < a2^2 + 2*a0*a4 + 2*a1*a3 + (carry from c3) + // < 2^(102 + 2*b + lg(5)) + 2^64. + // + // When b < 3 we get + // + // c4 < 2^110.33 so that carry < 2^59.33 + // + // so that + // + // a[0] + carry * 19 < 2^51 + 19 * 2^59.33 < 2^63.58 + // + // and there is no overflow. + a[0] = a[0] + carry * 19; + + // Now a[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + eps). + a[1] += a[0] >> 51; + a[0] &= LOW_51_BIT_MASK; + + // Now all a[i] < 2^(51 + eps) and a = self^(2^k). + + k = k - 1; + if k == 0 { + break; + } + } + + FieldElement64(a) } /// Returns the square of this field element. pub fn square(&self) -> FieldElement64 { - FieldElement64::reduce(self.square_inner()) + self.pow2k(1) } /// Returns 2 times the square of this field element. pub fn square2(&self) -> FieldElement64 { - let mut limbs = self.square_inner(); - // For this to work, need to have 1 extra bit of headroom after carry - // --> max 53 bit inputs, not 54 - // - // XXX check that this is correct; I think it isn't -- hdevalence - limbs[0] *= 2; - limbs[1] *= 2; - limbs[2] *= 2; - limbs[3] *= 2; - limbs[4] *= 2; - FieldElement64::reduce(limbs) + let mut square = self.pow2k(1); + for i in 0..5 { + square.0[i] *= 2; + } + + square } } diff --git a/src/field.rs b/src/field.rs index e6f2587..756423e 100644 --- a/src/field.rs +++ b/src/field.rs @@ -124,15 +124,6 @@ impl FieldElement { byte_is_nonzero(x) } - #[inline] - #[allow(dead_code)] - /// Requires k > 0; raise self to the 2^(2^k)-th power. - fn pow2k(&self, k: u32) -> FieldElement { - let mut z = self.square(); - for _ in 1..k { z = z.square(); } - z - } - /// Compute (self^(2^250-1), self^11), used as a helper function /// within invert() and pow22523(). ///