Move pow2k into the backends and use it to implement square()

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`.

src/backend/u64/field.rs

@@ -337,73 +337,124 @@ impl FieldElement64 {
         s
     }
 
-    #[inline(always)]
+    /// Given `k > 0`, return `self^(2^k)`.
-    fn square_inner(&self) -> [u64; 5] {
+    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 mut a: [u64; 5] = self.0;
-        let a: &[u64; 5] = &self.0;
 
-        // Precomputation: 64-bit multiply by 19
+        loop {
-        let a3_19 = 19 * a[3];
+            // Precondition: assume input limbs a[i] are bounded as
-        let a4_19 = 19 * a[4];
+            //
+            // 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
+            // Precomputation: 64-bit multiply by 19.
-        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) );
+            // This fits into a u64 whenever 51 + b + lg(19) < 64.
-        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]) );
+            // Since 51 + b + lg(19) < 51 + 4.25 + b
-        let mut c4: u128 = m(a[2],  a[2]) + 2*( m(a[0],  a[4]) + m(a[1],  a[3]) );
+            //                       = 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:
+            // Multiply to get 128-bit coefficients of output.
-        //    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 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
-        // The carry (c[i] >> 51) fits into a u64 iff 2b+6.27 < 64+51 iff b <= 54.
+            // on the 64-bit inputs.
-        // After the first carry pass, all c[i] fit into u64.
+            let     c0: u128 = m(a[0],  a[0]) + 2*( m(a[1], a4_19) + m(a[2], a3_19) );
-        debug_assert!(a[0] < (1 << 54));
+            let mut c1: u128 = m(a[3], a3_19) + 2*( m(a[0],  a[1]) + m(a[2], a4_19) );
-        debug_assert!(a[1] < (1 << 54));
+            let mut c2: u128 = m(a[1],  a[1]) + 2*( m(a[0],  a[2]) + m(a[4], a3_19) );
-        debug_assert!(a[2] < (1 << 54));
+            let mut c3: u128 = m(a[4], a4_19) + 2*( m(a[0],  a[3]) + m(a[1],  a[2]) );
-        debug_assert!(a[3] < (1 << 54));
+            let mut c4: u128 = m(a[2],  a[2]) + 2*( m(a[0],  a[4]) + m(a[1],  a[3]) );
-        debug_assert!(a[4] < (1 << 54));
 
-        // The 128-bit output limbs are stored in two 64-bit registers (low/high part).
+            // Same bound as in multiply:
-        // By rebinding the names after carrying, we free the upper registers for reuse.
+            //    c[i] < 2^(102 + 2*b) * (1+i + (4-i)*19) 
-        let low_51_bit_mask = (1u64 << 51) - 1;
+            //         < 2^(102 + lg(1 + 4*19) + 2*b)
-        c1 +=  (c0 >> 51) as u128;
+            //         < 2^(108.27 + 2*b)
-        let mut c0: u64 = (c0 as u64) & low_51_bit_mask;
+            //
-        c2 +=  (c1 >> 51) as u128;
+            // The carry (c[i] >> 51) fits into a u64 when
-        let c1: u64 = (c1 as u64) & low_51_bit_mask;
+            //    108.27 + 2*b - 51 < 64
-        c3 +=  (c2 >> 51) as u128;
+            //    2*b < 6.73
-        let c2: u64 = (c2 as u64) & low_51_bit_mask;
+            //    b < 3.365.
-        c4 +=  (c3 >> 51) as u128;
+            //
-        let c3: u64 = (c3 as u64) & low_51_bit_mask;
+            // So we require b < 3 to ensure this fits.
-        c0 += ((c4 >> 51) as u64) * 19;
+            debug_assert!(a[0] < (1 << 54));
-        let c4: u64 = (c4 as u64) & low_51_bit_mask;
+            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.
+            const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1;
-        [c0,c1,c2,c3,c4]
+
+            // 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();
+        let mut square = self.pow2k(1);
-        // For this to work, need to have 1 extra bit of headroom after carry
+        for i in 0..5 {
-        // --> max 53 bit inputs, not 54
+            square.0[i] *= 2;
-        //
+        }
-        // XXX check that this is correct; I think it isn't -- hdevalence
+
-        limbs[0] *= 2;
+        square
-        limbs[1] *= 2;
-        limbs[2] *= 2;
-        limbs[3] *= 2;
-        limbs[4] *= 2;
-        FieldElement64::reduce(limbs)
     }
 }

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().
     ///