Merge branch 'feature/montgomery-arithmetic_r1' into develop

Cargo.toml

@@ -33,7 +33,7 @@ version = "0.3"
 version = "0.6"
 
 [dependencies.subtle]
-version = "^0.2"
+version = "^0.3"
 default-features = false
 
 [dependencies.generic-array]

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 the Montgomery ladder.)
+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 the Montgomery ladder.)
+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/decaf.rs

@@ -189,7 +189,7 @@ impl<'de> Deserialize<'de> for DecafPoint {
 
 /// A point in a prime-order group.
 ///
-/// XXX think about how this API should work
+// XXX think about how this API should work
 #[derive(Copy, Clone)]
 pub struct DecafPoint(pub ExtendedPoint);
 
@@ -752,7 +752,7 @@ mod test {
     fn decaf_decompress_id() {
         let compressed_id = CompressedDecaf::identity();
         let id = compressed_id.decompress().unwrap();
-        assert_eq!(id.0.compress_edwards(), CompressedEdwardsY::identity());
+        assert_eq!(id.0.compress(), CompressedEdwardsY::identity());
     }
 
     #[test]
@@ -768,7 +768,7 @@ mod test {
         // Check that bp_recaf differs from bp by a point of order 4
         let diff = &constants::ED25519_BASEPOINT_POINT - &bp_recaf;
         let diff4 = diff.mult_by_pow_2(4); // XXX this is wrong
-        assert_eq!(diff4.compress_edwards(), CompressedEdwardsY::identity());
+        assert_eq!(diff4.compress(), CompressedEdwardsY::identity());
     }
 
     #[test]

src/edwards.rs

@@ -89,7 +89,7 @@ use core::ops::Index;
 use constants;
 use field::FieldElement;
 use scalar::Scalar;
-use montgomery::CompressedMontgomeryU;
+use montgomery::MontgomeryPoint;
 
 use subtle::slices_equal;
 use subtle::bytes_equal;
@@ -123,7 +123,6 @@ impl CompressedEdwardsY {
     }
 
     /// Copy this `CompressedEdwardsY` to an array of bytes.
-    /// XXX is this useful?
     pub fn to_bytes(&self) -> [u8; 32] {
         self.0
     }
@@ -169,7 +168,7 @@ impl Serialize for ExtendedPoint {
     fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
         where S: Serializer
     {
-        serializer.serialize_bytes(self.compress_edwards().as_bytes())
+        serializer.serialize_bytes(self.compress().as_bytes())
     }
 }
 
@@ -393,8 +392,8 @@ impl ConditionallyAssignable for ExtendedPoint {
 
 impl Equal for ExtendedPoint {
     fn ct_eq(&self, other: &ExtendedPoint) -> u8 {
-        slices_equal(self.compress_edwards().as_bytes(),
+        slices_equal(self.compress().as_bytes(),
-                     other.compress_edwards().as_bytes())
+                     other.compress().as_bytes())
     }
 }
 
@@ -436,7 +435,7 @@ impl ProjectivePoint {
     }
 
     /// Convert this point to a `CompressedEdwardsY`
-    pub fn compress_edwards(&self) -> CompressedEdwardsY {
+    pub fn compress(&self) -> CompressedEdwardsY {
         let recip = self.Z.invert();
         let x = &self.X * &recip;
         let y = &self.Y * &recip;
@@ -447,32 +446,67 @@ impl ProjectivePoint {
         CompressedEdwardsY(s)
     }
 
-    /// Convert this point to a `CompressedMontgomeryU`.
+    /// Convert this projective point in the Edwards model to its equivalent
-    /// Note that this discards the sign.
+    /// projective point on the Montgomery form of the curve.
     ///
-    /// # Return
+    /// Taking the Montgomery curve equation in affine coordinates:
-    /// - `None` if `self` is the identity point;
-    /// - `Some(CompressedMontgomeryU)` otherwise.
     ///
-    pub fn compress_montgomery(&self) -> Option<CompressedMontgomeryU> {
+    ///     E_(A,B) = Bv² = u³ + Au² + u   <span style="float: right">(1)</span>
-        // u = (1 + y) /  (1 - y)
+    ///
-        // v = sqrt(-486664) * u / x
+    /// and given its relations to the coordinates of the Edwards model:
-        //
+    ///
-        // since y = Y/Z, x = X/Z,
+    ///     u = (1+y)/(1-y)                <span style="float: right">(2)</span>
-        //
+    ///     v = (λu)/(x)
-        // u = (1 + Y/Z) / (1 - Y/Z);
+    ///
-        //   =   (Z + Y) / (Z - Y);
+    /// Converting from affine to projective coordinates in the Montgomery
-        //
+    /// model, we arrive at:
-        // exceptional points:
+    ///
-        // y = 1 <=> Y/Z = 1 <=> Z - Y = 0
+    ///     u = (Z+Y)/(Z-Y)                <span style="float: right">(3)</span>
-        let Z_plus_Y   = &self.Z + &self.Y;
+    ///     v = λ * ((Z+Y)/(Z-Y)) * (Z/X)
-        let Z_minus_Y  = &self.Z - &self.Y;
+    ///
-        let u = &Z_plus_Y * &Z_minus_Y.invert();
+    /// The transition between affine and projective is given by
-
+    ///
-        if Z_minus_Y.is_zero() == 0u8 {
+    ///      u → U/W                       <span style="float: right">(4)</span>
-            Some(CompressedMontgomeryU(u.to_bytes()))
+    ///      v → V/W
-        } else {
+    ///
-            None
+    /// thus the Montgomery curve equation (1) becomes
+    ///
+    ///      E_(A,B) : BV²W = U³ + AU²W + UW² ⊆ 𝗣^2  <span style="float: right">(5)</span>
+    ///
+    /// 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)`.
+    ///
+    /// Therefore, the direct translation between projective Montgomery points
+    /// and projective twisted Edwards points is
+    ///
+    ///      (U:W) = (Z+Y:Z-Y)             <span style="float: right">(6)</span>
+    ///
+    /// Note, however, that there appears to be an exception where `Z=Y`,
+    /// since—from equation 2—this would imply that `y=1` (thus causing the
+    /// denominator to be zero).  If this is the case, then it follows from the
+    /// twisted Edwards curve equation
+    ///
+    ///      -x² + y² = 1 + dx²y²          <span style="float: right">(7)</span>
+    ///
+    /// that
+    ///
+    ///      -x² + 1 = 1 + dx²
+    ///
+    /// and, assuming that `d ≠ -1`,
+    ///
+    ///      -x² = x²
+    ///       x  = 0
+    ///
+    /// Therefore, the only valid point with `y=1` is the twisted Edwards
+    /// identity point, which correctly becomes `(1:0)`, that is, the identity,
+    /// in the Montgomery model.
+    pub fn to_montgomery(&self) -> MontgomeryPoint {
+        MontgomeryPoint{
+            U: &self.Z + &self.Y,
+            W: &self.Z - &self.Y,
         }
     }
 }
@@ -515,20 +549,15 @@ impl ExtendedPoint {
         }
     }
 
-    /// Compress this point to `CompressedEdwardsY` format.
+    /// Convert this point to its equivalent on the Montgomery form of the
-    pub fn compress_edwards(&self) -> CompressedEdwardsY {
+    /// curve.
-        self.to_projective().compress_edwards()
+    pub fn to_montgomery(&self) -> MontgomeryPoint {
+        self.to_projective().to_montgomery()
     }
 
-    /// Convert this point to a `CompressedMontgomeryU`.
+    /// Compress this point to `CompressedEdwardsY` format.
-    /// Note that this discards the sign.
+    pub fn compress(&self) -> CompressedEdwardsY {
-    ///
+        self.to_projective().compress()
-    /// # Return
-    /// - `None` if `self` is the identity point;
-    /// - `Some(CompressedMontgomeryU)` otherwise.
-    ///
-    pub fn compress_montgomery(&self) -> Option<CompressedMontgomeryU> {
-        self.to_projective().compress_montgomery()
     }
 }
 
@@ -1305,7 +1334,7 @@ mod test {
         assert!(bp.is_valid());
         // Check that decompression actually gives the correct X coordinate
         assert_eq!(base_X, bp.X);
-        assert_eq!(bp.compress_edwards(), constants::BASE_CMPRSSD);
+        assert_eq!(bp.compress(), constants::BASE_CMPRSSD);
     }
 
     /// Test sign handling in decompression
@@ -1328,7 +1357,7 @@ mod test {
     #[test]
     fn basepoint_mult_one_vs_basepoint() {
         let bp = &constants::ED25519_BASEPOINT_TABLE * &Scalar::one();
-        let compressed = bp.compress_edwards();
+        let compressed = bp.compress();
         assert_eq!(compressed, constants::BASE_CMPRSSD);
     }
 
@@ -1336,7 +1365,7 @@ mod test {
     #[test]
     fn basepoint_table_basepoint_function_correct() {
         let bp = constants::ED25519_BASEPOINT_TABLE.basepoint();
-        assert_eq!(bp.compress_edwards(), constants::BASE_CMPRSSD);
+        assert_eq!(bp.compress(), constants::BASE_CMPRSSD);
     }
 
     /// Test `impl Add<ExtendedPoint> for ExtendedPoint`
@@ -1345,7 +1374,7 @@ mod test {
     fn basepoint_plus_basepoint_vs_basepoint2() {
         let bp = constants::ED25519_BASEPOINT_POINT;
         let bp_added = &bp + &bp;
-        assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD);
+        assert_eq!(bp_added.compress(), BASE2_CMPRSSD);
     }
 
     /// Test `impl Add<ProjectiveNielsPoint> for ExtendedPoint`
@@ -1354,7 +1383,7 @@ mod test {
     fn basepoint_plus_basepoint_projective_niels_vs_basepoint2() {
         let bp = constants::ED25519_BASEPOINT_POINT;
         let bp_added = (&bp + &bp.to_projective_niels()).to_extended();
-        assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD);
+        assert_eq!(bp_added.compress(), BASE2_CMPRSSD);
     }
 
     /// Test `impl Add<AffineNielsPoint> for ExtendedPoint`
@@ -1364,7 +1393,7 @@ mod test {
         let bp = constants::ED25519_BASEPOINT_POINT;
         let bp_affine_niels = bp.to_affine_niels();
         let bp_added = (&bp + &bp_affine_niels).to_extended();
-        assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD);
+        assert_eq!(bp_added.compress(), BASE2_CMPRSSD);
     }
 
     /// Check that equality of `ExtendedPoints` handles projective
@@ -1389,15 +1418,15 @@ mod test {
         let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
         let aB_affine_niels = aB.to_affine_niels();
         let also_aB = (&ExtendedPoint::identity() + &aB_affine_niels).to_extended();
-        assert_eq!(     aB.compress_edwards(),
+        assert_eq!(     aB.compress(),
-                   also_aB.compress_edwards());
+                   also_aB.compress());
     }
 
     /// Test basepoint_mult versus a known scalar multiple from ed25519.py
     #[test]
     fn basepoint_mult_vs_ed25519py() {
         let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
-        assert_eq!(aB.compress_edwards(), A_TIMES_BASEPOINT);
+        assert_eq!(aB.compress(), A_TIMES_BASEPOINT);
     }
 
     /// Test that multiplication by the basepoint order kills the basepoint
@@ -1415,20 +1444,20 @@ mod test {
         let table = EdwardsBasepointTable::create(&constants::ED25519_BASEPOINT_POINT);
         let aB_1 = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
         let aB_2 = &table * &A_SCALAR;
-        assert_eq!(aB_1.compress_edwards(), aB_2.compress_edwards());
+        assert_eq!(aB_1.compress(), aB_2.compress());
     }
 
     /// Test scalar_mult versus a known scalar multiple from ed25519.py
     #[test]
     fn scalar_mult_vs_ed25519py() {
         let aB = &constants::ED25519_BASEPOINT_POINT * &A_SCALAR;
-        assert_eq!(aB.compress_edwards(), A_TIMES_BASEPOINT);
+        assert_eq!(aB.compress(), A_TIMES_BASEPOINT);
     }
 
     /// Test basepoint.double() versus the 2*basepoint constant.
     #[test]
     fn basepoint_double_vs_basepoint2() {
-        assert_eq!(constants::ED25519_BASEPOINT_POINT.double().compress_edwards(),
+        assert_eq!(constants::ED25519_BASEPOINT_POINT.double().compress(),
                    BASE2_CMPRSSD);
     }
 
@@ -1437,14 +1466,14 @@ mod test {
     fn basepoint_mult_two_vs_basepoint2() {
         let mut two_bytes = [0u8; 32]; two_bytes[0] = 2;
         let bp2 = &constants::ED25519_BASEPOINT_TABLE * &Scalar(two_bytes);
-        assert_eq!(bp2.compress_edwards(), BASE2_CMPRSSD);
+        assert_eq!(bp2.compress(), BASE2_CMPRSSD);
     }
 
     /// Check that converting to projective and then back to extended round-trips.
     #[test]
     fn basepoint_projective_extended_round_trip() {
         assert_eq!(constants::ED25519_BASEPOINT_POINT
-                       .to_projective().to_extended().compress_edwards(),
+                       .to_projective().to_extended().compress(),
                    constants::BASE_CMPRSSD);
     }
 
@@ -1452,7 +1481,7 @@ mod test {
     #[test]
     fn basepoint16_vs_mult_by_pow_2_4() {
         let bp16 = constants::ED25519_BASEPOINT_POINT.mult_by_pow_2(4);
-        assert_eq!(bp16.compress_edwards(), BASE16_CMPRSSD);
+        assert_eq!(bp16.compress(), BASE16_CMPRSSD);
     }
 
     /// Test that the conditional assignment trait works for AffineNielsPoints.
@@ -1480,7 +1509,7 @@ mod test {
 
     #[test]
     fn compressed_identity() {
-        assert_eq!(ExtendedPoint::identity().compress_edwards(),
+        assert_eq!(ExtendedPoint::identity().compress(),
                    CompressedEdwardsY::identity());
     }
 
@@ -1518,7 +1547,7 @@ mod test {
         let P1 = &G * &s;
         let P2 = &s * &G;
 
-        assert!(P1.compress_edwards().to_bytes() == P2.compress_edwards().to_bytes());
+        assert!(P1.compress().to_bytes() == P2.compress().to_bytes());
     }
 
     #[test]
@@ -1542,7 +1571,7 @@ mod test {
         fn double_scalar_mult_basepoint_vs_ed25519py() {
             let A = A_TIMES_BASEPOINT.decompress().unwrap();
             let result = vartime::double_scalar_mult_basepoint(&A_SCALAR, &A, &B_SCALAR);
-            assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT);
+            assert_eq!(result.compress(), DOUBLE_SCALAR_MULT_RESULT);
         }
 
         #[test]
@@ -1552,7 +1581,7 @@ mod test {
                 &[A_SCALAR, B_SCALAR],
                 &[A, constants::ED25519_BASEPOINT_POINT]
             );
-            assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT);
+            assert_eq!(result.compress(), DOUBLE_SCALAR_MULT_RESULT);
         }
 
         #[test]
@@ -1567,7 +1596,7 @@ mod test {
                 &[A, constants::ED25519_BASEPOINT_POINT]
             );
 
-            assert_eq!(result_vartime.compress_edwards(), result_consttime.compress_edwards());
+            assert_eq!(result_vartime.compress(), result_consttime.compress());
         }
     }
 
@@ -1579,7 +1608,7 @@ mod test {
     fn serde_cbor_basepoint_roundtrip() {
         let output = serde_cbor::to_vec(&constants::ED25519_BASEPOINT_POINT).unwrap();
         let parsed: ExtendedPoint = serde_cbor::from_slice(&output).unwrap();
-        assert_eq!(parsed.compress_edwards(), constants::BASE_CMPRSSD);
+        assert_eq!(parsed.compress(), constants::BASE_CMPRSSD);
     }
 
     #[test]
@@ -1615,7 +1644,7 @@ mod bench {
     #[bench]
     fn edwards_compress(b: &mut Bencher) {
         let B = &constants::ED25519_BASEPOINT_POINT;
-        b.iter(|| B.compress_edwards());
+        b.iter(|| B.compress());
     }
 
     #[bench]

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,22 @@
 // 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 "group" module?  At the same time, we should
+// XXX probably make a `trait Group` once const generics are implemented in
+// XXX Rust. —isis
+use edwards::{Identity, ValidityCheck};
 
 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
@@ -33,13 +55,16 @@ use subtle::ConditionallyAssignable;
 /// coordinates.  For Montgomery curves, it is possible to compute the
 /// `u`-coordinate of `n(u,v)` just from `n` and `u`, so it is not
 /// necessary to use `v` for a Diffie-Hellman key exchange.
-///
-/// XXX add note on monty, twist security, edwards impl of x25519, rfc7748
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
 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 +90,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 +109,22 @@ 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(&self) -> MontgomeryPoint {
+        MontgomeryPoint{
+            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,34 +191,259 @@ 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 for the purposes of scalar multiplication, 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 {
+        // (U_P:W_P) = (U_Q:W_Q) iff U_P * W_Q == U_Q * W_P,
+        // since U_P/W_P == U_Q/W_Q.
+        (&self.U * &that.W).ct_eq(&(&self.W * &that.U))
+    }
+}
+
+/// Determine if this `MontgomeryPoint` is valid.
+///
+/// # Note
+///
+/// All projective 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), thus `(0:0)` is the only element in Fₚ² which is not a
+/// projective point.
+///
+/// # 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(&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).
+    ///
+    /// The doubling case is degenerate, in that `P ⦵ Q ∉ {O,T}`, where `T` is
+    /// the two torsion point.
+    fn differential_add(&self, that: &MontgomeryPoint,
+                        difference: &MontgomeryPoint) -> MontgomeryPoint {
+        // XXX Do we want these debug assertions? We would need to implement
+        // XXX is_two_torsion_point(). —isis
+        // 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()
+        }
+    }
+
+    /// Pseudo-doubling for single-coordinate Montgomery points.
+    ///
+    /// Given a Montgomery U-coordinate of a point `P`, compute the
+    /// U-coordinate given by
+    ///
+    ///      differential_double: x(P) ⟼ x([2]P)
+    ///
+    /// # Returns
+    ///
+    /// A Montgomery point equal to doubling this one.
+    ///
+    // XXX It seems possible that combining the differential_add() and
+    // XXX differential_double() methods would save a non-trivial amount of
+    // XXX computation in the ladder. —isis
+    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`.
+///
+/// The reader is refered to §5.3 of ["Montgomery Curves and Their Arithmetic"
+/// by Craig Costello and Benjamin Smith](https://eprint.iacr.org/2017/212.pdf)
+/// for an overview of side-channel-free Montgomery laddering algorithms.
+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::ED25519_BASEPOINT_TABLE;
+    use constants::BASE_COMPRESSED_MONTGOMERY;
     use edwards::Identity;
     use super::*;
 
-    /// The X25519 basepoint, in compressed Montgomery form.
+    use rand::OsRng;
-    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]);
 
     /// Test Montgomery conversion against the X25519 basepoint.
     #[test]
     fn basepoint_to_montgomery() {
-        assert_eq!(constants::ED25519_BASEPOINT_POINT.compress_montgomery().unwrap(),
+        assert_eq!(constants::ED25519_BASEPOINT_POINT.to_montgomery().compress(),
-                   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,
-                   constants::BASE_CMPRSSD);
+                   constants::BASE_CMPRSSD.decompress().unwrap().to_montgomery().compress());
     }
 
     /// If u = -1, then v^2 = u*(u^2+486662*u+1) = 486660.
@@ -189,14 +455,146 @@ 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
+    /// Montgomery compression of the identity point should not fail (since the
-    /// fail (it's sent to infinity).
+    /// mapping in `ProjectivePoint.to_montgomery()` should be valid for the
+    /// identity.
     #[test]
     fn identity_to_monty() {
         let id = ExtendedPoint::identity();
-        assert!(id.compress_montgomery().is_none());
+        assert_eq!(id.to_montgomery().compress(), MontgomeryPoint::identity().compress());
+    }
+
+    #[test]
+    fn projective_to_affine_roundtrips() {
+        assert_eq!(BASE_COMPRESSED_MONTGOMERY.decompress().compress(),
+                   BASE_COMPRESSED_MONTGOMERY);
+
+    }
+
+    #[test]
+    fn differential_double_matches_double() {
+        let p: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double();
+        let q: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress().differential_double();
+
+        assert_eq!(p.to_montgomery().compress(), q.compress());
+    }
+
+    #[test]
+    fn montgomery_ct_eq_ne() {
+        let mut csprng: OsRng = OsRng::new().unwrap();
+        let s1: Scalar = Scalar::random(&mut csprng);
+        let s2: Scalar = Scalar::random(&mut csprng);
+        let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery();
+        let p2: MontgomeryPoint = (&s2 * &ED25519_BASEPOINT_TABLE).to_montgomery();
+
+        assert_eq!(p1.ct_eq(&p2), 0);
+    }
+
+    #[test]
+    fn montgomery_ct_eq_eq() {
+        let mut csprng: OsRng = OsRng::new().unwrap();
+        let s1: Scalar = Scalar::random(&mut csprng);
+        let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery();
+
+        assert_eq!(p1.ct_eq(&p1), 1);
+    }
+
+    #[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();
+        let p2m: MontgomeryPoint = p2.to_montgomery();
+        let diffm: MontgomeryPoint = diff.to_montgomery();
+
+        let result = p1m.differential_add(&p2m, &diffm);
+
+        assert_eq!(result.compress(), (&p1 + &p2).to_montgomery().compress());
+    }
+
+    #[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();
+
+        let expected = &s * &p_edwards;
+        let result   = &s * &p_montgomery;
+
+        assert_eq!(result.compress(), expected.to_montgomery().compress())
+    }
+
+    #[test]
+    fn ladder_basepoint_times_two_matches_double() {
+        let two: Scalar = Scalar::from_u64(2u64);
+        let result: MontgomeryPoint = &BASE_COMPRESSED_MONTGOMERY.decompress() * &two;
+        let expected: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double();
+
+        assert_eq!(result.compress(), expected.to_montgomery().compress());
+    }
+
+    #[test]
+    #[should_panic(expected = "assertion failed: self[31] <= 127")]
+    fn ladder_matches_scalarmult_with_scalar_high_bit_set() {
+        let mut s: Scalar = Scalar::one();
+
+        s[31] = 255;
+
+        let result: MontgomeryPoint = &BASE_COMPRESSED_MONTGOMERY.decompress() * &s;
+        let expected: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s;
+
+        assert_eq!(result.compress(), expected.to_montgomery().compress())
+    }
+}
+
+#[cfg(all(test, feature = "bench"))]
+mod bench {
+    use rand::OsRng;
+    use constants::ED25519_BASEPOINT_TABLE;
+    use constants::BASE_COMPRESSED_MONTGOMERY;
+    use test::Bencher;
+    use super::*;
+
+    #[bench]
+    fn montgomery_ct_eq(b: &mut Bencher) {
+        let mut csprng: OsRng = OsRng::new().unwrap();
+        let s1: Scalar = Scalar::random(&mut csprng);
+        let s2: Scalar = Scalar::random(&mut csprng);
+        let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery();
+        let p2: MontgomeryPoint = (&s2 * &ED25519_BASEPOINT_TABLE).to_montgomery();
+
+        b.iter(| | p1.ct_eq(&p2))
+    }
+
+    #[bench]
+    fn montgomery_decompress(b: &mut Bencher) {
+        b.iter(| | BASE_COMPRESSED_MONTGOMERY.decompress());
+    }
+
+    #[bench]
+    fn montgomery_compress(b: &mut Bencher) {
+        let p: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress();
+
+        b.iter(| | p.compress());
+    }
+
+    #[bench]
+    fn montgomery_ladder(b: &mut Bencher) {
+        let mut csprng: OsRng = OsRng::new().unwrap();
+        let s: Scalar = Scalar::random(&mut csprng);
+        let p: MontgomeryPoint = (&Scalar::random(&mut csprng) * &ED25519_BASEPOINT_TABLE).to_montgomery();
+
+        b.iter(| | &s * &p);
     }
 }