Merge remote-tracking branch 'dalek/feature/optimize-64-bit' into develop

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

@@ -78,8 +78,7 @@ impl<'a, 'b> Sub<&'b FieldElement64> for &'a FieldElement64 {
         // just bigger than _rhs and avoid having to do a reduction.
         //
         // Since we don't yet have type-level integers to do this, we
-        // have to add an explicit reduction call here, which is a
+        // have to add an explicit reduction call here.
-        // somewhat significant cost.
         FieldElement64::reduce([
             (self.0[0] + 36028797018963664u64) - _rhs.0[0],
             (self.0[1] + 36028797018963952u64) - _rhs.0[1],
@@ -109,7 +108,19 @@ impl<'a, 'b> Mul<&'b FieldElement64> for &'a FieldElement64 {
         let a: &[u64; 5] = &self.0;
         let b: &[u64; 5] = &_rhs.0;
 
-        // 64-bit precomputations to avoid 128-bit multiplications
+        // Precondition: assume input limbs a[i], b[i] are bounded as
+        //
+        // a[i], b[i] < 2^(51 + b)
+        //
+        // where b is a real parameter measuring the "bit excess" of the limbs.
+
+        // 64-bit precomputations to avoid 128-bit multiplications.
+        //
+        // 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 b1_19 = b[1] * 19;
         let b2_19 = b[2] * 19;
         let b3_19 = b[3] * 19;
@@ -122,34 +133,69 @@ impl<'a, 'b> Mul<&'b FieldElement64> for &'a FieldElement64 {
         let mut c3: u128 = m(a[3],b[0]) + m(a[2],b[1])  + m(a[1],b[2])  + m(a[0],b[3])  + m(a[4],b4_19);
         let mut c4: u128 = m(a[4],b[0]) + m(a[3],b[1])  + m(a[2],b[2])  + m(a[1],b[3])  + m(a[0],b[4]);
 
-        // Now c[i] < 2^2b * (1+i + (4-i)*19) < 2^(2b + lg(1+4*19)) < 2^(2b + 6.27)
+        // How big are the c[i]? We have
-        // where b is the bitlength of the input limbs.
+        //
-
+        //    c[i] < 2^(102 + 2*b) * (1+i + (4-i)*19) 
-        // The carry (c[i] >> 51) fits into a u64 iff 2b+6.27 < 64+51 iff b <= 54.
+        //         < 2^(102 + lg(1 + 4*19) + 2*b)
-        // After the first carry pass, all c[i] fit into u64.
+        //         < 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!(b[0] < (1 << 54));
         debug_assert!(a[1] < (1 << 54)); debug_assert!(b[1] < (1 << 54));
         debug_assert!(a[2] < (1 << 54)); debug_assert!(b[2] < (1 << 54));
         debug_assert!(a[3] < (1 << 54)); debug_assert!(b[3] < (1 << 54));
         debug_assert!(a[4] < (1 << 54)); debug_assert!(b[4] < (1 << 54));
 
-        // The 128-bit output limbs are stored in two 64-bit registers
+        // Casting to u64 and back tells the compiler that the carry is
-        // (low/high part).  By rebinding the names after carrying, we
+        // bounded by 2^64, so that the addition is a u128 + u64 rather
-        // inform LLVM that the values have shrunk, so it can
+        // than u128 + u128.
-        // efficiently allocate registers.
-        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;
 
-        FieldElement64::reduce([c0,c1,c2,c3,c4])
+        const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1;
+        let mut out = [0u64; 5];
+
+        c1 += ((c0 >> 51) as u64) as u128;
+        out[0] = (c0 as u64) & LOW_51_BIT_MASK;
+
+        c2 += ((c1 >> 51) as u64) as u128;
+        out[1] = (c1 as u64) & LOW_51_BIT_MASK;
+
+        c3 += ((c2 >> 51) as u64) as u128;
+        out[2] = (c2 as u64) & LOW_51_BIT_MASK;
+
+        c4 += ((c3 >> 51) as u64) as u128;
+        out[3] = (c3 as u64) & LOW_51_BIT_MASK;
+
+        let carry: u64 = (c4 >> 51) as u64;
+        out[4] = (c4 as u64) & LOW_51_BIT_MASK;
+
+        // To see that this does not overflow, we need out[0] + carry * 19 < 2^64.
+        //
+        // c4 < a0*b4 + a1*b3 + a2*b2 + a3*b1 + a4*b0 + (carry from c3)
+        //    < 5*(2^(51 + b) * 2^(51 + b)) + (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
+        //
+        // out[0] + carry * 19 < 2^51 + 19 * 2^59.33 < 2^63.58
+        //
+        // and there is no overflow.
+        out[0] = out[0] + carry * 19;
+
+        // Now out[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + epsilon).
+        out[1] += out[0] >> 51;
+        out[0] &= LOW_51_BIT_MASK;
+
+        // Now out[i] < 2^(51 + epsilon) for all i.
+        FieldElement64(out)
     }
 }
 
@@ -200,20 +246,39 @@ impl FieldElement64 {
         FieldElement64([2251799813685228, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247])
     }
 
-    /// Given 64-bit limbs, reduce to enforce the bound c_i < 2^51.
+    /// Given 64-bit input limbs, reduce to enforce the bound 2^(51 + epsilon).
     #[inline(always)]
     fn reduce(mut limbs: [u64; 5]) -> FieldElement64 {
-        let low_51_bit_mask = (1u64 << 51) - 1;
+        const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1;
-        limbs[1] +=  limbs[0] >> 51;
+
-        limbs[0] = limbs[0] & low_51_bit_mask;
+        // Since the input limbs are bounded by 2^64, the biggest
-        limbs[2] +=  limbs[1] >> 51;
+        // carry-out is bounded by 2^13.
-        limbs[1] = limbs[1] & low_51_bit_mask;
+        //
-        limbs[3] +=  limbs[2] >> 51;
+        // The biggest carry-in is c4 * 19, resulting in
-        limbs[2] = limbs[2] & low_51_bit_mask;
+        //
-        limbs[4] +=  limbs[3] >> 51;
+        // 2^51 + 19*2^13 < 2^51.0000000001
-        limbs[3] = limbs[3] & low_51_bit_mask;
+        //
-        limbs[0] += (limbs[4] >> 51) * 19;
+        // Because we don't need to canonicalize, only to reduce the
-        limbs[4] = limbs[4] & low_51_bit_mask;
+        // limb sizes, it's OK to do a "weak reduction", where we
+        // compute the carry-outs in parallel.
+        
+        let c0 = limbs[0] >> 51;
+        let c1 = limbs[1] >> 51;
+        let c2 = limbs[2] >> 51;
+        let c3 = limbs[3] >> 51;
+        let c4 = limbs[4] >> 51;
+        
+        limbs[0] &= LOW_51_BIT_MASK;
+        limbs[1] &= LOW_51_BIT_MASK;
+        limbs[2] &= LOW_51_BIT_MASK;
+        limbs[3] &= LOW_51_BIT_MASK;
+        limbs[4] &= LOW_51_BIT_MASK;
+        
+        limbs[0] += c4 * 19;
+        limbs[1] += c0;
+        limbs[2] += c1;
+        limbs[3] += c2;
+        limbs[4] += c3;
 
         FieldElement64(limbs)
     }
@@ -260,18 +325,24 @@ impl FieldElement64 {
     /// Serialize this `FieldElement64` to a 32-byte array.  The
     /// encoding is canonical.
     pub fn to_bytes(&self) -> [u8; 32] {
-        // This reduces to the range [0,2^255), but we need [0,2^255-19).
-        let mut limbs = FieldElement64::reduce(self.0).0;
-
         // Let h = limbs[0] + limbs[1]*2^51 + ... + limbs[4]*2^204.
         //
-        // Write h = pq + r with 0 <= r < p.  We want to compute r = h mod p.
+        // Write h = pq + r with 0 <= r < p.
         //
-        // Since h < 2^255, q = 0 or 1, with q = 0 when h < p and q = 1 when h >= p.
+        // We want to compute r = h mod p.
+        //
+        // If h < 2*p = 2^256 - 38,
+        // then q = 0 or 1,
+        //
+        // with q = 0 when h < p
+        //  and q = 1 when h >= p.
         //
         // Notice that h >= p <==> h + 19 >= p + 19 <==> h + 19 >= 2^255.
         // Therefore q can be computed as the carry bit of h + 19.
 
+        // First, reduce the limbs to ensure h < 2*p.
+        let mut limbs = FieldElement64::reduce(self.0).0;
+
         let mut q = (limbs[0] + 19) >> 51;
         q = (limbs[1] + q) >> 51;
         q = (limbs[2] + q) >> 51;
@@ -337,73 +408,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 + epsilon).
+            a[1] += a[0] >> 51;
+            a[0] &= LOW_51_BIT_MASK;
+
+            // Now all a[i] < 2^(51 + epsilon) 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().
     ///
@@ -178,6 +169,64 @@ impl FieldElement {
         (t19, t3)
     }
 
+    /// Given a slice of public `FieldElements`, replace each with its inverse.
+    ///
+    /// All input `FieldElements` **MUST** be nonzero.
+    ///
+    /// This function is most efficient when the batch size (slice
+    /// length) is a power of 2.
+    #[cfg(any(feature = "alloc", feature = "std"))]
+    pub fn batch_invert(inputs: &mut [FieldElement]) {
+        // First, compute the product of all inputs using a product
+        // tree:
+        //
+        // Inputs: [x_0, x_1, x_2]
+        //
+        // Tree:
+        //
+        //                 x_0*x_1*x_2*1         tree[1]
+        //                   /       \
+        //               x_0*x_1     x_2*1       tree[2,3]
+        //                / \        / \
+        //              x_0  x_1   x_2  1        tree[4,5,6,7]
+        //
+        //  The leaves of the tree are the inputs.  We store the tree in
+        //  an array of length 2*n, similar to a binary heap.
+        //
+        //  To initialize the tree, set every node to 1, then fill in
+        //  the leaf nodes with the input variables.  Finally, set every
+        //  non-leaf node to be the product of its children.
+
+        let n = inputs.len().next_power_of_two();
+        let mut tree = vec![FieldElement::one(); 2*n];
+        tree[n..n+inputs.len()].copy_from_slice(inputs);
+        for i in (1..n).rev() {
+            tree[i] = &tree[2*i] * &tree[2*i+1];
+        }
+
+        // The root of the tree is the product of all inputs, and is
+        // stored at index 1.  Compute its inverse.
+        let allinv = tree[1].invert();
+
+        // To compute y_i = 1/x_i, start at the i-th leaf node of the
+        // tree, and walk up to the root of the tree, multiplying
+        // `allinv` by each sibling.  This computes
+        //
+        // y_i = y * (all x_j except x_i)
+        //
+        // using lg(n) multiplications for each y_i, taking n*lg(n) in
+        // total.
+        for i in 0..inputs.len() {
+            let mut inv = allinv;
+            let mut node = n + i;
+            while node > 1 {
+                inv *= &tree[node ^ 1];
+                node = node >> 1;
+            }
+            inputs[i] = inv;
+        }
+    }
+
     /// Given a nonzero field element, compute its inverse.
     ///
     /// The inverse is computed as self^(p-2), since
@@ -375,6 +424,21 @@ mod test {
         assert_eq!(FieldElement::one(), &a * &should_be_inverse);
     }
 
+    #[test]
+    fn batch_invert_a_matches_nonbatched() {
+        let a    = FieldElement::from_bytes(&A_BYTES);
+        let ap58 = FieldElement::from_bytes(&AP58_BYTES);
+        let asq  = FieldElement::from_bytes(&ASQ_BYTES);
+        let ainv = FieldElement::from_bytes(&AINV_BYTES);
+        let a2   = &a + &a;
+        let a_list = vec![a, ap58, asq, ainv, a2];
+        let mut ainv_list = a_list.clone();
+        FieldElement::batch_invert(&mut ainv_list[..]);
+        for i in 0..5 {
+            assert_eq!(a_list[i].invert(), ainv_list[i]);
+        }
+    }
+
     #[test]
     fn a_p58_vs_ap58_constant() {
         let a    = FieldElement::from_bytes(&A_BYTES);
@@ -470,4 +534,25 @@ mod bench {
         let a = FieldElement::from_bytes(&A_BYTES);
         b.iter(|| a.invert());
     }
+
+    #[bench]
+    fn batch_16_inv(b: &mut Bencher) {
+        let a = FieldElement::from_bytes(&A_BYTES);
+        let mut a_vec = vec![a; 16];
+        b.iter(|| FieldElement::batch_invert(&mut a_vec));
+    }
+
+    #[bench]
+    fn batch_128_inv(b: &mut Bencher) {
+        let a = FieldElement::from_bytes(&A_BYTES);
+        let mut a_vec = vec![a; 128];
+        b.iter(|| FieldElement::batch_invert(&mut a_vec));
+    }
+
+    #[bench]
+    fn batch_1024_inv(b: &mut Bencher) {
+        let a = FieldElement::from_bytes(&A_BYTES);
+        let mut a_vec = vec![a; 1024];
+        b.iter(|| FieldElement::batch_invert(&mut a_vec));
+    }
 }

src/ristretto.rs

@@ -666,6 +666,87 @@ impl RistrettoPoint {
         CompressedRistretto(s.to_bytes())
     }
 
+    /// Double-and-compress a batch of points.
+    #[cfg(any(feature = "alloc", feature = "std"))]
+    pub fn double_and_compress_batch<'a, I>(points: I) -> Vec<CompressedRistretto> 
+        where I: IntoIterator<Item = &'a RistrettoPoint>
+    {
+        #[derive(Copy, Clone, Debug)]
+        struct BatchCompressState {
+            e: FieldElement,
+            f: FieldElement,
+            g: FieldElement,
+            h: FieldElement,
+            eg: FieldElement,
+            fh: FieldElement,
+        }
+
+        impl BatchCompressState {
+            fn efgh(&self) -> FieldElement {
+                &self.eg * &self.fh
+            }
+        }
+
+        impl<'a> From<&'a RistrettoPoint> for BatchCompressState {
+            fn from(P: &'a RistrettoPoint) -> BatchCompressState {
+                let XX = P.0.X.square();
+                let YY = P.0.Y.square();
+                let ZZ = P.0.Z.square();
+                let dTT = &P.0.T.square() * &constants::EDWARDS_D;
+
+                let e = &P.0.X * &(&P.0.Y + &P.0.Y); // = 2*X*Y
+                let f = &ZZ + &dTT;                  // = Z^2 + d*T^2
+                let g = &YY + &XX;                   // = Y^2 - a*X^2
+                let h = &ZZ - &dTT;                  // = Z^2 - d*T^2
+
+                let eg = &e * &g;
+                let fh = &f * &h;
+
+                BatchCompressState{ e: e, f: f, g: g, h: h, eg: eg, fh: fh }
+            }
+        }
+
+        let states: Vec<BatchCompressState> = points.into_iter().map(|P| BatchCompressState::from(P)).collect();
+
+        let mut invs: Vec<FieldElement> = states.iter().map(|state| state.efgh()).collect();
+
+        FieldElement::batch_invert(&mut invs[..]);
+
+        states.iter().zip(invs.iter()).map(|(state, inv): (&BatchCompressState, &FieldElement)| {
+            let Zinv = &state.eg * &inv;
+            let Tinv = &state.fh * &inv;
+
+            let mut magic = constants::INVSQRT_A_MINUS_D;
+
+            let negcheck1 = (&state.eg * &Zinv).is_negative();
+
+            let mut e = state.e;
+            let mut g = state.g;
+            let mut h = state.h;
+
+            let minus_e = -&e;
+            let f_times_sqrta = &state.f * &constants::SQRT_M1;
+
+            e.conditional_assign(&state.g,       negcheck1);
+            g.conditional_assign(&minus_e,       negcheck1);
+            h.conditional_assign(&f_times_sqrta, negcheck1);
+
+            magic.conditional_assign(&constants::SQRT_M1, negcheck1);
+
+            let negcheck2 = (&(&h * &e) * &Zinv).is_negative();
+
+            g.conditional_negate(negcheck2);
+            
+            let mut s = &(&h - &g) * &(&magic * &(&g * &Tinv));
+
+            let s_is_negative = s.is_negative();
+            s.conditional_negate(s_is_negative);
+
+            CompressedRistretto(s.to_bytes())
+        }).collect()
+    }
+
+
     /// Return the coset self + E[4], for debugging.
     fn coset4(&self) -> [ExtendedPoint; 4] {
         [  self.0
@@ -1216,6 +1297,20 @@ mod test {
         }
     }
 
+    #[test]
+    fn double_and_compress_1024_random_points() {
+        let mut rng = OsRng::new().unwrap();
+
+        let points: Vec<RistrettoPoint> = 
+            (0..1024).map(|_| RistrettoPoint::random(&mut rng)).collect();
+
+        let compressed = RistrettoPoint::double_and_compress_batch(&points);
+        
+        for (P, P2_compressed) in points.iter().zip(compressed.iter()) {
+            assert_eq!(*P2_compressed, (P + P).compress());
+        }
+    }
+
     #[test]
     fn random_is_valid() {
         let mut rng = OsRng::new().unwrap();
@@ -1254,4 +1349,28 @@ mod bench {
         let P = B * &Scalar::random(&mut rng);
         b.iter(|| P.compress());
     }
+
+    fn double_and_compress_n_random_points(n: usize, b: &mut Bencher) {
+        let mut rng = OsRng::new().unwrap();
+
+        let points: Vec<RistrettoPoint> = 
+            (0..n).map(|_| RistrettoPoint::random(&mut rng)).collect();
+
+        b.iter(|| RistrettoPoint::double_and_compress_batch(&points) );
+    }
+
+    #[bench]
+    fn double_and_compress_16_random_points(b: &mut Bencher) {
+        double_and_compress_n_random_points(16, b);
+    }
+
+    #[bench]
+    fn double_and_compress_128_random_points(b: &mut Bencher) {
+        double_and_compress_n_random_points(128, b);
+    }
+
+    #[bench]
+    fn double_and_compress_1024_random_points(b: &mut Bencher) {
+        double_and_compress_n_random_points(1024, b);
+    }
 }

vendor/ristretto.sage

@@ -49,7 +49,11 @@ def isqrt(x,exn=InvalidEncodingException("Not on curve")):
     """Return 1/sqrt(x)"""
     if x==0: return 0
     if not is_square(x): raise exn
-    return 1/sqrt(x)
+    s = sqrt(x)
+    #if negative(s): s=-s
+    return 1/s
+
+def inv0(x): return 1/x if x != 0 else 0
 
 def isqrt_i(x):
     """Return 1/sqrt(x) or 1/sqrt(zeta * x)"""
@@ -117,16 +121,20 @@ class QuotientEdwardsPoint(object):
         else:
             return self.__class__(-self.x, -self.y)
     
+    def doubleAndEncodeSpec(self):
+        return (self+self).encode()
 
     # Utility functions
     @classmethod
-    def bytesToGf(cls,bytes,mustBeProper=True,mustBePositive=False):
+    def bytesToGf(cls,bytes,mustBeProper=True,mustBePositive=False,maskHiBits=False):
         """Convert little-endian bytes to field element, sanity check length"""
         if len(bytes) != cls.encLen:
             raise InvalidEncodingException("wrong length %d" % len(bytes))
         s = dec_le(bytes)
-        if mustBeProper and s >= cls.F.modulus():
+        if mustBeProper and s >= cls.F.order():
             raise InvalidEncodingException("%d out of range!" % s)
+        bitlen = int(ceil(log(cls.F.order())/log(2)))
+        if maskHiBits: s &= 2^bitlen-1
         s = cls.F(s)
         if mustBePositive and negative(s):
             raise InvalidEncodingException("%d is negative!" % s)
@@ -197,7 +205,42 @@ class RistrettoPoint(QuotientEdwardsPoint):
             if negative(isr^2*num*y*t): y = -y
             s = isr*y*(z-y)
             
+        return self.gfToBytes(s,mustBePositive=True)
+      
+    @optimized_version_of("doubleAndEncodeSpec")
+    def doubleAndEncode(self):
+        X,Y,Z,T = self.xyzt()
+        a,d,mneg = self.a,self.d,self.mneg
+
+        if self.cofactor==8:
+            e = 2*X*Y
+            f = Z^2+d*T^2
+            g = Y^2-a*X^2
+            h = Z^2-d*T^2
+            
+            inv1 = 1/(e*f*g*h)
+            z_inv = inv1*e*g # 1 / (f*h)
+            t_inv = inv1*f*h
         
+            if negative(e*g*z_inv):
+                if a==-1: sqrta = self.i
+                else:     sqrta = -1
+                e,f,g,h = g,h,-e,f*sqrta
+                factor = self.i
+            else:
+                factor = self.magic
+            
+            if negative(h*e*z_inv): g=-g
+            s = (h-g)*factor*g*t_inv
+            
+        else:
+            foo = Y^2+a*X^2
+            bar = X*Y
+            den = 1/(foo*bar)
+            if negative(2*bar^2*den): tmp = a*X^2
+            else: tmp = Y^2
+            s = self.magic*(Z^2-tmp)*foo*den
+            
         return self.gfToBytes(s,mustBePositive=True)
         
     @classmethod
@@ -238,8 +281,9 @@ class RistrettoPoint(QuotientEdwardsPoint):
     @classmethod
     def elligatorSpec(cls,r0):
         a,d = cls.a,cls.d
-        r = cls.qnr * cls.bytesToGf(r0)^2
+        r = cls.qnr * cls.bytesToGf(r0,mustBeProper=False,maskHiBits=True)^2
         den = (d*r-a)*(a*r-d)
+        if den == 0: return cls()
         n1 = cls.a*(r+1)*(a+d)*(d-a)/den
         n2 = r*n1
         if is_square(n1):
@@ -253,7 +297,7 @@ class RistrettoPoint(QuotientEdwardsPoint):
     @optimized_version_of("elligatorSpec")
     def elligator(cls,r0):
         a,d = cls.a,cls.d
-        r0 = cls.bytesToGf(r0)
+        r0 = cls.bytesToGf(r0,mustBeProper=False,maskHiBits=True)
         r = cls.qnr * r0^2
         den = (d*r-a)*(a*r-d)
         num = cls.a*(r+1)*(a+d)*(d-a)
@@ -278,15 +322,11 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
         
         if self.cofactor==8 and negative(x*y*self.isoMagic):
             x,y = self.torque()
-        
+            
-        isr2 = isqrt(a*(y^2-1)) * sqrt(a*d-1)
-        
         sr = xsqrt(1-a*x^2)
-        assert sr in [isr2*x*y,-isr2*x*y]
+        altx = x*y*self.isoMagic / sr
-        
+        if negative(altx): s = (1+sr)/x
-        altx = 1/isr2*self.isoMagic
+        else:              s = (1-sr)/x
-        if negative(altx): s = (1+x*y*isr2)/(a*x)
-        else:              s = (1-x*y*isr2)/(a*x)
         
         return self.gfToBytes(s,mustBePositive=True)
         
@@ -297,52 +337,141 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
         s = cls.bytesToGf(s,mustBePositive=True)
         
         if s==0: return cls()
-        isr = isqrt(s^4 + 2*(a-2*d)*s^2 + 1)
+        t = xsqrt(s^4 + 2*(a-2*d)*s^2 + 1)
-        altx = 2*s*isr*cls.isoMagic
+        altx = 2*s*cls.isoMagic/t
-        if negative(altx): isr = -isr
+        if negative(altx): t = -t
         x = 2*s / (1+a*s^2)
-        y = (1-a*s^2) * isr
+        y = (1-a*s^2) / t
         
         if cls.cofactor==8 and (negative(x*y*cls.isoMagic) or y==0):
             raise InvalidEncodingException("x*y is invalid: %d, %d" % (x,y))
         
         return cls(x,y)
 
-    @optimized_version_of("encodeSpec")
+    def toJacobiQuartic(self,toggle_rotation=False,toggle_altx=False,toggle_s=False):
-    def encode(self):
+        "Return s,t on jacobi curve"
-        """Encode, optimized version"""
         a,d = self.a,self.d
         x,y,z,t = self.xyzt()
         
         if self.cofactor == 8:
             # Cofactor 8 version
+            # Simulate IMAGINE_TWIST because that's how libdecaf does it
+            x = self.i*x
+            t = self.i*t
+            a = -a
+            d = -d
+            
+            # OK, the actual libdecaf code should be here
             num = (z+y)*(z-y)
             den = x*y
-            tmp = isqrt(num*(a-d)*den^2)
+            isr = isqrt(num*(a-d)*den^2)
-        
+    
-            if negative(tmp^2*den*num*(a-d)*t^2*self.isoMagic):
+            iden = isr * den * self.isoMagic # 1/sqrt((z+y)(z-y)) = 1/sqrt(1-Y^2) / z
-                den,num = num,den
+            inum = isr * num # sqrt(1-Y^2) * z / xysqrt(a-d) ~ 1/sqrt(1-ax^2)/z
-                tmp *= sqrt(a-d) # witness that cofactor is 8
-                yisr = x*sqrt(a)
-                toggle = (a==1)
-            else:
-                yisr = y*(a*d-1)
-                toggle = False
             
-            tiisr = tmp*num
+            if negative(iden*inum*self.i*t^2*(d-a)) != toggle_rotation:
-            altx = tiisr*t*self.isoMagic
+                iden,inum = inum,iden
-            if negative(altx) != toggle: tiisr =- tiisr
+                fac = x*sqrt(a)
-            s = tmp*den*yisr*(tiisr*z - 1)
+                toggle=(a==-1)
+            else:
+                fac = y
+                toggle=False
+            
+            imi = self.isoMagic * self.i
+            altx = inum*t*imi
+            neg_altx = negative(altx) != toggle_altx
+            if neg_altx != toggle: inum =- inum
+            
+            tmp = fac*(inum*z + 1)
+            s = iden*tmp*imi
+            
+            negm1 = (negative(s) != toggle_s) != neg_altx
+            if negm1: m1 = a*fac + z
+            else:     m1 = a*fac - z
+            
+            swap = toggle_s
         
         else:
             # Much simpler cofactor 4 version
             num = (x+t)*(x-t)
             isr = isqrt(num*(a-d)*x^2)
-            ratio = isr*num
+            ratio = isr*num 
-            if negative(ratio*self.isoMagic): ratio=-ratio
+            altx = ratio*self.isoMagic
-            s = (a-d)*isr*x*(ratio*z - t)
+            
+            neg_altx = negative(altx) != toggle_altx
+            if neg_altx: ratio =- ratio
+                
+            tmp = ratio*z - t
+            s = (a-d)*isr*x*tmp
+            
+            negx = (negative(s) != toggle_s) != neg_altx
+            if negx: m1 = -a*t + x
+            else:    m1 = -a*t - x
+            
+            swap = toggle_s
+            
+        if negative(s): s = -s
         
-        return self.gfToBytes(s,mustBePositive=True)
+        return s,m1,a*tmp,swap
+    
+    def invertElligator(self,toggle_r=False,*args,**kwargs):
+        "Produce preimage of self under elligator, or None"
+        a,d = self.a,self.d
+        
+        rets = []
+        
+        tr = [False,True] if self.cofactor == 8 else [False]
+        for toggle_rotation in tr:
+            for toggle_altx in [False,True]:
+                for toggle_s in [False,True]:
+                    for toggle_r in [False,True]:
+                        s,m1,m12,swap = self.toJacobiQuartic(toggle_rotation,toggle_altx,toggle_s)
+
+                        #print
+                        #print toggle_rotation,toggle_altx,toggle_s
+                        #print m1
+                        #print m12
+                    
+                    
+                        if self == self.__class__():
+                            if self.cofactor == 4:
+                                # Hacks for identity!
+                                if toggle_altx: m12 = 1
+                                elif toggle_s: m1 = 1
+                                elif toggle_r: continue
+                                ## BOTH???
+                                
+                            else:
+                                m12 = 1
+                                imi = self.isoMagic * self.i
+                                if toggle_rotation:
+                                    if toggle_altx: m1 = -imi
+                                    else:           m1 = +imi
+                                else:
+                                    if toggle_altx: m1 = 0
+                                    else: m1 = a-d
+                    
+                        rnum = (d*a*m12-m1)
+                        rden = ((d*a-1)*m12+m1)
+                        if swap: rnum,rden = rden,rnum
+                    
+                        ok,sr = isqrt_i(rnum*rden*self.qnr)
+                        if not ok: continue
+                        sr *= rnum
+                        #print "Works! %d %x" % (swap,sr)
+                    
+                        if negative(sr) != toggle_r: sr = -sr
+                        ret = self.gfToBytes(sr)
+                        if self.elligator(ret) != self and self.elligator(ret) != -self:
+                            print "WRONG!",[toggle_rotation,toggle_altx,toggle_s]
+                        if self.elligator(ret) == -self and self != -self: print "Negated!",[toggle_rotation,toggle_altx,toggle_s]
+                        rets.append(bytes(ret))
+        return rets
+
+    @optimized_version_of("encodeSpec")
+    def encode(self):
+        """Encode, optimized version"""    
+        return self.gfToBytes(self.toJacobiQuartic()[0])
         
     @classmethod
     @optimized_version_of("decodeSpec")
@@ -351,7 +480,7 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
         a,d = cls.a,cls.d
         s = cls.bytesToGf(s,mustBePositive=True)
         
-        if s==0: return cls()
+        #if s==0: return cls()
         s2 = s^2
         den = 1+a*s2
         num = den^2 - 4*d*s2
@@ -374,13 +503,63 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
         x = 2*s / (1+a*s^2)
         y = (1-a*s^2) / t
         return cls(x,sgn*y)
+
+    @optimized_version_of("doubleAndEncodeSpec")
+    def doubleAndEncode(self):
+        X,Y,Z,T = self.xyzt()
+        a,d = self.a,self.d
+        
+        if self.cofactor == 8:
+            # Cofactor 8 version
+            # Simulate IMAGINE_TWIST because that's how libdecaf does it
+            X = self.i*X
+            T = self.i*T
+            a = -a
+            d = -d
+            # TODO: This is only being called for a=-1, so could
+            # be wrong for a=1
+            
+            e = 2*X*Y
+            f = Y^2+a*X^2
+            g = Y^2-a*X^2
+            h = Z^2-d*T^2
+            
+            eim = e*self.isoMagic
+            inv = 1/(eim*g*f*h)
+            fh_inv = eim*g*inv*self.i
+            
+            if negative(eim*g*fh_inv):
+                idf = g*self.isoMagic*self.i
+                bar = f
+                foo = g
+                test = eim*f
+            else:
+                idf = eim
+                bar = h
+                foo = -eim
+                test = g*h
+            
+            if negative(test*fh_inv): bar =- bar
+            s = idf*(foo+bar)*inv*f*h
+        
+        else:
+            xy = X*Y
+            h = Z^2-d*T^2
+            inv = 1/(xy*h)
+            if negative(inv*2*xy^2*self.isoMagic): tmp = Y
+            else: tmp = X
+            s = tmp^2*h*inv # = X/Y or Y/X, interestingly
+            
+        return self.gfToBytes(s,mustBePositive=True)
             
     @classmethod
-    def elligatorSpec(cls,r0):
+    def elligatorSpec(cls,r0,fromR=False):
         a,d = cls.a,cls.d
-        r = cls.qnr * cls.bytesToGf(r0)^2
+        if fromR: r = r0
+        else: r = cls.qnr * cls.bytesToGf(r0,mustBeProper=False,maskHiBits=True)^2
         
         den = (d*r-(d-a))*((d-a)*r-d)
+        if den == 0: return cls()
         n1 = (r+1)*(a-2*d)/den
         n2 = r*n1
         if is_square(n1):
@@ -394,7 +573,7 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
     @optimized_version_of("elligatorSpec")
     def elligator(cls,r0):
         a,d = cls.a,cls.d
-        r0 = cls.bytesToGf(r0)
+        r0 = cls.bytesToGf(r0,mustBeProper=False,maskHiBits=True)
         r = cls.qnr * r0^2
         den = (d*r-(d-a))*((d-a)*r-d)
         num = (r+1)*(a-2*d)
@@ -408,6 +587,40 @@ class Decaf_1_1_Point(QuotientEdwardsPoint):
         if negative(s) == iss: s = -s
         return cls.fromJacobiQuartic(s,t)
             
+    def elligatorInverseBruteForce(self):
+        """Invert Elligator using SAGE's polynomial solver"""
+        a,d = self.a,self.d
+        R.<r0> = self.F[]
+        r = self.qnr * r0^2
+        den = (d*r-(d-a))*((d-a)*r-d)
+        n1 = (r+1)*(a-2*d)/den
+        n2 = r*n1
+        ret = set()
+        for s2,t in [(n1, -(r-1)*(a-2*d)^2 / den - 1),
+                     (n2,r*(r-1)*(a-2*d)^2 / den - 1)]:
+            x2 = 4*s2/(1+a*s2)^2
+            y = (1-a*s2) / t
+
+            selfT = self
+            for i in xrange(self.cofactor/2):
+                xT,yT = selfT
+                polyX = xT^2-x2
+                polyY = yT-y
+                sx = set(r for r,_ in polyX.numerator().roots())
+                sy = set(r for r,_ in polyY.numerator().roots())
+                ret = ret.union(sx.intersection(sy))
+            
+                selfT = selfT.torque()
+
+        ret = [self.gfToBytes(r) for r in ret]
+        
+        for r in ret:
+            assert self.elligator(r) in [self,-self]
+            
+        ret = [r for r in ret if self.elligator(r) == self]
+
+        return ret
+            
 class Ed25519Point(RistrettoPoint):
     F = GF(2^255-19)
     d = F(-121665/121666)
@@ -455,7 +668,7 @@ class IsoEd448Point(RistrettoPoint):
     @classmethod
     def base(cls):
         return cls(  # RFC has it wrong
-         -345397493039729516374008604150537410266655260075183290216406970281645695073672344430481787759340633221708391583424041788924124567700732,
+         345397493039729516374008604150537410266655260075183290216406970281645695073672344430481787759340633221708391583424041788924124567700732,
             -363419362147803445274661903944002267176820680343659030140745099590306164083365386343198191849338272965044442230921818680526749009182718
         )
             
@@ -464,7 +677,6 @@ class TwistedEd448GoldilocksPoint(Decaf_1_1_Point):
     d = F(-39082)
     a = F(-1)
     qnr = -1
-    magic = isqrt(a*d-1)
     cofactor = 4
     encLen = 56
     isoMagic = IsoEd448Point.magic
@@ -478,14 +690,13 @@ class Ed448GoldilocksPoint(Decaf_1_1_Point):
     d = F(-39081)
     a = F(1)
     qnr = -1
-    magic = isqrt(a*d-1)
     cofactor = 4
     encLen = 56
     isoMagic = IsoEd448Point.magic
     
     @classmethod
     def base(cls):
-        return -2*cls( # FIXME: make not negative
+        return 2*cls(
  224580040295924300187604334099896036246789641632564134246125461686950415467406032909029192869357953282578032075146446173674602635247710, 298819210078481492676017930443930673437544040154080242095928241372331506189835876003536878655418784733982303233503462500531545062832660
         )
 
@@ -532,19 +743,29 @@ def test(cls,n):
         
     
     P = cls.base()
-    print "base", list(P.encode())
-    for i in xrange(16):
-	Q = P*i
-	print i, list(Q.encode())
-
     Q = cls()
     for i in xrange(n):
-        #print i, binascii.hexlify(Q.encode())
+        #print binascii.hexlify(Q.encode())
-        QQ = cls.decode(Q.encode())
+        QE = Q.encode()
+        QQ = cls.decode(QE)
         if QQ != Q: raise TestFailedException("Round trip %s != %s" % (str(QQ),str(Q)))
+    
+        # Testing s -> 1/s: encodes -point on cofactor 
+        s = cls.bytesToGf(QE)
+        if s != 0:
+            ss = cls.gfToBytes(1/s,mustBePositive=True)
+            try:
+                QN = cls.decode(ss)
+                if cls.cofactor == 8:
+                    raise TestFailedException("1/s shouldnt work for cofactor 8")
+                if QN != -Q:
+                    raise TestFailedException("s -> 1/s should negate point for cofactor 4")
+            except InvalidEncodingException as e:
+                # Should be raised iff cofactor==8
+                if cls.cofactor == 4:
+                    raise TestFailedException("s -> 1/s should work for cofactor 4")
         
         QT = Q
-        QE = Q.encode()
         for h in xrange(cls.cofactor):
             QT = QT.torque()
             if QT.encode() != QE:
@@ -559,27 +780,26 @@ def test(cls,n):
         Q2 = Q0*(r+1)
         if Q1 + Q0 != Q2: raise TestFailedException("Scalarmul doesn't work")
         Q = Q1
-
-test(Ed25519Point,100)
-#test(NegEd25519Point,100)
-#test(IsoEd25519Point,100)
-#test(IsoEd448Point,100)
-#test(TwistedEd448GoldilocksPoint,100)
-#test(Ed448GoldilocksPoint,100)
         
-   
 def testElligator(cls,n):
     print "Testing elligator on %s" % cls.__name__
     for i in xrange(n):
         r = randombytes(cls.encLen)
-        Q = cls.elligator(r)
+        P = cls.elligator(r)
-	print list(r), list(Q.encode())
+        if hasattr(P,"invertElligator"):
-
+            iv = P.invertElligator()
-testElligator(Ed25519Point,100)
+            modr = bytes(cls.gfToBytes(cls.bytesToGf(r,mustBeProper=False,maskHiBits=True)))
-#testElligator(NegEd25519Point,100)
+            iv2 = P.torque().invertElligator()
-#testElligator(IsoEd448Point,100)
+            if modr not in iv: print "Failed to invert Elligator!"
-#testElligator(Ed448GoldilocksPoint,100)
+            if len(iv) != len(set(iv)):
-#testElligator(TwistedEd448GoldilocksPoint,100)
+                print "Elligator inverses not unique!", len(set(iv)), len(iv)
+            if iv != iv2:
+                print "Elligator is untorqueable!"
+                #print [binascii.hexlify(j) for j in iv]
+                #print [binascii.hexlify(j) for j in iv2]
+                #break
+        else:
+            pass # TODO
 
 def gangtest(classes,n):
     print "Gang test",[cls.__name__ for cls in classes]
@@ -607,5 +827,31 @@ def gangtest(classes,n):
             for c,ret in zip(classes,rets):
                 print c,binascii.hexlify(ret)
             print
-gangtest([IsoEd448Point,TwistedEd448GoldilocksPoint,Ed448GoldilocksPoint],100)
+
-gangtest([Ed25519Point,IsoEd25519Point],100)
+def testDoubleAndEncode(cls,n):
+    print "Testing doubleAndEncode on %s" % cls.__name__
+    for i in xrange(n):
+        r1 = randombytes(cls.encLen)
+        r2 = randombytes(cls.encLen)
+        u = cls.elligator(r1) + cls.elligator(r2)
+        u.doubleAndEncode()
+
+testDoubleAndEncode(Ed25519Point,100)
+testDoubleAndEncode(NegEd25519Point,100)
+testDoubleAndEncode(IsoEd25519Point,100)
+testDoubleAndEncode(IsoEd448Point,100)
+testDoubleAndEncode(TwistedEd448GoldilocksPoint,100)
+#test(Ed25519Point,100)
+#test(NegEd25519Point,100)
+#test(IsoEd25519Point,100)
+#test(IsoEd448Point,100)
+#test(TwistedEd448GoldilocksPoint,100)
+#test(Ed448GoldilocksPoint,100)
+#testElligator(Ed25519Point,100)
+#testElligator(NegEd25519Point,100)
+#testElligator(IsoEd25519Point,100)
+#testElligator(IsoEd448Point,100)
+#testElligator(Ed448GoldilocksPoint,100)
+#testElligator(TwistedEd448GoldilocksPoint,100)
+#gangtest([IsoEd448Point,TwistedEd448GoldilocksPoint,Ed448GoldilocksPoint],100)
+#gangtest([Ed25519Point,IsoEd25519Point],100)