Implement Montgomery arithmetic and laddering.

* ADDs part of https://github.com/isislovecruft/curve25519-dalek/issues/47

src/constants.rs

@@ -21,6 +21,7 @@
 use edwards::CompressedEdwardsY;
 #[cfg(feature = "yolocrypto")]
 use decaf::{DecafPoint, DecafBasepointTable};
+use montgomery::CompressedMontgomeryU;
 use scalar::Scalar;
 
 #[cfg(feature="radix_51")]
@@ -52,6 +53,14 @@ pub const BASE_CMPRSSD: CompressedEdwardsY =
                         0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66,
                         0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66]);
 
+/// The X25519 basepoint, in compressed Montgomery form.
+pub const BASE_COMPRESSED_MONTGOMERY: CompressedMontgomeryU =
+    CompressedMontgomeryU([0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+                           0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+                           0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+                           0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]);
+
+
 /// The Ed25519 basepoint, as a `DecafPoint`.  This is called `_POINT` to distinguish it from
 /// `_TABLE`, which provides fast scalar multiplication.
 #[cfg(feature = "yolocrypto")] pub const DECAF_ED25519_BASEPOINT_POINT: DecafPoint =

src/constants_32bit.rs

@@ -75,6 +75,9 @@ pub const HALF: FieldElement32 = FieldElement32([
 pub const A: FieldElement32       = FieldElement32([
     486662, 0, 0, 0, 0, 0, 0, 0, 0, 0, ]);
 
+/// `APLUS2_OVER_FOUR` is (A+2)/4. (This is used internally within Montgomery laddering.)
+pub const APLUS2_OVER_FOUR: FieldElement32 = FieldElement32([121666, 0, 0, 0, 0, 0, 0, 0, 0, 0]);
+
 /// `SQRT_MINUS_A` is sqrt(-486662)
 // XXX I think that this was used in Adam's code for his elligator
 // implementation, but that should maybe be using sqrt(-486664)

src/constants_64bit.rs

@@ -54,6 +54,9 @@ pub const HALF: FieldElement64 = FieldElement64([2251799813685239, 2251799813685
 /// In Montgomery form y² = x³+Ax²+x, Curve25519 has A=486662.
 pub const A: FieldElement64 = FieldElement64([486662, 0, 0, 0, 0]);
 
+/// `APLUS2_OVER_FOUR` is (A+2)/4. (This is used internally within Montgomery laddering.)
+pub const APLUS2_OVER_FOUR: FieldElement64 = FieldElement64([121666, 0, 0, 0, 0]);
+
 /// `SQRT_MINUS_A` is sqrt(-486662)
 // XXX I think that this was used in Adam's code for his elligator
 // implementation, but that should maybe be using sqrt(-486664)

src/edwards.rs

@@ -90,6 +90,7 @@ use constants;
 use field::FieldElement;
 use scalar::Scalar;
 use montgomery::CompressedMontgomeryU;
+use montgomery::MontgomeryPoint;
 
 use subtle::slices_equal;
 use subtle::bytes_equal;
@@ -447,14 +448,14 @@ impl ProjectivePoint {
         CompressedEdwardsY(s)
     }
 
-    /// Convert this point to a `CompressedMontgomeryU`.
+    /// Convert this point to a Montgomery u-coordinate (affine).
     /// Note that this discards the sign.
     ///
     /// # Return
     /// - `None` if `self` is the identity point;
-    /// - `Some(CompressedMontgomeryU)` otherwise.
+    /// - `Some(FieldElement)` otherwise.
     ///
-    pub fn compress_montgomery(&self) -> Option<CompressedMontgomeryU> {
+    fn convert_to_montgomery(&self) -> Option<FieldElement> {
         // u = (1 + y) /  (1 - y)
         // v = sqrt(-486664) * u / x
         //
@@ -470,7 +471,38 @@ impl ProjectivePoint {
         let u = &Z_plus_Y * &Z_minus_Y.invert();
 
         if Z_minus_Y.is_zero() == 0u8 {
-            Some(CompressedMontgomeryU(u.to_bytes()))
+            Some(u)
+        } else {
+            None
+        }
+    }
+
+    /// Convert this point to a `CompressedMontgomeryU`.
+    /// Note that this discards the sign.
+    ///
+    /// # Return
+    /// - `None` if `self` is the identity point;
+    /// - `Some(CompressedMontgomeryU)` otherwise.
+    ///
+    pub fn compress_montgomery(&self) -> Option<CompressedMontgomeryU> {
+        let u: Option<FieldElement> = self.convert_to_montgomery();
+
+        if u.is_some() {
+            Some(CompressedMontgomeryU(u.unwrap().to_bytes()))
+        } else {
+            None
+        }
+    }
+
+    /// Convert this point to its equivalent on the Montgomery form of
+    /// the curve, without compressing.
+    ///
+    /// DOCDOC
+    pub fn to_montgomery(&self) -> Option<MontgomeryPoint> {
+        let u: Option<FieldElement> = self.convert_to_montgomery();
+
+        if u.is_some() {
+            Some(MontgomeryPoint{ U: u.unwrap(), Z: FieldElement::one() })
         } else {
             None
         }
@@ -515,6 +547,11 @@ impl ExtendedPoint {
         }
     }
 
+    /// DOCDOC
+    pub fn to_montgomery(&self) -> Option<MontgomeryPoint> {
+        self.to_projective().to_montgomery()
+    }
+
     /// Compress this point to `CompressedEdwardsY` format.
     pub fn compress_edwards(&self) -> CompressedEdwardsY {
         self.to_projective().compress_edwards()

src/field.rs

@@ -197,11 +197,15 @@ impl FieldElement {
     }
 
     /// Given a nonzero field element, compute its inverse.
+    ///
     /// The inverse is computed as self^(p-2), since
     /// x^(p-2)x = x^(p-1) = 1 (mod p).
-    ///
-    /// XXX should we add a debug_assert that self is nonzero?
+    //
+    // XXX do we want the debug assertion to check for zero? it breaks behaviour
+    // such as that such as in curve25519_dalek::montgomery::test::identity_to_monty.
     pub fn invert(&self) -> FieldElement {
+        // debug_assert!(*self != FieldElement::zero());
+
         // The bits of p-2 = 2^255 -19 -2 are 11010111111...11.
         //
         //                                 nonzero bits of exponent

src/montgomery.rs

@@ -8,7 +8,19 @@
 // - Isis Agora Lovecruft <isis@patternsinthevoid.net>
 // - Henry de Valence <hdevalence@hdevalence.ca>
 
-//! Montgomery arithmetic prototype, subject to revision.
+//! Montgomery arithmetic.
+//!
+//! Apart from the compressed point implementation
+//! (i.e. `CompressedMontgomeryU`), this module is a "clean room" implementation
+//! of the Montgomery arithmetic described in the following papers:
+//!
+//! * Costello, Craig, and Benjamin Smith. "Montgomery curves and their
+//!   arithmetic." Journal of Cryptographic Engineering (2017): 1-14.
+//!   [PDF](http://eprint.iacr.org/2017/212.pdf)
+//!
+//! * Montgomery, Peter L. "Speeding the Pollard and elliptic curve methods of
+//!   factorization." Mathematics of computation 48.177 (1987): 243-264.
+//!   [PDF](http://www.ams.org/mcom/1987-48-177/S0025-5718-1987-0866113-7/)
 
 // We allow non snake_case names because coordinates in projective space are
 // traditionally denoted by the capitalisation of their respective
@@ -16,12 +28,21 @@
 // affine and projective cakes and eat both of them too.
 #![allow(non_snake_case)]
 
+use core::ops::{Mul, MulAssign};
 
 use constants;
 use field::FieldElement;
 use edwards::{ExtendedPoint, CompressedEdwardsY};
+use scalar::Scalar;
 
+// XXX move these to a common "traits" or "group" module? —isis
+use edwards::{Identity, ValidityCheck};
+
+use subtle::slices_equal;
 use subtle::ConditionallyAssignable;
+use subtle::ConditionallySwappable;
+use subtle::Equal;
+use subtle::Mask;
 
 /// In "Montgomery u" format, as used in X25519, a point `(u,v)` on
 /// the Montgomery curve
@@ -40,6 +61,11 @@ pub struct CompressedMontgomeryU(pub [u8; 32]);
 
 impl CompressedMontgomeryU {
     /// View this `CompressedMontgomeryU` as an array of bytes.
+    pub fn as_bytes<'a>(&'a self) -> &'a [u8; 32] {
+        &self.0
+    }
+
+    /// Convert this `CompressedMontgomeryU` to an array of bytes.
     pub fn to_bytes(&self) -> [u8; 32] {
         self.0
     }
@@ -65,7 +91,7 @@ impl CompressedMontgomeryU {
     /// * `v` is not square.
     //
     // XXX any other exceptional points for the birational map?
-    pub fn decompress(&self) -> Option<ExtendedPoint> {
+    pub fn decompress_edwards(&self) -> Option<ExtendedPoint> {
         let u: FieldElement = FieldElement::from_bytes(&self.0);
 
         // If u = -1, then v^2 = u*(u^2+486662*u+1) = 486660.
@@ -84,6 +110,23 @@ impl CompressedMontgomeryU {
         CompressedEdwardsY(y.to_bytes()).decompress()
     }
 
+    /// Decompress this `CompressedMontgomeryU` to a `MontgomeryPoint`.
+    ///
+    /// Going from affine to projective coordinates, we have:
+    ///
+    ///      u → U/W
+    ///
+    /// # Returns
+    ///
+    /// A projective `MontgomeryPoint` corresponding to this compressed point.
+    pub fn decompress_montgomery(&self) -> MontgomeryPoint {
+        MontgomeryPoint{
+            // XXX is it a problem here if we're not using a canonical encoding? —isis
+            U: FieldElement::from_bytes(&self.0),
+            W: FieldElement::one(),
+        }
+    }
+
     /// Given a Montgomery `u` coordinate, compute an Edwards `y` via
     /// `y = (u-1)/(u+1)`.
     ///
@@ -150,33 +193,246 @@ impl CompressedMontgomeryU {
     }
 }
 
+/// A point on the Montgomery form of the curve, in projective 𝗣^2 coordinates.
+///
+/// The transition between affine and projective is given by
+///
+///      u → U/W
+///      v → V/W
+///
+/// thus the Montgomery curve equation
+///
+///      E_(A,B) : Bv² = u(u² + Au + 1)
+///
+/// becomes
+///
+///      E_(A,B) : BV²W = U(U² + AUW + W²) ⊆ 𝗣^2
+///
+/// Here, again, to differentiate from points in the twisted Edwards model, we
+/// call the point `(x,y)` in affine coordinates `(u,v)` and similarly in projective
+/// space we use `(U:V:W)`.  However, since (as per Montgomery's original work) the
+/// v-coordinate is superfluous to the definition of the group law, we merely
+/// use `(U:W)`.
+#[derive(Copy, Clone, Debug)]
+#[allow(missing_docs)]
+pub struct MontgomeryPoint{
+    pub U: FieldElement,
+    pub W: FieldElement,
+}
+
+/// The identity point is a unique point (the only where `W = 0`) on the curve.
+///
+/// In projective coordinates, the quotient map `x : E (A,B) → E/<⦵> = 𝗣¹` is
+///
+///               ⎧ (x_P:1) if P = (x_P:y_P:1) ,
+///      x : P ↦  ⎨
+///               ⎩   (1:0) if P = O = (0:1:0) .
+///
+/// We emphasize that the formula `x((U: V : W)) = (U : W)` only holds on the
+/// open subset of `E_(A,B)` where `W ≠ 0`; it does not extend to the point
+/// `O = (0:1:0)` at infinity, because `(0:0)` is not a projective point.
+///
+/// # Returns
+///
+/// The (exceptional) point at infinity in the Montgomery model.
+impl Identity for MontgomeryPoint {
+    fn identity() -> MontgomeryPoint {
+        MontgomeryPoint {
+            U: FieldElement::one(),
+            W: FieldElement::zero(),
+        }
+    }
+}
+
+/// Determine if two `MontgomeryPoint`s are equal, in constant time.
+///
+/// # Note
+///
+/// Because a compressed point on the Montgomery form of the curve doesn't
+/// include the sign bit, there's two points here (if translated from the
+/// Edwards form) which will equate.
+///
+/// # Returns
+///
+/// `1` if the points are equal, and `0` otherwise.
+impl Equal for MontgomeryPoint {
+    fn ct_eq(&self, that: &MontgomeryPoint) -> u8 {
+        slices_equal(self.compress_montgomery().as_bytes(),
+                     that.compress_montgomery().as_bytes())
+    }
+}
+
+/// Determine if this `MontgomeryPoint` is valid.
+///
+/// # Note
+///
+/// All points, except for `(X:W) = (0:0)`, are valid, since the projective
+/// model is linear through the origin and is comprised by all `X` in
+/// ℤ/(2²⁵⁵-19).
+///
+/// # Returns
+///
+/// `true` if it is valid, and `false` otherwise.
+impl ValidityCheck for MontgomeryPoint {
+    fn is_valid(&self) -> bool {
+        let zero = FieldElement::zero();
+
+        if (self.U.ct_eq(&zero) & self.W.ct_eq(&zero)) == 1 {
+            return true;
+        }
+        false
+    }
+}
+
+/// Conditionally assign another `MontgomeryPoint` to this point, in constant time.
+///
+/// If `choice == 1`, assign `that` to `self`.  Otherwise, leave `self`
+/// unchanged.
+impl ConditionallyAssignable for MontgomeryPoint {
+    fn conditional_assign(&mut self, that: &MontgomeryPoint, choice: Mask) {
+        self.U.conditional_assign(&that.U, choice);
+        self.W.conditional_assign(&that.W, choice);
+    }
+}
+
+impl MontgomeryPoint {
+    /// Compress this point to only its u-coordinate (note: affine).
+    ///
+    /// # Returns
+    ///
+    /// A `CompressedMontgomeryU`.
+    pub fn compress_montgomery(&self) -> CompressedMontgomeryU {
+        let u_affine: FieldElement = &self.U * &self.W.invert();
+
+        CompressedMontgomeryU(u_affine.to_bytes())
+    }
+
+    /// Differential addition for single-coordinate Montgomery points.
+    ///
+    /// Montgomery coordinates in projective 𝗣¹ space are odd in that 𝗣¹
+    /// inherits none of the group structure from E_(A,B).  Hence, the mapping
+    /// of the group operation, `⊕`, is undefined for the pair `(x(P), x(Q))`;
+    /// that is, given `x(P)` and `x(Q)`, we cannot derive `x(P ⊕ Q)`.  This is
+    /// due to the fact that, in Montgomery coordinates, `x(P)` determines `P`
+    /// only up to a sign, and thus we cannot differentiate `x(P ⊕ Q)` from
+    /// `x(P ⊖ Q)`.  However, via differential addition, any three of the values
+    /// `{x(P), x(Q), x(P ⊕ Q), x(P ⊖ Q)}` determines the forth, so we can
+    /// define *pseudo-addition* for a singular coordinate.
+    ///
+    /// # Warning
+    ///
+    /// If the `difference` is the identity point, or a two torsion point, the
+    /// results of this method are not correct, but instead result in `(0:0)`
+    /// (an invalid projective point in the Montgomery model).
+    ///
+    // XXX API-wise, do we care that doubling is degenerate, or should we allow
+    //     the user to do a stupid and inefficient (albeit not incorrect) thing?
+    fn differential_add(&self, that: &MontgomeryPoint,
+                        difference: &MontgomeryPoint) -> MontgomeryPoint {
+        // debug_assert!(self.ct_eq(that) != 1);     // The doubling case is degenerate
+        // debug_assert!(!difference.is_identity()); // P ⦵ Q ∉ {O,T}
+        // debug_assert!(!difference.is_two_torsion_point());
+
+        let v1: FieldElement = &(&self.U + &self.W) * &(&that.U - &that.W);
+        let v2: FieldElement = &(&self.U - &self.W) * &(&that.U + &that.W);
+
+        MontgomeryPoint {
+            U: &difference.W * &(&v1 + &v2).square(),  // does reduction on square()
+            W: &difference.U * &(&v1 - &v2).square(),  // does reduction on square()
+        }
+    }
+
+    /// Differential doubling for single-coordinate Montgomery points.
+    ///
+    /// DOCDOC
+    ///
+    /// # Returns
+    ///
+    /// A Montgomery point.
+    fn differential_double(&self) -> MontgomeryPoint {
+        let mut v1: FieldElement;
+        let v2: FieldElement;
+        let v3: FieldElement;
+
+        v1 = (&self.U + &self.W).square();
+        v2 = (&self.U - &self.W).square();
+
+        let U: FieldElement = &v1 * &v2;
+
+        v1 -= &v2;
+        v3  = &(&constants::APLUS2_OVER_FOUR * &v1) + &v2;
+
+        let W: FieldElement = &v1 * &v3;
+
+        MontgomeryPoint{ U: U, W: W }
+    }
+}
+
+/// Multiply this `MontgomeryPoint` by a `Scalar`.
+///
+/// DOCDOC
+/// explain montgomery laddering
+impl<'a, 'b> Mul<&'b Scalar> for &'a MontgomeryPoint {
+    type Output = MontgomeryPoint;
+
+    fn mul(self, scalar: &'b Scalar) -> MontgomeryPoint {
+        let mut x0: MontgomeryPoint = MontgomeryPoint::identity();
+        let mut x1: MontgomeryPoint = *self;
+
+        let bits: [i8; 256] = scalar.bits();
+
+        for i in (0..255).rev() {
+            let mask: u8 = (bits[i+1] ^ bits[i]) as u8;
+
+            debug_assert!(mask == 0 || mask == 1);
+
+            x0.conditional_swap(&mut x1, mask);
+            x1 = x0.differential_add(&x1, &self);
+            x0 = x0.differential_double();
+        }
+        x0.conditional_swap(&mut x1, bits[0] as u8);
+        x0
+    }
+}
+
+impl<'b> MulAssign<&'b Scalar> for MontgomeryPoint {
+    fn mul_assign(&mut self, scalar: &'b Scalar) {
+        let result = (self as &MontgomeryPoint) * scalar;
+        *self = result;
+    }
+}
+
+impl<'a, 'b> Mul<&'b MontgomeryPoint> for &'a Scalar {
+    type Output = MontgomeryPoint;
+
+    fn mul(self, point: &'b MontgomeryPoint) -> MontgomeryPoint {
+        point * &self
+    }
+}
+
 // ------------------------------------------------------------------------
 // Tests
 // ------------------------------------------------------------------------
 
 #[cfg(test)]
 mod test {
+    use constants::BASE_COMPRESSED_MONTGOMERY;
     use edwards::Identity;
     use super::*;
 
-    /// The X25519 basepoint, in compressed Montgomery form.
-    static BASE_CMPRSSD_MONTY: CompressedMontgomeryU =
-        CompressedMontgomeryU([0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-                               0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-                               0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-                               0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]);
+    use rand::OsRng;
 
     /// Test Montgomery conversion against the X25519 basepoint.
     #[test]
     fn basepoint_to_montgomery() {
         assert_eq!(constants::ED25519_BASEPOINT_POINT.compress_montgomery().unwrap(),
-                   BASE_CMPRSSD_MONTY);
+                   BASE_COMPRESSED_MONTGOMERY);
     }
 
     /// Test Montgomery conversion against the X25519 basepoint.
     #[test]
     fn basepoint_from_montgomery() {
-        assert_eq!(BASE_CMPRSSD_MONTY.decompress().unwrap().compress_edwards(),
+        assert_eq!(BASE_COMPRESSED_MONTGOMERY.decompress_edwards().unwrap().compress_edwards(),
                    constants::BASE_CMPRSSD);
     }
 
@@ -189,7 +445,7 @@ mod test {
         let minus_one = FieldElement::minus_one();
         let minus_one_bytes = minus_one.to_bytes();
         let div_by_zero_u = CompressedMontgomeryU(minus_one_bytes);
-        assert!(div_by_zero_u.decompress().is_none());
+        assert!(div_by_zero_u.decompress_edwards().is_none());
     }
 
     /// Montgomery compression of the identity point should
@@ -199,4 +455,60 @@ mod test {
         let id = ExtendedPoint::identity();
         assert!(id.compress_montgomery().is_none());
     }
+
+    #[test]
+    fn projective_to_affine_roundtrips() {
+        let p = BASE_COMPRESSED_MONTGOMERY.decompress_montgomery();
+
+    }
+
+    #[test]
+    fn differential_double_matches_double() {
+        let p: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double();
+        let q: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress_montgomery().differential_double();
+
+        assert_eq!(p.compress_montgomery().unwrap(), q.compress_montgomery());
+    }
+
+    #[test]
+    fn differential_add_matches_edwards_model() {
+        let mut csprng: OsRng = OsRng::new().unwrap();
+
+        let s1: Scalar = Scalar::random(&mut csprng);
+        let s2: Scalar = Scalar::random(&mut csprng);
+        let p1: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s1;
+        let p2: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s2;
+        let diff: ExtendedPoint = &p1 - &p2;
+
+        let p1m: MontgomeryPoint = p1.to_montgomery().unwrap();
+        let p2m: MontgomeryPoint = p2.to_montgomery().unwrap();
+        let diffm: MontgomeryPoint = diff.to_montgomery().unwrap();
+
+        let result = p1m.differential_add(&p2m, &diffm);
+
+        assert_eq!(result.compress_montgomery(), (&p1 + &p2).compress_montgomery().unwrap());
+    }
+
+    #[test]
+    fn ladder_matches_scalarmult() {
+        let mut csprng: OsRng = OsRng::new().unwrap();
+
+        let s: Scalar = Scalar::random(&mut csprng);
+        let p_edwards: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s;
+        let p_montgomery: MontgomeryPoint = p_edwards.to_montgomery().unwrap();
+
+        let expected = &s * &p_edwards;
+        let result   = &s * &p_montgomery;
+
+        assert_eq!(result.compress_montgomery(), expected.compress_montgomery().unwrap())
+    }
+
+    #[test]
+    fn ladder_basepoint_times_two_matches_double() {
+        let two: Scalar = Scalar::from_u64(2u64);
+        let result: MontgomeryPoint = &BASE_COMPRESSED_MONTGOMERY.decompress_montgomery() * &two;
+        let mut expected: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double();
+
+        assert_eq!(result.compress_montgomery(), expected.compress_montgomery().unwrap());
+    }
 }