Merge remote-tracking branch 'hdevalence/feature/decaf_r15' into develop

Cargo.toml

@@ -25,6 +25,7 @@ version = "0.3"
 [features]
 default = ["std"]
 std = ["rand"]
+yolocrypto = []
 
 # The development profile, used for `cargo build`.
 [profile.dev]

src/constants.rs

@@ -17,8 +17,10 @@
 #![allow(non_snake_case)]
 #![allow(non_upper_case_globals)]
 #![allow(missing_docs)]
+#![allow(non_snake_case)]
 
 use field::FieldElement;
+use curve::ExtendedPoint;
 use curve::PreComputedPoint;
 use curve::CompressedEdwardsY;
 use scalar::Scalar;
@@ -29,10 +31,35 @@ pub const d: FieldElement       = FieldElement([
 pub const d2: FieldElement      = FieldElement([
     -21827239,  -5839606, -30745221,  13898782,    229458,
     15978800,  -12551817,  -6495438,  29715968,   9444199, ]);
+pub const d4: FieldElement      = FieldElement([
+    23454405,  -11679213,   5618422,  -5756869,    458917,
+    -1596832,  -25103633, -12990876,  -7676928, -14666033  ]);
+pub const a_minus_d: FieldElement    = FieldElement([
+    10913609,  -13857413,  15372611,  -6949391,   -114729,
+     8787816,    6275908,   3247719,  18696448,  12055116, ]);
+
+/// (p-1)/2, in little-endian bytes.
+pub const HALF_P_MINUS_1_BYTES: [u8; 32] =
+    [0xf6, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+     0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+     0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
+     0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x3f];
+
+/// Precomputed value of one of the square roots of -1 (mod p)
 pub const SQRT_M1: FieldElement = FieldElement([
     -32595792,  -7943725,   9377950,   3500415,  12389472,
     -272473,   -25146209,  -2005654,    326686,  11406482, ]);
 
+/// Precomputed value of the other square root of -1 (mod p),
+/// i.e., MSQRT_M1 = -SQRT_M1.
+pub const MSQRT_M1: FieldElement = FieldElement([
+    32595792,    7943725,  -9377950,  -3500415, -12389472,
+    272473,     25146209,   2005654,   -326686, -11406482, ]);
+
+/// Precomputed value of 1/2 (mod p).
+pub const HALF: FieldElement = FieldElement([
+    10, 0, 0, 0, 0, 0, 0, 0, 0, -16777216, ]);
+
 /// In Montgomery form y² = x³+Ax²+x, Curve25519 has A=486662.
 pub const A: FieldElement       = FieldElement([
     486662, 0, 0, 0, 0, 0, 0, 0, 0, 0, ]);
@@ -67,6 +94,14 @@ pub const BASE_CMPRSSD: CompressedEdwardsY =
                         0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66,
                         0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66]);
 
+/// Basepoint has y = 4/5.
+pub const BASEPOINT: ExtendedPoint = ExtendedPoint{
+        X: FieldElement([-14297830, -7645148, 16144683, -16471763, 27570974, -2696100, -26142465, 8378389, 20764389, 8758491]),
+        Y: FieldElement([-26843541, -6710886, 13421773, -13421773, 26843546, 6710886, -13421773, 13421773, -26843546, -6710886]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([28827062, -6116119, -27349572, 244363, 8635006, 11264893, 19351346, 13413597, 16611511, -6414980]),
+};
+
 /// `l` is the order of base point, i.e. 2^252 +
 /// 27742317777372353535851937790883648493, in little-endian form
 pub const l: Scalar = Scalar([ 0xed, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58,
@@ -80,6 +115,63 @@ pub const lminus1: Scalar = Scalar([ 0xec, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0
                                      0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14,
                                      0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
                                      0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10 ]);
+/// The 8-torsion subgroup Ɛ[8].
+///
+/// In the case of Curve25519, it is cyclic; the `i`th element of the
+/// array is `i*P`, where `P` is a point of order 8 generating Ɛ[8].
+///
+/// Thus Ɛ[4] is the points indexed by 0,2,4,6 and Ɛ[2] is the points
+/// indexed by 0,4. 
+pub const EIGHT_TORSION: [ExtendedPoint; 8] = [
+    ExtendedPoint{
+        X: FieldElement([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        Y: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+    },
+    ExtendedPoint{
+        X: FieldElement([21352778, 5345713, 4660180, -8347857, 24143090, 14568123, 30185756, -12247770, -33528939, 8345319]),
+        Y: FieldElement([6952922, 1265500, -6862341, 7057498, 4037696, 5447722, -31680899, 15325402, 19365852, -1569102]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([-25262188, -11972680, 11716002, -5869612, -18193162, 16297739, 20670665, -8559098, 3541543, -5011181])
+    },
+    ExtendedPoint{
+        X: FieldElement([32595792, 7943725, -9377950, -3500415, -12389472, 272473, 25146209, 2005654, -326686, -11406482]),
+        Y: FieldElement([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+    },
+    ExtendedPoint{
+        X: FieldElement([21352778, 5345713, 4660180, -8347857, 24143090, 14568123, 30185756, -12247770, -33528939, 8345319]),
+        Y: FieldElement([-6952922, -1265500, 6862341, -7057498, -4037696, -5447722, 31680899, -15325402, -19365852, 1569102]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([25262188, 11972680, -11716002, 5869612, 18193162, -16297739, -20670665, 8559098, -3541543, 5011181])
+    },
+    ExtendedPoint{
+        X: FieldElement([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        Y: FieldElement([-1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+    },
+    ExtendedPoint{
+        X: FieldElement([-21352778, -5345713, -4660180, 8347857, -24143090, -14568123, -30185756, 12247770, 33528939, -8345319]),
+        Y: FieldElement([-6952922, -1265500, 6862341, -7057498, -4037696, -5447722, 31680899, -15325402, -19365852, 1569102]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([-25262188, -11972680, 11716002, -5869612, -18193162, 16297739, 20670665, -8559098, 3541543, -5011181])
+    },
+    ExtendedPoint{
+        X: FieldElement([-32595792, -7943725, 9377950, 3500415, 12389472, -272473, -25146209, -2005654, 326686, 11406482]),
+        Y: FieldElement([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+    },
+    ExtendedPoint{
+        X: FieldElement([-21352778, -5345713, -4660180, 8347857, -24143090, -14568123, -30185756, 12247770, 33528939, -8345319]),
+        Y: FieldElement([6952922, 1265500, -6862341, 7057498, 4037696, 5447722, -31680899, 15325402, 19365852, -1569102]),
+        Z: FieldElement([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
+        T: FieldElement([25262188, 11972680, -11716002, 5869612, 18193162, -16297739, -20670665, 8559098, -3541543, 5011181])
+    },
+];
 
 pub const bi: [PreComputedPoint; 8] = [
     PreComputedPoint{
@@ -1478,14 +1570,77 @@ pub const base: [[PreComputedPoint; 8]; 32] = [
 mod test {
     use field::FieldElement;
     use curve::PreComputedPoint;
+    use curve::CompressedEdwardsY;
+    use curve::ExtendedPoint;
+    use curve::Identity;
+    use curve::IsIdentity;
+    use curve::ValidityCheck;
     use constants;
 
     #[test]
-    /// Test that SQRT_M1 is a square root of -1
+    fn test_eight_torsion() {
+        let mut bytes = [0;32];
+        bytes[0] = 1;
+        let compressed_id = CompressedEdwardsY(bytes);
+        for i in 0..8 {
+            let Q = constants::EIGHT_TORSION[i].mult_by_pow_2(3);
+            assert!(Q.is_valid());
+            assert!(Q.is_identity());
+        }
+    }
+
+    #[test]
+    fn test_four_torsion() {
+        let mut bytes = [0;32];
+        bytes[0] = 1;
+        let compressed_id = CompressedEdwardsY(bytes);
+        for i in (0..8).filter(|i| i % 2 == 0) {
+            let Q = constants::EIGHT_TORSION[i].mult_by_pow_2(2);
+            assert!(Q.is_valid());
+            assert!(Q.is_identity());
+        }
+    }
+
+    #[test]
+    fn test_two_torsion() {
+        let mut bytes = [0;32];
+        bytes[0] = 1;
+        let compressed_id = CompressedEdwardsY(bytes);
+        for i in (0..8).filter(|i| i % 4 == 0) {
+            let Q = constants::EIGHT_TORSION[i].mult_by_pow_2(1);
+            assert!(Q.is_valid());
+            assert!(Q.is_identity());
+        }
+    }
+
+    #[test]
+    fn test_half() {
+        let one = FieldElement([1,0,0,0,0,0,0,0,0,0]);
+        let two = FieldElement([2,0,0,0,0,0,0,0,0,0]);
+        assert_eq!(one, &two * &constants::HALF);
+    }
+
+    #[test]
+    /// Test that SQRT_M1 and MSQRT_M1 are square roots of -1
     fn test_sqrt_minus_one() {
         let minus_one = FieldElement([-1,0,0,0,0,0,0,0,0,0]);
         let sqrt_m1_sq = &constants::SQRT_M1 * &constants::SQRT_M1;
-        assert_eq!(minus_one, sqrt_m1_sq);
+        let msqrt_m1_sq = &constants::MSQRT_M1 * &constants::MSQRT_M1;
+        assert_eq!(minus_one,  sqrt_m1_sq);
+        assert_eq!(minus_one, msqrt_m1_sq);
+    }
+
+    #[test]
+    fn test_sqrt_constants_sign() {
+        let one       = FieldElement([ 1,0,0,0,0,0,0,0,0,0]);
+        let minus_one = FieldElement([-1,0,0,0,0,0,0,0,0,0]);
+        let invsqrt_m1 = minus_one.invsqrt().unwrap();
+        let sign_test_sqrt  = &invsqrt_m1 * &constants::SQRT_M1;
+        let sign_test_msqrt = &invsqrt_m1 * &constants::MSQRT_M1;
+        // XXX it seems we have flipped the sign relative to
+        // the invsqrt function?
+        assert_eq!(sign_test_sqrt, minus_one);
+        assert_eq!(sign_test_msqrt, one);
     }
 
     #[test]
@@ -1499,6 +1654,20 @@ mod test {
         assert_eq!(d2, constants::d2);
     }
 
+    #[test]
+    fn test_d4() {
+        let mut four = FieldElement::zero();
+        four[0] = 4;
+        assert_eq!(&constants::d * &four, constants::d4);
+    }
+
+    #[test]
+    fn test_a_minus_d() {
+        let a = FieldElement([-1,0,0,0,0,0,0,0,0,0]);
+        let a_minus_d = &a - &constants::d;
+        assert_eq!(a_minus_d, constants::a_minus_d);
+    }
+
     /// Test the values in the lookup table of precomputed multiples
     /// of the basepoint.
     #[test]

src/curve.rs

@@ -80,6 +80,7 @@
 use core::fmt::Debug;
 use core::iter::Iterator;
 use core::ops::{Add, Sub, Neg, Index};
+use core::cmp::{PartialEq, Eq};
 
 use constants;
 use field::FieldElement;
@@ -159,7 +160,7 @@ impl CompressedEdwardsY {
             X *= &constants::SQRT_M1;
         }
 
-        if X.is_negative() != (self[31] >> 7) as i32 {
+        if X.is_negative_ed25519() != (self[31] >> 7) as i32 {
             X = X.neg();
         }
         T = &X * &Y;
@@ -174,12 +175,16 @@ impl CompressedEdwardsY {
 
 /// An `ExtendedPoint` is a point on the curve in 𝗣³(𝔽ₚ).
 /// A point (x,y) in the affine model corresponds to (x:y:1:xy).
+// XXX members should not be public, but that's needed for the
+// constants module. Fix when RFC #1422 lands:
+// https://github.com/rust-lang/rust/issues/32409
 #[derive(Copy, Clone)]
+#[allow(missing_docs)]
 pub struct ExtendedPoint {
-    X: FieldElement,
+    pub X: FieldElement,
-    Y: FieldElement,
+    pub Y: FieldElement,
-    Z: FieldElement,
+    pub Z: FieldElement,
-    T: FieldElement,
+    pub T: FieldElement,
 }
 
 /// A `ProjectivePoint` is a point on the curve in 𝗣²(𝔽ₚ).
@@ -271,6 +276,38 @@ impl Identity for PreComputedPoint {
     }
 }
 
+// ------------------------------------------------------------------------
+// Validity checks (for debugging, not CT)
+// ------------------------------------------------------------------------
+
+/// Trait for checking whether a point is on the curve
+pub trait ValidityCheck {
+    /// Checks whether the point is on the curve. Not CT.
+    fn is_valid(&self) -> bool;
+}
+
+impl ValidityCheck for ProjectivePoint {
+    fn is_valid(&self) -> bool {
+        // Curve equation is    -x^2 + y^2 = 1 + d*x^2*y^2,
+        // homogenized as (-X^2 + Y^2)*Z^2 = Z^4 + d*X^2*Y^2
+        let XX = self.X.square();
+        let YY = self.Y.square();
+        let ZZ = self.Z.square();
+        let ZZZZ = ZZ.square();
+        let lhs = &(&YY - &XX) * &ZZ;
+        let rhs = &ZZZZ + &(&constants::d * &(&XX * &YY));
+
+        lhs == rhs
+    }
+}
+
+impl ValidityCheck for ExtendedPoint {
+    // XXX this should also check that T is correct
+    fn is_valid(&self) -> bool {
+        self.to_projective().is_valid()
+    }
+}
+
 // ------------------------------------------------------------------------
 // Constant-time assignment
 // ------------------------------------------------------------------------
@@ -355,7 +392,7 @@ impl ProjectivePoint {
         let mut s: [u8; 32];
 
         s      =  y.to_bytes();
-        s[31] ^= (x.is_negative() << 7) as u8;
+        s[31] ^= (x.is_negative_ed25519() << 7) as u8;
         CompressedEdwardsY(s)
     }
 }
@@ -431,7 +468,7 @@ impl CompletedPoint {
 
 impl ProjectivePoint {
     /// Double this point: return self + self
-    fn double(&self) -> CompletedPoint { // Double()
+    pub fn double(&self) -> CompletedPoint { // Double()
         let XX          = self.X.square();
         let YY          = self.Y.square();
         let ZZ2         = self.Z.square2();
@@ -451,7 +488,7 @@ impl ProjectivePoint {
 
 impl ExtendedPoint {
     /// Add this point to itself.
-    fn double(&self) -> ExtendedPoint {
+    pub fn double(&self) -> ExtendedPoint {
         self.to_projective().double().to_extended()
     }
 }
@@ -599,18 +636,24 @@ impl<'a> Neg for &'a PreComputedPoint {
 // Scalar multiplication
 // ------------------------------------------------------------------------
 
-impl ExtendedPoint {
+/// Trait for scalar multiplication of an arbitrary point.
-    /// Scalar multiplication: compute `a * self`.
+pub trait ScalarMult<S> {
+    /// Compute `scalar * self`.
+    fn scalar_mult(&self, scalar: &S) -> Self;
+}
+
+impl ScalarMult<Scalar> for ExtendedPoint {
+    /// Scalar multiplication: compute `scalar * self`.
     ///
     /// Uses a window of size 4.  Note: for scalar multiplication of
     /// the basepoint, `basepoint_mult` is approximately 4x faster.
-    pub fn scalar_mult(&self, a: &Scalar) -> ExtendedPoint {
+    fn scalar_mult(&self, scalar: &Scalar) -> ExtendedPoint {
         let A = self.to_cached();
         let mut As: [CachedPoint; 8] = [A; 8];
         for i in 0..7 {
             As[i+1] = (self + &As[i]).to_extended().to_cached();
         }
-        let e = a.to_radix_16();
+        let e = scalar.to_radix_16();
         let mut h = ExtendedPoint::identity();
         let mut t: CompletedPoint;
         for i in (0..64).rev() {
@@ -620,8 +663,22 @@ impl ExtendedPoint {
         }
         h
     }
+}
 
-    /// Construct an `ExtendedPoint` from a `Scalar`, `a`, by
+/// Trait for scalar multiplication of a distinguished basepoint.
+pub trait BasepointMult<S> {
+    /// Return the basepoint `B`.
+    fn basepoint() -> Self;
+    /// Compute `scalar * B`.
+    fn basepoint_mult(scalar: &S) -> Self;
+}
+
+impl BasepointMult<Scalar> for ExtendedPoint {
+    fn basepoint() -> ExtendedPoint {
+        constants::BASEPOINT
+    }
+
+    /// Construct an `ExtendedPoint` from a `Scalar`, `scalar`, by
     /// computing the multiple `aB` of the basepoint `B`.
     ///
     /// Precondition: the scalar must be reduced.
@@ -646,8 +703,8 @@ impl ExtendedPoint {
     /// We then use the `select_precomputed_point` function, which
     /// takes `-8 ≤ x < 8` and `[16^2i * B, ..., 8 * 16^2i * B]`,
     /// and returns `x * 16^2i * B` in constant time.
-    pub fn basepoint_mult(a: &Scalar) -> ExtendedPoint { //GeScalarMultBase
+    fn basepoint_mult(scalar: &Scalar) -> ExtendedPoint { //GeScalarMultBase
-        let e = a.to_radix_16();
+        let e = scalar.to_radix_16();
         let mut h = ExtendedPoint::identity();
         let mut t: CompletedPoint;
 
@@ -665,7 +722,9 @@ impl ExtendedPoint {
 
         h
     }
+}
 
+impl ExtendedPoint {
     /// Multiply by the cofactor: compute `8 * self`.
     ///
     /// Convenience wrapper around `mult_by_pow_2`.
@@ -862,6 +921,8 @@ impl Debug for CachedPoint {
 #[cfg(test)]
 mod test {
     use test::Bencher;
+    use rand::OsRng;
+
     use field::FieldElement;
     use scalar::Scalar;
     use subtle::CTAssignable;
@@ -922,6 +983,8 @@ mod test {
         let base_X = FieldElement::from_bytes(&BASE_X_COORD_BYTES);
         let bp  =  BASE_CMPRSSD.decompress().unwrap();
         let bp2 = BASE2_CMPRSSD.decompress().unwrap();
+        assert!( bp.is_valid());
+        assert!(bp2.is_valid());
         let compressed  =  bp.compress();
         let compressed2 = bp2.compress();
         // Check that decompression actually gives the correct X coordinate

src/decaf.rs

@@ -0,0 +1,389 @@
+// -*- mode: rust; -*-
+//
+// To the extent possible under law, the authors have waived all copyright and
+// related or neighboring rights to curve25519-dalek, using the Creative
+// Commons "CC0" public domain dedication.  See
+// <http://creativecommons.org/publicdomain/zero/.0/> for full details.
+//
+// Authors:
+// - Isis Agora Lovecruft <isis@patternsinthevoid.net>
+// - Henry de Valence <hdevalence@hdevalence.ca>
+
+//! An implementation of Mike Hamburg's Decaf cofactor-eliminating
+//! point-compression scheme, providing a prime-order group on top of
+//! a non-prime-order elliptic curve.
+//!
+//! Note: this code is currently feature-gated with the `yolocrypto`
+//! feature flag, because our implementation is still unfinished.
+
+// We allow non snake_case names because coordinates in projective space are
+// traditionally denoted by the capitalisation of their respective
+// counterparts in affine space.  Yeah, you heard me, rustc, I'm gonna have my
+// affine and projective cakes and eat both of them too.
+#![allow(non_snake_case)]
+
+use core::fmt::Debug;
+
+use constants;
+use field::FieldElement;
+use subtle::CTAssignable;
+use subtle::CTNegatable;
+
+use core::ops::{Add, Sub, Neg};
+
+use curve::ExtendedPoint;
+use curve::BasepointMult;
+use curve::ScalarMult;
+use curve::Identity;
+use scalar::Scalar;
+
+// ------------------------------------------------------------------------
+// Compressed points
+// ------------------------------------------------------------------------
+
+/// A point serialized using Mike Hamburg's Decaf scheme.
+///
+/// XXX think about how this API should work
+#[derive(Copy, Clone, Eq, PartialEq)]
+pub struct CompressedDecaf(pub [u8; 32]);
+
+/// The result of compressing a `DecafPoint`.
+impl CompressedDecaf {
+    /// View this `CompressedDecaf` as an array of bytes.
+    pub fn to_bytes(&self) -> [u8;32] {
+        self.0
+    }
+
+    /// Attempt to decompress to an `DecafPoint`.
+    pub fn decompress(&self) -> Option<DecafPoint> {
+        // XXX should decoding be CT ?
+        // XXX need to check that xy is nonnegative and reject otherwise
+        let s = FieldElement::from_bytes(&self.0);
+
+        // Check that s = |s| and reject otherwise.
+        let mut abs_s = s;
+        let neg = abs_s.is_negative_decaf();
+        abs_s.conditional_negate(neg);
+        if abs_s != s { return None; }
+
+        let ss = s.square();
+        let X = &s + &s;                    // X = 2s
+        let Z = &FieldElement::one() - &ss; // Z = 1+as^2
+        let u = &(&Z * &Z) - &(&constants::d4 * &ss); // u = Z^2 - 4ds^2
+        let uss = &u * &ss;
+        let mut v = match uss.invsqrt() {
+            Some(v) => v,
+            None => return None,
+        };
+        // Now v = 1/sqrt(us^2) if us^2 is a nonzero square, 0 if us^2 is zero.
+        let uv = &v * &u;
+        if uv.is_negative_decaf() == 1u8 {
+            v.negate();
+        }
+        let mut two_minus_Z = -&Z; two_minus_Z[0] += 2;
+        let mut w = &v * &(&s * &two_minus_Z);
+        w.conditional_assign(&FieldElement::one(), s.is_zero());
+        let Y = &w * &Z;
+        let T = &w * &X;
+
+        // "To decode the point, one must decode it to affine form
+        // instead of projective, and check that xy is non-negative."
+        //
+        // XXX can we merge this inversion with the one above?
+        let xy = &T * &Z.invert();
+        if (Y.is_nonzero() & xy.is_nonnegative_decaf()) == 1u8 {
+            Some(DecafPoint(ExtendedPoint{ X: X, Y: Y, Z: Z, T: T }))
+        } else {
+            None
+        }
+    }
+}
+
+impl Identity for CompressedDecaf {
+    fn identity() -> CompressedDecaf {
+        CompressedDecaf([0u8;32])
+    }
+}
+
+/// A point in a prime-order group.
+///
+/// XXX think about how this API should work
+#[derive(Copy, Clone)]
+pub struct DecafPoint(pub ExtendedPoint);
+
+impl DecafPoint {
+    /// Compress in Decaf format.
+    pub fn compress(&self) -> CompressedDecaf {
+        // Q: Do we want to encode twisted or untwisted?
+        //
+        // Notes: 
+        // Recall that the twisted Edwards curve E_{a,d} is of the form
+        //
+        //     ax^2 + y^2 = 1 + dx^2y^2. 
+        //
+        // Internally, we operate on the curve with a = -1, d =
+        // -121665/121666, a.k.a., the twist.  But maybe we would like
+        // to use Decaf on the untwisted curve with a = 1, d =
+        // 121665/121666.  (why? interop?)
+        //
+        // Fix i, a square root of -1 (mod p).
+        //
+        // The map x -> ix is an isomorphism from E_{a,d} to E_{-a,-d}. 
+        // Its inverse is x -> -ix.
+        // let untwisted_X = &self.X * &constants::MSQRT_M1;
+        // etc.
+
+        // Step 0: pre-rotation, needed for Decaf with E[8] = Z/8
+
+        let mut X = self.0.X;
+        let mut Y = self.0.Y;
+        let mut T = self.0.T;
+
+        // If y nonzero and xy nonnegative, continue.
+        // Otherwise, add Q_6 = (i,0) = constants::EIGHT_TORSION[6]
+        // (x,y) + Q_6 = (iy,ix)
+        // (X:Y:Z:T) + Q_6 = (iY:iX:Z:-T)
+
+        // XXX it should be possible to avoid this inversion, but
+        // let's make sure the code is correct first
+        let xy = &T * &self.0.Z.invert();
+        let is_neg_mask = 1u8 & !(Y.is_nonzero() & xy.is_nonnegative_decaf());
+        let iX = &X * &constants::SQRT_M1;
+        let iY = &Y * &constants::SQRT_M1;
+        X.conditional_assign(&iY, is_neg_mask);
+        Y.conditional_assign(&iX, is_neg_mask);
+        T.conditional_negate(is_neg_mask);
+
+        // Step 1: Compute r = 1/sqrt((a-d)(Z+Y)(Z-Y))
+        let Z_plus_Y  = &self.0.Z + &Y;
+        let Z_minus_Y = &self.0.Z - &Y;
+        let t = &constants::a_minus_d * &(&Z_plus_Y * &Z_minus_Y);
+        // t should always be square (why?)
+        // XXX is it safe to use option types here?
+        let mut r = t.invsqrt().unwrap();
+
+        // Step 2: Compute u = (a-d)r
+        let u = &constants::a_minus_d * &r;
+
+        // Step 3: Negate r if -2uZ is negative.
+        let uZ = &u * &self.0.Z;
+        let m2uZ = -&(&uZ + &uZ);
+        r.conditional_negate(m2uZ.is_negative_decaf());
+
+        // Step 4: Compute s = | u(r(aZX - dYT)+Y)/a|
+        //                   = |u(r(-ZX - dYT)+Y)| since a = -1
+        let minus_ZX = -&(&self.0.Z * &X);
+        let dYT = &constants::d * &(&Y * &T);
+        // Compute s = u(r(aZX - dYT)+Y) and cnegate for abs
+        let mut s = &u * &(&(&r * &(&minus_ZX - &dYT)) + &Y);
+        let neg = s.is_negative_decaf();
+        s.conditional_negate(neg);
+        CompressedDecaf(s.to_bytes())
+    }
+
+    /// Return the coset self + E[4], for debugging.
+    fn coset4(&self) -> [ExtendedPoint; 4] {
+        [  self.0
+        , &self.0 + &constants::EIGHT_TORSION[2]
+        , &self.0 + &constants::EIGHT_TORSION[4]
+        , &self.0 + &constants::EIGHT_TORSION[6]
+        ]
+    }
+}
+
+impl Identity for DecafPoint {
+    fn identity() -> DecafPoint {
+        DecafPoint(ExtendedPoint::identity())
+    }
+}
+
+// ------------------------------------------------------------------------
+// Equality
+// ------------------------------------------------------------------------
+
+/// XXX check whether there's a simple way to do equality checking
+/// with cofactor 8, not just cofactor 4, and add a CT equality function?
+impl PartialEq for DecafPoint {
+    fn eq(&self, other: &DecafPoint) -> bool {
+        let  self_compressed =  self.compress();
+        let other_compressed = other.compress();
+        self_compressed == other_compressed
+    }
+}
+
+impl Eq for DecafPoint {}
+
+// ------------------------------------------------------------------------
+// Arithmetic
+// ------------------------------------------------------------------------
+
+impl<'a, 'b> Add<&'b DecafPoint> for &'a DecafPoint {
+    type Output = DecafPoint;
+
+    fn add(self, other: &'b DecafPoint) -> DecafPoint {
+        DecafPoint(&self.0 + &other.0)
+    }
+}
+
+impl<'a, 'b> Sub<&'b DecafPoint> for &'a DecafPoint {
+    type Output = DecafPoint;
+
+    fn sub(self, other: &'b DecafPoint) -> DecafPoint {
+        DecafPoint(&self.0 - &other.0)
+    }
+}
+
+impl<'a> Neg for &'a DecafPoint {
+    type Output = DecafPoint;
+
+    fn neg(self) -> DecafPoint {
+        DecafPoint(-&self.0)
+    }
+}
+
+impl ScalarMult<Scalar> for DecafPoint {
+    fn scalar_mult(&self, scalar: &Scalar) -> DecafPoint {
+        DecafPoint(self.0.scalar_mult(scalar))
+    }
+}
+
+impl BasepointMult<Scalar> for DecafPoint {
+    // XXX is this actually in the image of the isogeny,
+    // or do we need a different basepoint?
+    fn basepoint() -> DecafPoint {
+        DecafPoint(constants::BASEPOINT)
+    }
+
+    fn basepoint_mult(scalar: &Scalar) -> DecafPoint {
+        DecafPoint(ExtendedPoint::basepoint_mult(scalar))
+    }
+}
+
+// ------------------------------------------------------------------------
+// Debug traits
+// ------------------------------------------------------------------------
+
+impl Debug for CompressedDecaf {
+    fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result {
+        write!(f, "CompressedDecaf: {:?}", &self.0[..])
+    }
+}
+
+impl Debug for DecafPoint {
+    fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result {
+        let coset = self.coset4();
+        write!(f, "DecafPoint: coset \n{:?}\n{:?}\n{:?}\n{:?}",
+               coset[0], coset[1], coset[2], coset[3])
+    }
+}
+
+// ------------------------------------------------------------------------
+// Tests
+// ------------------------------------------------------------------------
+
+#[cfg(test)]
+mod test {
+    use rand::OsRng;
+
+    use scalar::Scalar;
+    use constants;
+    use constants::BASE_CMPRSSD;
+    use curve::CompressedEdwardsY;
+    use curve::ExtendedPoint;
+    use curve::BasepointMult;
+    use curve::Identity;
+    use super::*;
+
+    #[test]
+    #[should_panic]
+    fn test_decaf_decompress_negative_s_fails() {
+        // constants::d is neg, so decompression should fail as |d| != d.
+        let bad_compressed = CompressedDecaf(constants::d.to_bytes());
+        bad_compressed.decompress().unwrap();
+    }
+
+    #[test]
+    fn test_decaf_decompress_id() {
+        let compressed_id = CompressedDecaf::identity();
+        let id = compressed_id.decompress().unwrap();
+        // This should compress (as ed25519) to the following:
+        let mut bytes = [0u8; 32]; bytes[0] = 1;
+        assert_eq!(id.0.compress(), CompressedEdwardsY(bytes));
+    }
+
+    #[test]
+    fn test_decaf_compress_id() {
+        let id = DecafPoint::identity();
+        assert_eq!(id.compress(), CompressedDecaf::identity());
+    }
+
+    #[test]
+    fn test_decaf_basepoint_roundtrip() {
+        let bp_compressed_decaf = DecafPoint::basepoint().compress();
+        let bp_recaf = bp_compressed_decaf.decompress().unwrap().0;
+        // Check that bp_recaf differs from bp by a point of order 4
+        let diff = &ExtendedPoint::basepoint() - &bp_recaf;
+        let diff4 = diff.mult_by_pow_2(4);
+        assert_eq!(diff4.compress(), ExtendedPoint::identity().compress());
+    }
+
+    #[test]
+    fn test_decaf_four_torsion_basepoint() {
+        let bp = DecafPoint::basepoint();
+        let bp_coset = bp.coset4();
+        for i in 0..4 {
+            assert_eq!(bp, DecafPoint(bp_coset[i]));
+        }
+    }
+
+    #[test]
+    fn test_decaf_four_torsion_random() {
+        let mut rng = OsRng::new().unwrap();
+        let s = Scalar::random(&mut rng);
+        let P = DecafPoint::basepoint_mult(&s);
+        let P_coset = P.coset4();
+        for i in 0..4 {
+            assert_eq!(P, DecafPoint(P_coset[i]));
+        }
+    }
+
+    #[test]
+    fn test_decaf_random_roundtrip() {
+        let mut rng = OsRng::new().unwrap();
+        for j in 0..100 {
+        let s = Scalar::random(&mut rng);
+        let P = DecafPoint::basepoint_mult(&s);
+        let compressed_P = P.compress();
+        let Q = compressed_P.decompress().unwrap();
+        for i in 0..4 {
+            assert_eq!(P, Q);
+        }
+        }
+    }
+}
+
+#[cfg(test)]
+mod bench {
+    use rand::OsRng;
+    use test::Bencher;
+
+    use super::*;
+
+    #[bench]
+    fn decompression(b: &mut Bencher) {
+        let mut rng = OsRng::new().unwrap();
+        let s = Scalar::random(&mut rng);
+        let P = DecafPoint::basepoint_mult(&s);
+        let P_compressed = P.compress();
+        b.iter(|| P_compressed.decompress().unwrap());
+    }
+
+    #[bench]
+    fn compression(b: &mut Bencher) {
+        let mut rng = OsRng::new().unwrap();
+        let s = Scalar::random(&mut rng);
+        let P = DecafPoint::basepoint_mult(&s);
+        b.iter(|| P.compress());
+    }
+}
+

src/field.rs

@@ -31,6 +31,8 @@ use subtle::CTEq;
 
 use utils::{load3, load4};
 
+use constants;
+
 /// FieldElements are represented as an array of ten "Limbs", which are radix
 /// 25.5, that is, each Limb of a FieldElement alternates between being
 /// represented as a factor of 2^25 or 2^26 more than the last corresponding
@@ -510,17 +512,60 @@ impl FieldElement {
         (!equal_so_far & 1 & greater) as u8
     }
 
-    /// Determine if this `FieldElement` is negative.
+    /// Determine if this `FieldElement` is negative, in the
+    /// sense used in the ed25519 paper.
     ///
     /// # Return
     ///
     /// If negative, return `1i32`.  Otherwise, return `0i32`.
     // XXX should return u8
-    pub fn is_negative(&self) -> i32 { //FeIsNegative
+    pub fn is_negative_ed25519(&self) -> i32 { //FeIsNegative
         let bytes = self.to_bytes();
         (bytes[0] & 1) as i32
     }
 
+    /// Determine if this `FieldElement` is negative, in the
+    /// sense used by Decaf: `x` is nonnegative if the least
+    /// absolute residue for `x` lies in `[0, (p-1)/2]`, and
+    /// is negative otherwise.
+    ///
+    /// # Return
+    ///
+    /// Returns `1u8` if negative, `0u8` if nonnegative.
+    ///
+    /// # Implementation
+    ///
+    /// Uses a trick borrowed from Mike Hamburg's code.  Let `x \in
+    /// F_p` and let `y \in Z` be the least absolute residue for `x`.
+    /// Suppose `y ≤ (p-1)/2`.  Then `2y < p` so `2y = 2y mod p` and
+    /// `2y mod p` is even.  On the other hand, if `y > (p-1)/2` then
+    /// `2y ≥ p`; since `y < p`, `2y \in [p, 2p)`, so `2y mod p =
+    /// 2y-p`, which is odd.
+    ///
+    /// Thus we can test whether `y ≤ (p-1)/2` by checking whether `2y
+    /// mod p` is even.
+    pub fn is_negative_decaf(&self) -> u8 {
+        let y = self + self;
+        (y.to_bytes()[0] & 1) as u8
+    }
+
+    /// Determine if this `FieldElement` is nonnegative, in the
+    /// sense used by Decaf: `x` is nonnegative if the least
+    /// absolute residue for `x` lies in `[0, (p-1)/2]`, and
+    /// is negative otherwise.
+    pub fn is_nonnegative_decaf(&self) -> u8 {
+        1u8 & (!self.is_negative_decaf())
+    }
+
+    /// Determine if this `FieldElement` is zero.
+    ///
+    /// # Return
+    ///
+    /// If zero, return `1u8`.  Otherwise, return `0u8`.
+    pub fn is_zero(&self) -> u8 {
+        return 1u8 & (!self.is_nonzero());
+    }
+
     /// Determine if this `FieldElement` is non-zero.
     ///
     /// # Return
@@ -772,6 +817,48 @@ impl FieldElement {
         t21
     }
 
+    /// Try to compute 1/sqrt(self).
+    ///
+    /// # Return
+    ///
+    /// * If `self` is zero, returns zero.
+    /// * If `self` is square, returns 1/sqrt(self).
+    /// * If `self` is nonsquare, returns `None`.
+    pub fn invsqrt(&self) -> Option<FieldElement> {
+        // We are to compute v as:
+        //     / 1/sqrt(self)  if self is square, nonzero;
+        // v = |            0  if self is zero;
+        //     \     [reject]  if self is nonsquare.
+        //
+        // Using the same trick as in ed25519 decoding, we merge the
+        // inversion, the square root, and the square test as follows.
+        //
+        // To compute sqrt(α), we can compute β = α^((p+3)/8).
+        // Then β^2 = ±α, so multiplying β by sqrt(-1) if necessary
+        // gives sqrt(α).
+        //
+        // To compute 1/sqrt(α), we observe that
+        //    1/β = α^(p-1 - (p+3)/8) = α^((7p-11)/8)
+        //                            = α^3 * (α^7)^((p-5)/8).
+        //
+        // If α is square, then (1/β)^2 = ±(1/α), so that (1/β)^2 α = ±1.
+        let a3 = &self.square() * self;        // α^3
+        let a7 =   &a3.square() * self;        // α^7
+        let mut v = &a3 * &a7.pow_p58();       // α^(p-1-(p+3)/8)
+        let check =  self * &v.square();       // ±1 if α is square
+
+        if v.is_zero() == 1u8 {
+            return Some(v);  // α was zero all along
+        } else if check == FieldElement::one() {
+            return Some(v);  // computed the correct sqrt
+        } else if check == -&FieldElement::one() {
+                             // wrong sign, multiply by sqrt(-1)
+            return Some(&v * &constants::SQRT_M1);
+        } else {
+            return None;     // input was nonsquare
+        }
+    }
+
     /// chi calculates `self^((p-1)/2)`.
     ///
     /// # Return

src/lib.rs

@@ -48,8 +48,12 @@ extern crate rand;
 // Modules for low-level operations directly on field elements and curve points.
 
 pub mod field;
-pub mod curve;
 pub mod scalar;
+pub mod curve;
+
+// Feature gate decaf while our implementation is unfinished and probably incorrect.
+#[cfg(feature = "yolocrypto")]
+pub mod decaf;
 
 // Constant-time functions and other miscelaneous utilities.