From 6747133519aed35239394d96bed8f7b44623f4ac Mon Sep 17 00:00:00 2001 From: khyperia Date: Mon, 21 Aug 2017 21:35:26 -0700 Subject: [PATCH 01/23] Optimize scalar inversion by implementing square() This shows an 11% speedup for invert() --- src/scalar.rs | 51 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 50 insertions(+), 1 deletion(-) diff --git a/src/scalar.rs b/src/scalar.rs index 876061d..b306aee 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -568,7 +568,7 @@ impl UnpackedScalar { let mut y = UnpackedScalar::one(); // Run through bits of l-2 from highest to least for bit in constants::l_minus_2.bits().iter().rev() { - y = UnpackedScalar::multiply_add(&y, &y, &UnpackedScalar::zero()); + y = y.square(); if *bit == 1 { y = UnpackedScalar::multiply_add(&y, self, &UnpackedScalar::zero()); } @@ -576,6 +576,40 @@ impl UnpackedScalar { y } + /// Compute `a^2 (mod l)`. + pub fn square(&self) -> UnpackedScalar { + let a = self.0; + let mut result = [0i64; 24]; + + result[0] = a[0]*a[0]; + result[1] = 2i64 * a[0]*a[1]; + result[2] = 2i64 * (a[0]*a[2]) + a[1]*a[1]; + result[3] = 2i64 * (a[0]*a[3] + a[1]*a[2]); + result[4] = 2i64 * (a[0]*a[4] + a[1]*a[3]) + a[2]*a[2]; + result[5] = 2i64 * (a[0]*a[5] + a[1]*a[4] + a[2]*a[3]); + result[6] = 2i64 * (a[0]*a[6] + a[1]*a[5] + a[2]*a[4]) + a[3]*a[3]; + result[7] = 2i64 * (a[0]*a[7] + a[1]*a[6] + a[2]*a[5] + a[3]*a[4]); + result[8] = 2i64 * (a[0]*a[8] + a[1]*a[7] + a[2]*a[6] + a[3]*a[5]) + a[4]*a[4]; + result[9] = 2i64 * (a[0]*a[9] + a[1]*a[8] + a[2]*a[7] + a[3]*a[6] + a[4]*a[5]); + result[10] = 2i64 * (a[0]*a[10] + a[1]*a[9] + a[2]*a[8] + a[3]*a[7] + a[4]*a[6]) + a[5]*a[5]; + result[11] = 2i64 * (a[0]*a[11] + a[1]*a[10] + a[2]*a[9] + a[3]*a[8] + a[4]*a[7] + a[5]*a[6]); + result[12] = 2i64 * (a[1]*a[11] + a[2]*a[10] + a[3]*a[9] + a[4]*a[8] + a[5]*a[7]) + a[6]*a[6]; + result[13] = 2i64 * (a[2]*a[11] + a[3]*a[10] + a[4]*a[9] + a[5]*a[8] + a[6]*a[7]); + result[14] = 2i64 * (a[3]*a[11] + a[4]*a[10] + a[5]*a[9] + a[6]*a[8]) + a[7]*a[7]; + result[15] = 2i64 * (a[4]*a[11] + a[5]*a[10] + a[6]*a[9] + a[7]*a[8]); + result[16] = 2i64 * (a[5]*a[11] + a[6]*a[10] + a[7]*a[9]) + a[8]*a[8]; + result[17] = 2i64 * (a[6]*a[11] + a[7]*a[10] + a[8]*a[9]); + result[18] = 2i64 * (a[7]*a[11] + a[8]*a[10]) + a[9]*a[9]; + result[19] = 2i64 * (a[8]*a[11] + a[9]*a[10]); + result[20] = 2i64 * (a[9]*a[11]) + a[10]*a[10]; + result[21] = 2i64 * (a[10]*a[11]); + result[22] = a[11]*a[11]; + result[23] = 0i64; + + // Reduce limbs + UnpackedScalar::reduce_limbs(&mut result) + } + /// Compute `ab+c (mod l)`. pub fn multiply_add(a: &UnpackedScalar, b: &UnpackedScalar, @@ -838,6 +872,15 @@ mod test { } } + #[test] + fn square() { + let expected = Scalar::multiply_add(&X, &X, &Scalar::zero()); + let actual = X.unpack().square().pack(); + for i in 0..32 { + assert!(expected[i] == actual[i]); + } + } + #[test] fn scalar_reduce() { let mut bignum = [0u8; 64]; @@ -914,6 +957,12 @@ mod bench { b.iter(|| x.invert()); } + #[bench] + fn square(b: &mut Bencher) { + let x = X.unpack(); + b.iter(|| x.square()); + } + #[bench] fn scalar_unpacked_multiply_add(b: &mut Bencher) { let x = X.unpack(); From 91a7c641c26c79a653130f0a319439145026c1b6 Mon Sep 17 00:00:00 2001 From: Brian Smith Date: Sun, 3 Sep 2017 16:49:26 -1000 Subject: [PATCH 02/23] Use more efficient addition chain for scalar inversion. Use the addition chain from https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion. In my benchmarking, this consistently runs at least 20% faster. --- src/scalar.rs | 59 +++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 53 insertions(+), 6 deletions(-) diff --git a/src/scalar.rs b/src/scalar.rs index b306aee..5490a21 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -2,11 +2,13 @@ // // This file is part of curve25519-dalek. // Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence +// Portions Copyright 2017 Brian Smith // See LICENSE for licensing information. // // Authors: // - Isis Agora Lovecruft // - Henry de Valence +// - Brian Smith //! Arithmetic for scalar multiplication. //! @@ -565,14 +567,59 @@ impl UnpackedScalar { /// Compute the multiplicative inverse of this scalar. pub fn invert(&self) -> UnpackedScalar { - let mut y = UnpackedScalar::one(); - // Run through bits of l-2 from highest to least - for bit in constants::l_minus_2.bits().iter().rev() { - y = y.square(); - if *bit == 1 { - y = UnpackedScalar::multiply_add(&y, self, &UnpackedScalar::zero()); + // This is a direct transliteration of the addition chain from + // https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion + // as it was published on 2017-09-03. + + let _1 = *self; + let _10 = _1.square(); + let _11 = UnpackedScalar::multiply_add(&_10, &_1, &UnpackedScalar::zero()); + let _101 = UnpackedScalar::multiply_add(&_11, &_10, &UnpackedScalar::zero()); + let _111 = UnpackedScalar::multiply_add(&_101, &_10, &UnpackedScalar::zero()); + let _1001 = UnpackedScalar::multiply_add(&_111, &_10, &UnpackedScalar::zero()); + let _1011 = UnpackedScalar::multiply_add(&_1001, &_10, &UnpackedScalar::zero()); + let _1101 = UnpackedScalar::multiply_add(&_1011, &_10, &UnpackedScalar::zero()); + let _1111 = UnpackedScalar::multiply_add(&_1101, &_10, &UnpackedScalar::zero()); + + // 0b10000 + let mut y = UnpackedScalar::multiply_add(&_1111, &_1, &UnpackedScalar::zero()); + + #[inline] + fn square_multiply(y: &mut UnpackedScalar, squarings: usize, x: &UnpackedScalar) { + for _ in 0..squarings { + *y = y.square(); } + *y = UnpackedScalar::multiply_add(y, x, &UnpackedScalar::zero()); } + + square_multiply(&mut y, 123 + 3, &_101); + square_multiply(&mut y, 2 + 2, &_11); + square_multiply(&mut y, 1 + 4, &_1111); + square_multiply(&mut y, 1 + 4, &_1111); + square_multiply(&mut y, 4, &_1001); + square_multiply(&mut y, 4, &_1101); + square_multiply(&mut y, 3, &_111); + square_multiply(&mut y, 1 + 3, &_101); + square_multiply(&mut y, 3 + 3, &_101); + square_multiply(&mut y, 3, &_111); + square_multiply(&mut y, 1 + 4, &_1111); + square_multiply(&mut y, 2 + 3, &_111); + square_multiply(&mut y, 2 + 2, &_11); + square_multiply(&mut y, 1 + 4, &_1011); + square_multiply(&mut y, 2 + 4, &_1011); + square_multiply(&mut y, 6 + 4, &_1001); + square_multiply(&mut y, 2 + 2, &_11); + square_multiply(&mut y, 3 + 2, &_11); + square_multiply(&mut y, 3 + 2, &_11); + square_multiply(&mut y, 1 + 4, &_1001); + square_multiply(&mut y, 1 + 3, &_111); + square_multiply(&mut y, 2 + 4, &_1111); + square_multiply(&mut y, 1 + 4, &_1011); + square_multiply(&mut y, 3, &_101); + square_multiply(&mut y, 2 + 4, &_1111); + square_multiply(&mut y, 3, &_101); + square_multiply(&mut y, 1 + 2, &_11); + y } From 028140bb334a1d49b0acee05fea8fe4f42bc210f Mon Sep 17 00:00:00 2001 From: Brian Smith Date: Mon, 4 Sep 2017 09:19:58 -1000 Subject: [PATCH 03/23] Reformat addition chain window building code to better show pattern. Make the 2 digit, `_10`, the first argument to more closely match the Haskell code in the source article. Align the code into columns to further clarify the patterns. --- src/scalar.rs | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/scalar.rs b/src/scalar.rs index 5490a21..8cf4c95 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -571,17 +571,17 @@ impl UnpackedScalar { // https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion // as it was published on 2017-09-03. - let _1 = *self; - let _10 = _1.square(); - let _11 = UnpackedScalar::multiply_add(&_10, &_1, &UnpackedScalar::zero()); - let _101 = UnpackedScalar::multiply_add(&_11, &_10, &UnpackedScalar::zero()); - let _111 = UnpackedScalar::multiply_add(&_101, &_10, &UnpackedScalar::zero()); - let _1001 = UnpackedScalar::multiply_add(&_111, &_10, &UnpackedScalar::zero()); - let _1011 = UnpackedScalar::multiply_add(&_1001, &_10, &UnpackedScalar::zero()); - let _1101 = UnpackedScalar::multiply_add(&_1011, &_10, &UnpackedScalar::zero()); - let _1111 = UnpackedScalar::multiply_add(&_1101, &_10, &UnpackedScalar::zero()); + let _1 = *self; + let _10 = _1.square(); + let _11 = UnpackedScalar::multiply_add(&_10, &_1, &UnpackedScalar::zero()); + let _101 = UnpackedScalar::multiply_add(&_10, &_11, &UnpackedScalar::zero()); + let _111 = UnpackedScalar::multiply_add(&_10, &_101, &UnpackedScalar::zero()); + let _1001 = UnpackedScalar::multiply_add(&_10, &_111, &UnpackedScalar::zero()); + let _1011 = UnpackedScalar::multiply_add(&_10, &_1001, &UnpackedScalar::zero()); + let _1101 = UnpackedScalar::multiply_add(&_10, &_1011, &UnpackedScalar::zero()); + let _1111 = UnpackedScalar::multiply_add(&_10, &_1101, &UnpackedScalar::zero()); - // 0b10000 + // _10000 let mut y = UnpackedScalar::multiply_add(&_1111, &_1, &UnpackedScalar::zero()); #[inline] From 7ed9eb8617c5de3b7af533f19c871bc7d6e4c930 Mon Sep 17 00:00:00 2001 From: Brian Smith Date: Mon, 4 Sep 2017 09:34:04 -1000 Subject: [PATCH 04/23] Replace one multiplication with a squaring in scalar inversion. This brings the code up to date with the 2017-09-04 version of the source article. --- src/scalar.rs | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/scalar.rs b/src/scalar.rs index 8cf4c95..00124c1 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -573,13 +573,13 @@ impl UnpackedScalar { let _1 = *self; let _10 = _1.square(); - let _11 = UnpackedScalar::multiply_add(&_10, &_1, &UnpackedScalar::zero()); - let _101 = UnpackedScalar::multiply_add(&_10, &_11, &UnpackedScalar::zero()); - let _111 = UnpackedScalar::multiply_add(&_10, &_101, &UnpackedScalar::zero()); - let _1001 = UnpackedScalar::multiply_add(&_10, &_111, &UnpackedScalar::zero()); - let _1011 = UnpackedScalar::multiply_add(&_10, &_1001, &UnpackedScalar::zero()); - let _1101 = UnpackedScalar::multiply_add(&_10, &_1011, &UnpackedScalar::zero()); - let _1111 = UnpackedScalar::multiply_add(&_10, &_1101, &UnpackedScalar::zero()); + let _100 = _10.square(); + let _11 = UnpackedScalar::multiply_add(&_10, &_1, &UnpackedScalar::zero()); + let _101 = UnpackedScalar::multiply_add(&_10, &_11, &UnpackedScalar::zero()); + let _111 = UnpackedScalar::multiply_add(&_10, &_101, &UnpackedScalar::zero()); + let _1001 = UnpackedScalar::multiply_add(&_10, &_111, &UnpackedScalar::zero()); + let _1011 = UnpackedScalar::multiply_add(&_10, &_1001, &UnpackedScalar::zero()); + let _1111 = UnpackedScalar::multiply_add(&_100, &_1011, &UnpackedScalar::zero()); // _10000 let mut y = UnpackedScalar::multiply_add(&_1111, &_1, &UnpackedScalar::zero()); @@ -597,8 +597,8 @@ impl UnpackedScalar { square_multiply(&mut y, 1 + 4, &_1111); square_multiply(&mut y, 1 + 4, &_1111); square_multiply(&mut y, 4, &_1001); - square_multiply(&mut y, 4, &_1101); - square_multiply(&mut y, 3, &_111); + square_multiply(&mut y, 2, &_11); + square_multiply(&mut y, 1 + 4, &_1111); square_multiply(&mut y, 1 + 3, &_101); square_multiply(&mut y, 3 + 3, &_101); square_multiply(&mut y, 3, &_111); From acd3826fe28d2dddcc97e6ead9bb1acb834e1c1c Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 7 Sep 2017 20:31:06 +0000 Subject: [PATCH 05/23] Implement Montgomery arithmetic and laddering. * ADDs part of https://github.com/isislovecruft/curve25519-dalek/issues/47 --- src/constants.rs | 9 ++ src/constants_32bit.rs | 3 + src/constants_64bit.rs | 3 + src/edwards.rs | 45 +++++- src/field.rs | 8 +- src/montgomery.rs | 334 +++++++++++++++++++++++++++++++++++++++-- 6 files changed, 385 insertions(+), 17 deletions(-) diff --git a/src/constants.rs b/src/constants.rs index c8d6b8d..4c08102 100644 --- a/src/constants.rs +++ b/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 = diff --git a/src/constants_32bit.rs b/src/constants_32bit.rs index 88db84d..02b629b 100644 --- a/src/constants_32bit.rs +++ b/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) diff --git a/src/constants_64bit.rs b/src/constants_64bit.rs index 046da20..8beee81 100644 --- a/src/constants_64bit.rs +++ b/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) diff --git a/src/edwards.rs b/src/edwards.rs index 2456954..9a03295 100644 --- a/src/edwards.rs +++ b/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 { + fn convert_to_montgomery(&self) -> Option { // 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 { + let u: Option = 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 { + let u: Option = 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 { + self.to_projective().to_montgomery() + } + /// Compress this point to `CompressedEdwardsY` format. pub fn compress_edwards(&self) -> CompressedEdwardsY { self.to_projective().compress_edwards() diff --git a/src/field.rs b/src/field.rs index 617adbc..730c3fa 100644 --- a/src/field.rs +++ b/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 diff --git a/src/montgomery.rs b/src/montgomery.rs index 29e1b20..0b2c16a 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -8,7 +8,19 @@ // - Isis Agora Lovecruft // - Henry de Valence -//! 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 { + pub fn decompress_edwards(&self) -> Option { 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()); + } } From 2939d26b5c42b414b7626f504868948dd1e69171 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 14 Sep 2017 01:54:28 +0000 Subject: [PATCH 06/23] Remove direct compression methods between points in curve models. compress_edwards() is now named compress() and works only on points which are in Edwards form. Similarly, compress_montgomery() is now also called compress(), and it only works on point already in Mongomery form. To switch between forms, use to_montgomery(). Conversion from Montgomery directly to Edwards is not yet implemented. * CHANGE the API requested in https://github.com/isislovecruft/curve25519-dalek/issues/47, hopefully for the better. --- src/decaf.rs | 13 +++- src/edwards.rs | 182 ++++++++++++++++++++++------------------------ src/montgomery.rs | 47 ++++++------ 3 files changed, 121 insertions(+), 121 deletions(-) diff --git a/src/decaf.rs b/src/decaf.rs index f8eb0c4..23511dd 100644 --- a/src/decaf.rs +++ b/src/decaf.rs @@ -194,6 +194,11 @@ impl<'de> Deserialize<'de> for DecafPoint { pub struct DecafPoint(pub ExtendedPoint); impl DecafPoint { + /// Convert this `DecafPoint` to its underlying `ExtendedPoint`. + pub fn to_edwards(&self) -> ExtendedPoint { + self.0 + } + /// Compress in Decaf format. pub fn compress(&self) -> CompressedDecaf { // Q: Do we want to encode twisted or untwisted? @@ -752,7 +757,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.to_edwards().compress(), CompressedEdwardsY::identity()); } #[test] @@ -764,11 +769,11 @@ mod test { #[test] fn decaf_basepoint_roundtrip() { let bp_compressed_decaf = constants::DECAF_ED25519_BASEPOINT_POINT.compress(); - let bp_recaf = bp_compressed_decaf.decompress().unwrap().0; + let bp_recaf = bp_compressed_decaf.decompress().unwrap().to_edwards(); // 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] @@ -838,7 +843,7 @@ mod test { for _ in 0..100 { let P = DecafPoint::random(&mut rng); // Check that P is on the curve - assert!(P.0.is_valid()); + assert!(P.to_edwards().is_valid()); // Check that P is in the image of the decaf map P.compress(); } diff --git a/src/edwards.rs b/src/edwards.rs index 9a03295..0b8b333 100644 --- a/src/edwards.rs +++ b/src/edwards.rs @@ -89,7 +89,6 @@ use core::ops::Index; use constants; use field::FieldElement; use scalar::Scalar; -use montgomery::CompressedMontgomeryU; use montgomery::MontgomeryPoint; use subtle::slices_equal; @@ -124,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 } @@ -170,7 +168,7 @@ impl Serialize for ExtendedPoint { fn serialize(&self, serializer: S) -> Result where S: Serializer { - serializer.serialize_bytes(self.compress_edwards().as_bytes()) + serializer.serialize_bytes(self.compress().as_bytes()) } } @@ -394,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(), - other.compress_edwards().as_bytes()) + slices_equal(self.compress().as_bytes(), + other.compress().as_bytes()) } } @@ -437,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; @@ -448,63 +446,67 @@ impl ProjectivePoint { CompressedEdwardsY(s) } - /// Convert this point to a Montgomery u-coordinate (affine). - /// Note that this discards the sign. + /// Convert this projective point in the Edwards model to its equivalent + /// projective point on the Montgomery form of the curve. /// - /// # Return - /// - `None` if `self` is the identity point; - /// - `Some(FieldElement)` otherwise. + /// Taking the Montgomery curve equation in affine coordinates: /// - fn convert_to_montgomery(&self) -> Option { - // u = (1 + y) / (1 - y) - // v = sqrt(-486664) * u / x - // - // since y = Y/Z, x = X/Z, - // - // u = (1 + Y/Z) / (1 - Y/Z); - // = (Z + Y) / (Z - Y); - // - // exceptional points: - // y = 1 <=> Y/Z = 1 <=> Z - Y = 0 - let Z_plus_Y = &self.Z + &self.Y; - let Z_minus_Y = &self.Z - &self.Y; - let u = &Z_plus_Y * &Z_minus_Y.invert(); - - if Z_minus_Y.is_zero() == 0u8 { - Some(u) - } else { - None - } - } - - /// Convert this point to a `CompressedMontgomeryU`. - /// Note that this discards the sign. + ///     E_(A,B) = Bv² = u³ + Au² + u   (1) /// - /// # Return - /// - `None` if `self` is the identity point; - /// - `Some(CompressedMontgomeryU)` otherwise. + /// and given its relations to the coordinates of the Edwards model: /// - pub fn compress_montgomery(&self) -> Option { - let u: Option = 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. + ///     u = (1+y)/(1-y)        (2) + ///     v = (λu)/(x) /// - /// DOCDOC - pub fn to_montgomery(&self) -> Option { - let u: Option = self.convert_to_montgomery(); - - if u.is_some() { - Some(MontgomeryPoint{ U: u.unwrap(), Z: FieldElement::one() }) - } else { - None + /// Converting from affine to projective coordinates in the Montgomery + /// model, we arrive at: + /// + ///     u = (Z+Y)/(Z-Y)        (3) + ///     v = λ * ((Z+Y)/(Z-Y)) * (Z/X) + /// + /// The transition between affine and projective is given by + /// + ///     u → U/W        (4) + ///     v → V/W + /// + /// thus the Montgomery curve equation (1) becomes + /// + ///     E_(A,B) : BV²W = U³ + AU²W + UW² ⊆ 𝗣^2  (5) + /// + /// 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) (6) + /// + /// 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² (7) + /// + /// 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, } } } @@ -547,25 +549,15 @@ impl ExtendedPoint { } } - /// DOCDOC - pub fn to_montgomery(&self) -> Option { + /// Convert this point to its equivalent on the Montgomery form of the + /// curve. + pub fn to_montgomery(&self) -> MontgomeryPoint { self.to_projective().to_montgomery() } /// Compress this point to `CompressedEdwardsY` format. - pub fn compress_edwards(&self) -> CompressedEdwardsY { - self.to_projective().compress_edwards() - } - - /// 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 { - self.to_projective().compress_montgomery() + pub fn compress(&self) -> CompressedEdwardsY { + self.to_projective().compress() } } @@ -1342,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 @@ -1365,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); } @@ -1373,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 for ExtendedPoint` @@ -1382,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 for ExtendedPoint` @@ -1391,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 for ExtendedPoint` @@ -1401,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 @@ -1426,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(), - also_aB.compress_edwards()); + assert_eq!( aB.compress(), + 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 @@ -1452,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); } @@ -1474,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); } @@ -1489,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. @@ -1517,7 +1509,7 @@ mod test { #[test] fn compressed_identity() { - assert_eq!(ExtendedPoint::identity().compress_edwards(), + assert_eq!(ExtendedPoint::identity().compress(), CompressedEdwardsY::identity()); } @@ -1555,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] @@ -1579,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] @@ -1589,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] @@ -1604,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()); } } @@ -1616,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] @@ -1652,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] diff --git a/src/montgomery.rs b/src/montgomery.rs index 0b2c16a..f2e0268 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -119,7 +119,7 @@ impl CompressedMontgomeryU { /// # Returns /// /// A projective `MontgomeryPoint` corresponding to this compressed point. - pub fn decompress_montgomery(&self) -> MontgomeryPoint { + pub fn decompress(&self) -> MontgomeryPoint { MontgomeryPoint{ // XXX is it a problem here if we're not using a canonical encoding? —isis U: FieldElement::from_bytes(&self.0), @@ -257,8 +257,8 @@ impl Identity for MontgomeryPoint { /// `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()) + slices_equal(self.compress().as_bytes(), + that.compress().as_bytes()) } } @@ -301,7 +301,7 @@ impl MontgomeryPoint { /// # Returns /// /// A `CompressedMontgomeryU`. - pub fn compress_montgomery(&self) -> CompressedMontgomeryU { + pub fn compress(&self) -> CompressedMontgomeryU { let u_affine: FieldElement = &self.U * &self.W.invert(); CompressedMontgomeryU(u_affine.to_bytes()) @@ -425,15 +425,15 @@ mod test { /// 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_COMPRESSED_MONTGOMERY); } /// Test Montgomery conversion against the X25519 basepoint. #[test] fn basepoint_from_montgomery() { - assert_eq!(BASE_COMPRESSED_MONTGOMERY.decompress_edwards().unwrap().compress_edwards(), - constants::BASE_CMPRSSD); + assert_eq!(BASE_COMPRESSED_MONTGOMERY, + constants::BASE_CMPRSSD.decompress().unwrap().to_montgomery().compress()); } /// If u = -1, then v^2 = u*(u^2+486662*u+1) = 486660. @@ -448,26 +448,28 @@ mod test { assert!(div_by_zero_u.decompress_edwards().is_none()); } - /// Montgomery compression of the identity point should - /// fail (it's sent to infinity). + /// Montgomery compression of the identity point should not fail (since the + /// 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() { - let p = BASE_COMPRESSED_MONTGOMERY.decompress_montgomery(); + 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_montgomery().differential_double(); + let q: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress().differential_double(); - assert_eq!(p.compress_montgomery().unwrap(), q.compress_montgomery()); + assert_eq!(p.to_montgomery().compress(), q.compress()); } #[test] @@ -480,13 +482,13 @@ mod test { 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 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_montgomery(), (&p1 + &p2).compress_montgomery().unwrap()); + assert_eq!(result.compress(), (&p1 + &p2).to_montgomery().compress()); } #[test] @@ -495,20 +497,21 @@ mod test { 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 p_montgomery: MontgomeryPoint = p_edwards.to_montgomery(); let expected = &s * &p_edwards; let result = &s * &p_montgomery; - assert_eq!(result.compress_montgomery(), expected.compress_montgomery().unwrap()) + 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_montgomery() * &two; - let mut expected: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double(); + let result: MontgomeryPoint = &BASE_COMPRESSED_MONTGOMERY.decompress() * &two; + let expected: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double(); + + assert_eq!(result.compress(), expected.to_montgomery().compress()); - assert_eq!(result.compress_montgomery(), expected.compress_montgomery().unwrap()); } } From 6be10341e2003648344661a72b28a39157ab41e7 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 14 Sep 2017 02:08:11 +0000 Subject: [PATCH 07/23] Add benchmarks for Mongomery point (de)compression and laddering. --- src/montgomery.rs | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/src/montgomery.rs b/src/montgomery.rs index f2e0268..1f4f52d 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -512,6 +512,35 @@ mod test { let expected: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double(); 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_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); } } From d86fcd4e17e8a6cdb3913d3a2fa82a91d8403801 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 14 Sep 2017 02:40:18 +0000 Subject: [PATCH 08/23] Use subtle branch with ConditionallySwappable trait for now. --- Cargo.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Cargo.toml b/Cargo.toml index 8c22506..c1cd666 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -35,6 +35,8 @@ version = "0.6" [dependencies.subtle] version = "^0.2" default-features = false +git = "https://github.com/isislovecruft/subtle" +branch = "feature/conditional-swap_r1" [dependencies.generic-array] # same version that digest depends on From 7e53499a1046a0c4fb79c0770ead17ec40c42109 Mon Sep 17 00:00:00 2001 From: Andrew Moon Date: Sun, 24 Sep 2017 22:24:35 -0500 Subject: [PATCH 09/23] optimized scalar implementations for 32/64 bit --- src/constants_32bit.rs | 21 ++ src/constants_64bit.rs | 13 + src/lib.rs | 7 + src/scalar.rs | 429 +++++-------------------------- src/scalar_32bit.rs | 559 +++++++++++++++++++++++++++++++++++++++++ src/scalar_64bit.rs | 475 ++++++++++++++++++++++++++++++++++ 6 files changed, 1143 insertions(+), 361 deletions(-) create mode 100644 src/scalar_32bit.rs create mode 100644 src/scalar_64bit.rs diff --git a/src/constants_32bit.rs b/src/constants_32bit.rs index 88db84d..dc7468b 100644 --- a/src/constants_32bit.rs +++ b/src/constants_32bit.rs @@ -19,6 +19,7 @@ #![allow(non_snake_case)] use field_32bit::FieldElement32; +use scalar_32bit::Scalar32; use edwards::ExtendedPoint; use edwards::AffineNielsPoint; use edwards::EdwardsBasepointTable; @@ -93,6 +94,26 @@ pub const SQRT_MINUS_HALF: FieldElement32 = FieldElement32([ // sqrtMinusHalf -17256545, 3971863, 28865457, -1750208, 27359696, -16640980, 12573105, 1002827, -163343, 11073975, ]); +/// `L` is the order of base point, i.e. 2^252 + +/// 27742317777372353535851937790883648493 +pub const L: Scalar32 = Scalar32([ 0x1cf5d3ed, 0x009318d2, 0x1de73596, 0x1df3bd45, + 0x0000014d, 0x00000000, 0x00000000, 0x00000000, + 0x00100000 ]); + +/// `L` * `LFACTOR` = -1 (mod 2^29) +pub const LFACTOR: u32 = 0x12547e1b; + +/// `R` = R % L where R = 2^261 +pub const R: Scalar32 = Scalar32([ 0x114df9ed, 0x1a617303, 0x0f7c098c, 0x16793167, + 0x1ffd656e, 0x1fffffff, 0x1fffffff, 0x1fffffff, + 0x000fffff ]); + +/// `RR` = (R^2) % L where R = 2^261 +pub const RR: Scalar32 = Scalar32([ 0x0b5f9d12, 0x1e141b17, 0x158d7f3d, 0x143f3757, + 0x1972d781, 0x042feb7c, 0x1ceec73d, 0x1e184d1e, + 0x0005046d ]); + + /// Basepoint has y = 4/5. This is called `_POINT` to distinguish it from `_TABLE`, which should /// be used for scalar multiplication (it's much faster). pub const ED25519_BASEPOINT_POINT: ExtendedPoint = ExtendedPoint{ diff --git a/src/constants_64bit.rs b/src/constants_64bit.rs index 046da20..d058bf0 100644 --- a/src/constants_64bit.rs +++ b/src/constants_64bit.rs @@ -19,6 +19,7 @@ #![allow(non_snake_case)] use field_64bit::FieldElement64; +use scalar_64bit::Scalar64; use edwards::ExtendedPoint; use edwards::AffineNielsPoint; use edwards::EdwardsBasepointTable; @@ -66,6 +67,18 @@ pub const SQRT_MINUS_APLUS2: FieldElement64 = FieldElement64([1693982333959686, /// `SQRT_MINUS_HALF` is sqrt(-1/2) pub const SQRT_MINUS_HALF: FieldElement64 = FieldElement64([266547196637087, 2134345371906993, 1135042577398223, 67298593331632, 743161882051057]); +/// `L` is the order of base point, i.e. 2^252 + 27742317777372353535851937790883648493 +pub const L: Scalar64 = Scalar64([ 0x0002631a5cf5d3ed, 0x000dea2f79cd6581, 0x000000000014def9, 0x0000000000000000, 0x0000100000000000 ]); + +/// `L` * `LFACTOR` = -1 (mod 2^51) +pub const LFACTOR: u64 = 0x51da312547e1b; + +/// `R` = R % L where R = 2^260 +pub const R: Scalar64 = Scalar64([ 0x000f48bd6721e6ed, 0x0003bab5ac67e45a, 0x000fffffeb35e51b, 0x000fffffffffffff, 0x00000fffffffffff ]); + +/// `RR` = (R^2) % L where R = 2^260 +pub const RR: Scalar64 = Scalar64([ 0x0009d265e952d13b, 0x000d63c715bea69f, 0x0005be65cb687604, 0x0003dceec73d217f, 0x000009411b7c309a ]); + /// Basepoint has y = 4/5. This is called `_POINT` to distinguish it from `_TABLE`, which should /// be used for scalar multiplication (it's much faster). pub const ED25519_BASEPOINT_POINT: ExtendedPoint = ExtendedPoint{ diff --git a/src/lib.rs b/src/lib.rs index 31c18ff..be018f2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -41,6 +41,8 @@ extern crate test; #[cfg(test)] extern crate sha2; +// this appears to only be used for serde support right now? +#[cfg(feature = "serde")] #[macro_use] extern crate arrayref; @@ -71,6 +73,11 @@ mod field_32bit; mod field_64bit; pub mod scalar; +#[cfg(not(feature="radix_51"))] +mod scalar_32bit; +#[cfg(feature="radix_51")] +mod scalar_64bit; + pub mod edwards; pub mod montgomery; diff --git a/src/scalar.rs b/src/scalar.rs index 00124c1..399d00d 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -22,12 +22,11 @@ //! //! The `Scalar` struct represents an element in ℤ/lℤ. //! -//! Arithmetic operations on `Scalar`s are done using 12 21-bit limbs. -//! However, in contrast to `FieldElement`s, `Scalar`s are stored in +//! In contrast to `FieldElement`s, `Scalar`s are stored in //! memory as bytes, allowing easy access to the bits of the `Scalar` //! when multiplying a point by a scalar. For efficient arithmetic -//! between two scalars, the `UnpackedScalar` struct is stored as -//! limbs. +//! between two scalars, the `UnpackedScalar` struct (internally +//! either `Scalar32` or `Scalar64`) is stored as limbs. use core::fmt::Debug; use core::ops::Neg; @@ -43,9 +42,6 @@ use rand::Rng; use digest::Digest; use generic_array::typenum::U64; -use constants; -use utils::{load3, load4}; - use subtle::slices_equal; use subtle::ConditionallyAssignable; use subtle::Equal; @@ -108,50 +104,47 @@ impl IndexMut for Scalar { impl<'b> MulAssign<&'b Scalar> for Scalar { fn mul_assign(&mut self, _rhs: &'b Scalar) { - let result = (self as &Scalar) * _rhs; - self.0 = result.0; + *self = Scalar::mul(self, _rhs) } } impl<'a, 'b> Mul<&'b Scalar> for &'a Scalar { type Output = Scalar; fn mul(self, _rhs: &'b Scalar) -> Scalar { - Scalar::multiply_add(self, _rhs, &Scalar::zero()) + Scalar::mul(self, _rhs) } } impl<'b> AddAssign<&'b Scalar> for Scalar { fn add_assign(&mut self, _rhs: &'b Scalar) { - *self = Scalar::multiply_add(&Scalar::one(), self, _rhs); + *self = Scalar::add(self, _rhs); } } impl<'a, 'b> Add<&'b Scalar> for &'a Scalar { type Output = Scalar; fn add(self, _rhs: &'b Scalar) -> Scalar { - Scalar::multiply_add(&Scalar::one(), self, _rhs) + Scalar::add(self, _rhs) } } impl<'b> SubAssign<&'b Scalar> for Scalar { fn sub_assign(&mut self, _rhs: &'b Scalar) { - // (l-1)*_rhs + self = self - _rhs - *self = Scalar::multiply_add(&constants::l_minus_1, _rhs, self); + *self = Scalar::sub(self, _rhs); } } impl<'a, 'b> Sub<&'b Scalar> for &'a Scalar { type Output = Scalar; fn sub(self, _rhs: &'b Scalar) -> Scalar { - // (l-1)*_rhs + self = self - _rhs - Scalar::multiply_add(&constants::l_minus_1, _rhs, self) + Scalar::sub(self, _rhs) } } impl<'a> Neg for &'a Scalar { type Output = Scalar; fn neg(self) -> Scalar { - self * &constants::l_minus_1 + Scalar::sub(&Scalar::zero(), self) } } @@ -234,6 +227,18 @@ impl<'de> Deserialize<'de> for Scalar { } } +/// An `UnpackedScalar` represents an element of the field GF(l), optimized for speed. +#[cfg(feature="radix_51")] +type UnpackedScalar = Scalar64; +#[cfg(feature="radix_51")] +use scalar_64bit::*; + +/// An `UnpackedScalar` represents an element of the field GF(l), optimized for speed. +#[cfg(not(feature="radix_51"))] +type UnpackedScalar = Scalar32; +#[cfg(not(feature="radix_51"))] +use scalar_32bit::*; + impl Scalar { /// Return a `Scalar` chosen uniformly at random using a user-provided RNG. /// @@ -384,26 +389,6 @@ impl Scalar { naf } - // Unpack a scalar into 12 21-bit limbs. - fn unpack(&self) -> UnpackedScalar { - let mask_21bits: i64 = (1 << 21) - 1; - let mut a = UnpackedScalar([0i64; 12]); - a[ 0] = mask_21bits & load3(&self.0[ 0..]) ; - a[ 1] = mask_21bits & (load4(&self.0[ 2..]) >> 5); - a[ 2] = mask_21bits & (load3(&self.0[ 5..]) >> 2); - a[ 3] = mask_21bits & (load4(&self.0[ 7..]) >> 7); - a[ 4] = mask_21bits & (load4(&self.0[10..]) >> 4); - a[ 5] = mask_21bits & (load3(&self.0[13..]) >> 1); - a[ 6] = mask_21bits & (load4(&self.0[15..]) >> 6); - a[ 7] = mask_21bits & (load3(&self.0[18..]) >> 3); - a[ 8] = mask_21bits & load3(&self.0[21..]) ; - a[ 9] = mask_21bits & (load4(&self.0[23..]) >> 5); - a[10] = mask_21bits & (load3(&self.0[26..]) >> 2); - a[11] = load4(&self.0[28..]) >> 7 ; - - a - } - /// Write this scalar in radix 16, with coefficients in `[-8,8)`, /// i.e., compute `a_i` such that /// @@ -442,127 +427,41 @@ impl Scalar { output } - /// Compute `ab+c (mod l)`. - /// XXX should this exist, or should we just have Mul, Add etc impls - /// that unpack and then call UnpackedScalar::multiply_add ? - pub fn multiply_add(a: &Scalar, b: &Scalar, c: &Scalar) -> Scalar { - // Unpack scalars into limbs - let al = a.unpack(); - let bl = b.unpack(); - let cl = c.unpack(); + /// Unpack this `Scalar` to an `UnpackedScalar` + pub fn unpack(&self) -> UnpackedScalar { + UnpackedScalar::from_bytes(&self.0) + } - // Multiply and repack - UnpackedScalar::multiply_add(&al, &bl, &cl).pack() + /// Compute `a + b` (mod l) + pub fn add(a: &Scalar, b: &Scalar) -> Scalar { + UnpackedScalar::add(&a.unpack(), &b.unpack()).pack() + } + + /// Compute `a - b` (mod l). + pub fn sub(a: &Scalar, b: &Scalar) -> Scalar { + UnpackedScalar::sub(&a.unpack(), &b.unpack()).pack() + } + + /// Compute `a * b` (mod l). + pub fn mul(a: &Scalar, b: &Scalar) -> Scalar { + UnpackedScalar::mul(&a.unpack(), &b.unpack()).pack() + } + + /// Compute `(a * b) + c` (mod l). + pub fn multiply_add(a: &Scalar, b: &Scalar, c: &Scalar) -> Scalar { + UnpackedScalar::add(&UnpackedScalar::mul(&a.unpack(), &b.unpack()), &c.unpack()).pack() } /// Reduce a 512-bit little endian number mod l pub fn reduce(input: &[u8; 64]) -> Scalar { - let mut s = [0i64; 24]; - - // XXX express this as two unpack_limbs - // some issues re: masking with the top byte of the 32byte input - let mask_21bits: i64 = (1 << 21) -1; - s[0] = mask_21bits & load3(&input[ 0..]) ; - s[1] = mask_21bits & (load4(&input[ 2..]) >> 5); - s[2] = mask_21bits & (load3(&input[ 5..]) >> 2); - s[3] = mask_21bits & (load4(&input[ 7..]) >> 7); - s[4] = mask_21bits & (load4(&input[10..]) >> 4); - s[5] = mask_21bits & (load3(&input[13..]) >> 1); - s[6] = mask_21bits & (load4(&input[15..]) >> 6); - s[7] = mask_21bits & (load3(&input[18..]) >> 3); - s[8] = mask_21bits & load3(&input[21..]) ; - s[9] = mask_21bits & (load4(&input[23..]) >> 5); - s[10] = mask_21bits & (load3(&input[26..]) >> 2); - s[11] = mask_21bits & (load4(&input[28..]) >> 7); - s[12] = mask_21bits & (load4(&input[31..]) >> 4); - s[13] = mask_21bits & (load3(&input[34..]) >> 1); - s[14] = mask_21bits & (load4(&input[36..]) >> 6); - s[15] = mask_21bits & (load3(&input[39..]) >> 3); - s[16] = mask_21bits & load3(&input[42..]) ; - s[17] = mask_21bits & (load4(&input[44..]) >> 5); - s[18] = mask_21bits & (load3(&input[47..]) >> 2); - s[19] = mask_21bits & (load4(&input[49..]) >> 7); - s[20] = mask_21bits & (load4(&input[52..]) >> 4); - s[21] = mask_21bits & (load3(&input[55..]) >> 1); - s[22] = mask_21bits & (load4(&input[57..]) >> 6); - s[23] = load4(&input[60..]) >> 3 ; - - // XXX replacing the previous code in this function with the - // call to reduce_limbs adds two extra carry passes (the ones - // at the top of the reduce_limbs function). Otherwise they - // are identical. The test seems to work OK but it would be - // good to check that this really is OK to add. - UnpackedScalar::reduce_limbs(&mut s).pack() - } -} - -/// The `UnpackedScalar` struct represents an element in ℤ/lℤ as 12 -/// 21-bit limbs. -#[derive(Copy,Clone)] -pub struct UnpackedScalar(pub [i64; 12]); - -impl Index for UnpackedScalar { - type Output = i64; - - fn index(&self, _index: usize) -> &i64 { - &(self.0[_index]) - } -} - -impl IndexMut for UnpackedScalar { - fn index_mut(&mut self, _index: usize) -> &mut i64 { - &mut (self.0[_index]) + UnpackedScalar::from_bytes_wide(input).pack() } } impl UnpackedScalar { /// Pack the limbs of this `UnpackedScalar` into a `Scalar`. fn pack(&self) -> Scalar { - let mut s = Scalar::zero(); - s[0] = (self.0[ 0] >> 0) as u8; - s[1] = (self.0[ 0] >> 8) as u8; - s[2] = ((self.0[ 0] >> 16) | (self.0[ 1] << 5)) as u8; - s[3] = (self.0[ 1] >> 3) as u8; - s[4] = (self.0[ 1] >> 11) as u8; - s[5] = ((self.0[ 1] >> 19) | (self.0[ 2] << 2)) as u8; - s[6] = (self.0[ 2] >> 6) as u8; - s[7] = ((self.0[ 2] >> 14) | (self.0[ 3] << 7)) as u8; - s[8] = (self.0[ 3] >> 1) as u8; - s[9] = (self.0[ 3] >> 9) as u8; - s[10] = ((self.0[ 3] >> 17) | (self.0[ 4] << 4)) as u8; - s[11] = (self.0[ 4] >> 4) as u8; - s[12] = (self.0[ 4] >> 12) as u8; - s[13] = ((self.0[ 4] >> 20) | (self.0[ 5] << 1)) as u8; - s[14] = (self.0[ 5] >> 7) as u8; - s[15] = ((self.0[ 5] >> 15) | (self.0[ 6] << 6)) as u8; - s[16] = (self.0[ 6] >> 2) as u8; - s[17] = (self.0[ 6] >> 10) as u8; - s[18] = ((self.0[ 6] >> 18) | (self.0[ 7] << 3)) as u8; - s[19] = (self.0[ 7] >> 5) as u8; - s[20] = (self.0[ 7] >> 13) as u8; - s[21] = (self.0[ 8] >> 0) as u8; - s[22] = (self.0[ 8] >> 8) as u8; - s[23] = ((self.0[ 8] >> 16) | (self.0[ 9] << 5)) as u8; - s[24] = (self.0[ 9] >> 3) as u8; - s[25] = (self.0[ 9] >> 11) as u8; - s[26] = ((self.0[ 9] >> 19) | (self.0[10] << 2)) as u8; - s[27] = (self.0[10] >> 6) as u8; - s[28] = ((self.0[10] >> 14) | (self.0[11] << 7)) as u8; - s[29] = (self.0[11] >> 1) as u8; - s[30] = (self.0[11] >> 9) as u8; - s[31] = (self.0[11] >> 17) as u8; - - s - } - - /// Return the zero scalar. - pub fn zero() -> UnpackedScalar { - UnpackedScalar([0,0,0,0,0,0,0,0,0,0,0,0]) - } - - /// Return the one scalar. - pub fn one() -> UnpackedScalar { - UnpackedScalar([1,0,0,0,0,0,0,0,0,0,0,0]) + Scalar(self.to_bytes()) } /// Compute the multiplicative inverse of this scalar. @@ -571,25 +470,25 @@ impl UnpackedScalar { // https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion // as it was published on 2017-09-03. - let _1 = *self; - let _10 = _1.square(); - let _100 = _10.square(); - let _11 = UnpackedScalar::multiply_add(&_10, &_1, &UnpackedScalar::zero()); - let _101 = UnpackedScalar::multiply_add(&_10, &_11, &UnpackedScalar::zero()); - let _111 = UnpackedScalar::multiply_add(&_10, &_101, &UnpackedScalar::zero()); - let _1001 = UnpackedScalar::multiply_add(&_10, &_111, &UnpackedScalar::zero()); - let _1011 = UnpackedScalar::multiply_add(&_10, &_1001, &UnpackedScalar::zero()); - let _1111 = UnpackedScalar::multiply_add(&_100, &_1011, &UnpackedScalar::zero()); + let _1 = self.to_montgomery(); + let _10 = _1.montgomery_square(); + let _100 = _10.montgomery_square(); + let _11 = UnpackedScalar::montgomery_mul(&_10, &_1); + let _101 = UnpackedScalar::montgomery_mul(&_10, &_11); + let _111 = UnpackedScalar::montgomery_mul(&_10, &_101); + let _1001 = UnpackedScalar::montgomery_mul(&_10, &_111); + let _1011 = UnpackedScalar::montgomery_mul(&_10, &_1001); + let _1111 = UnpackedScalar::montgomery_mul(&_100, &_1011); // _10000 - let mut y = UnpackedScalar::multiply_add(&_1111, &_1, &UnpackedScalar::zero()); + let mut y = UnpackedScalar::montgomery_mul(&_1111, &_1); #[inline] fn square_multiply(y: &mut UnpackedScalar, squarings: usize, x: &UnpackedScalar) { for _ in 0..squarings { - *y = y.square(); + *y = y.montgomery_square(); } - *y = UnpackedScalar::multiply_add(y, x, &UnpackedScalar::zero()); + *y = UnpackedScalar::montgomery_mul(y, x); } square_multiply(&mut y, 123 + 3, &_101); @@ -620,194 +519,14 @@ impl UnpackedScalar { square_multiply(&mut y, 3, &_101); square_multiply(&mut y, 1 + 2, &_11); - y - } - - /// Compute `a^2 (mod l)`. - pub fn square(&self) -> UnpackedScalar { - let a = self.0; - let mut result = [0i64; 24]; - - result[0] = a[0]*a[0]; - result[1] = 2i64 * a[0]*a[1]; - result[2] = 2i64 * (a[0]*a[2]) + a[1]*a[1]; - result[3] = 2i64 * (a[0]*a[3] + a[1]*a[2]); - result[4] = 2i64 * (a[0]*a[4] + a[1]*a[3]) + a[2]*a[2]; - result[5] = 2i64 * (a[0]*a[5] + a[1]*a[4] + a[2]*a[3]); - result[6] = 2i64 * (a[0]*a[6] + a[1]*a[5] + a[2]*a[4]) + a[3]*a[3]; - result[7] = 2i64 * (a[0]*a[7] + a[1]*a[6] + a[2]*a[5] + a[3]*a[4]); - result[8] = 2i64 * (a[0]*a[8] + a[1]*a[7] + a[2]*a[6] + a[3]*a[5]) + a[4]*a[4]; - result[9] = 2i64 * (a[0]*a[9] + a[1]*a[8] + a[2]*a[7] + a[3]*a[6] + a[4]*a[5]); - result[10] = 2i64 * (a[0]*a[10] + a[1]*a[9] + a[2]*a[8] + a[3]*a[7] + a[4]*a[6]) + a[5]*a[5]; - result[11] = 2i64 * (a[0]*a[11] + a[1]*a[10] + a[2]*a[9] + a[3]*a[8] + a[4]*a[7] + a[5]*a[6]); - result[12] = 2i64 * (a[1]*a[11] + a[2]*a[10] + a[3]*a[9] + a[4]*a[8] + a[5]*a[7]) + a[6]*a[6]; - result[13] = 2i64 * (a[2]*a[11] + a[3]*a[10] + a[4]*a[9] + a[5]*a[8] + a[6]*a[7]); - result[14] = 2i64 * (a[3]*a[11] + a[4]*a[10] + a[5]*a[9] + a[6]*a[8]) + a[7]*a[7]; - result[15] = 2i64 * (a[4]*a[11] + a[5]*a[10] + a[6]*a[9] + a[7]*a[8]); - result[16] = 2i64 * (a[5]*a[11] + a[6]*a[10] + a[7]*a[9]) + a[8]*a[8]; - result[17] = 2i64 * (a[6]*a[11] + a[7]*a[10] + a[8]*a[9]); - result[18] = 2i64 * (a[7]*a[11] + a[8]*a[10]) + a[9]*a[9]; - result[19] = 2i64 * (a[8]*a[11] + a[9]*a[10]); - result[20] = 2i64 * (a[9]*a[11]) + a[10]*a[10]; - result[21] = 2i64 * (a[10]*a[11]); - result[22] = a[11]*a[11]; - result[23] = 0i64; - - // Reduce limbs - UnpackedScalar::reduce_limbs(&mut result) - } - - /// Compute `ab+c (mod l)`. - pub fn multiply_add(a: &UnpackedScalar, - b: &UnpackedScalar, - c: &UnpackedScalar) -> UnpackedScalar { - let mut result = [0i64; 24]; - - // Multiply a and b, and add c - result[0] = c[0] + a[0]*b[0]; - result[1] = c[1] + a[0]*b[1] + a[1]*b[0]; - result[2] = c[2] + a[0]*b[2] + a[1]*b[1] + a[2]*b[0]; - result[3] = c[3] + a[0]*b[3] + a[1]*b[2] + a[2]*b[1] + a[3]*b[0]; - result[4] = c[4] + a[0]*b[4] + a[1]*b[3] + a[2]*b[2] + a[3]*b[1] + a[4]*b[0]; - result[5] = c[5] + a[0]*b[5] + a[1]*b[4] + a[2]*b[3] + a[3]*b[2] + a[4]*b[1] + a[5]*b[0]; - result[6] = c[6] + a[0]*b[6] + a[1]*b[5] + a[2]*b[4] + a[3]*b[3] + a[4]*b[2] + a[5]*b[1] + a[6]*b[0]; - result[7] = c[7] + a[0]*b[7] + a[1]*b[6] + a[2]*b[5] + a[3]*b[4] + a[4]*b[3] + a[5]*b[2] + a[6]*b[1] + a[7]*b[0]; - result[8] = c[8] + a[0]*b[8] + a[1]*b[7] + a[2]*b[6] + a[3]*b[5] + a[4]*b[4] + a[5]*b[3] + a[6]*b[2] + a[7]*b[1] + a[8]*b[0]; - result[9] = c[9] + a[0]*b[9] + a[1]*b[8] + a[2]*b[7] + a[3]*b[6] + a[4]*b[5] + a[5]*b[4] + a[6]*b[3] + a[7]*b[2] + a[8]*b[1] + a[9]*b[0]; - result[10] = c[10] + a[0]*b[10] + a[1]*b[9] + a[2]*b[8] + a[3]*b[7] + a[4]*b[6] + a[5]*b[5] + a[6]*b[4] + a[7]*b[3] + a[8]*b[2] + a[9]*b[1] + a[10]*b[0]; - result[11] = c[11] + a[0]*b[11] + a[1]*b[10] + a[2]*b[9] + a[3]*b[8] + a[4]*b[7] + a[5]*b[6] + a[6]*b[5] + a[7]*b[4] + a[8]*b[3] + a[9]*b[2] + a[10]*b[1] + a[11]*b[0]; - result[12] = a[1]*b[11] + a[2]*b[10] + a[3]*b[9] + a[4]*b[8] + a[5]*b[7] + a[6]*b[6] + a[7]*b[5] + a[8]*b[4] + a[9]*b[3] + a[10]*b[2] + a[11]*b[1]; - result[13] = a[2]*b[11] + a[3]*b[10] + a[4]*b[9] + a[5]*b[8] + a[6]*b[7] + a[7]*b[6] + a[8]*b[5] + a[9]*b[4] + a[10]*b[3] + a[11]*b[2]; - result[14] = a[3]*b[11] + a[4]*b[10] + a[5]*b[9] + a[6]*b[8] + a[7]*b[7] + a[8]*b[6] + a[9]*b[5] + a[10]*b[4] + a[11]*b[3]; - result[15] = a[4]*b[11] + a[5]*b[10] + a[6]*b[9] + a[7]*b[8] + a[8]*b[7] + a[9]*b[6] + a[10]*b[5] + a[11]*b[4]; - result[16] = a[5]*b[11] + a[6]*b[10] + a[7]*b[9] + a[8]*b[8] + a[9]*b[7] + a[10]*b[6] + a[11]*b[5]; - result[17] = a[6]*b[11] + a[7]*b[10] + a[8]*b[9] + a[9]*b[8] + a[10]*b[7] + a[11]*b[6]; - result[18] = a[7]*b[11] + a[8]*b[10] + a[9]*b[9] + a[10]*b[8] + a[11]*b[7]; - result[19] = a[8]*b[11] + a[9]*b[10] + a[10]*b[9] + a[11]*b[8]; - result[20] = a[9]*b[11] + a[10]*b[10] + a[11]*b[9]; - result[21] = a[10]*b[11] + a[11]*b[10]; - result[22] = a[11]*b[11]; - result[23] = 0i64; - - // Reduce limbs - UnpackedScalar::reduce_limbs(&mut result) - } - - /// Reduce 24 limbs to 12, consuming the input. Reduction is mod - /// - /// l = 2^252 + 27742317777372353535851937790883648493, - /// - /// so - /// - /// 2^252 = -27742317777372353535851937790883648493 (mod l). - /// - /// We can write the right-hand side in 21-bit limbs as - /// - /// rhs = 666643 * 2^0 - /// + 470296 * 2^21 - /// + 654183 * 2^42 - /// - 997805 * 2^63 - /// + 136657 * 2^84 - /// - 683901 * 2^105 - /// - /// The (12+k)-th limb of `limbs` is the coefficient of - /// - /// 2^(252 + 21*k) - /// - /// since 12*21 = 252. By the above, we have that - /// - /// c * 2^(252 + 21*k) = c * 666643 * 2^(21*k) - /// + c * 470296 * 2^(42*k) + ... - /// - /// so we can eliminate it by adding those values to the lower - /// limbs. Reduction mod l amounts to eliminating all of the - /// high limbs while carrying as appropriate to prevent - /// overflows in the lower limbs. - fn reduce_limbs(mut limbs: &mut [i64; 24]) -> UnpackedScalar { - #[inline] - #[allow(dead_code)] - fn do_reduction(limbs: &mut [i64; 24], i: usize) { - limbs[i - 12] += limbs[i] * 666643; - limbs[i - 11] += limbs[i] * 470296; - limbs[i - 10] += limbs[i] * 654183; - limbs[i - 9] -= limbs[i] * 997805; - limbs[i - 8] += limbs[i] * 136657; - limbs[i - 7] -= limbs[i] * 683901; - limbs[i] = 0; - } - /// Carry excess from the `i`-th limb into the `(i+1)`-th limb. - /// Postcondition: `0 <= limbs[i] < 2^21`. - #[inline] - #[allow(dead_code)] - fn do_carry_uncentered(limbs: &mut [i64; 24], i: usize) { - let carry: i64 = limbs[i] >> 21; - limbs[i+1] += carry; - limbs[i ] -= carry << 21; - } - #[inline] - #[allow(dead_code)] - /// Carry excess from the `i`-th limb into the `(i+1)`-th limb. - /// Postcondition: `-2^20 <= limbs[i] < 2^20`. - fn do_carry_centered(limbs: &mut [i64; 24], i: usize) { - let carry: i64 = (limbs[i] + (1<<20)) >> 21; - limbs[i+1] += carry; - limbs[i ] -= carry << 21; - } - - for i in 0..23 { - do_carry_centered(&mut limbs, i); - } - for i in (0..23).filter(|x| x % 2 == 1) { - do_carry_centered(&mut limbs, i); - } - - do_reduction(&mut limbs, 23); - do_reduction(&mut limbs, 22); - do_reduction(&mut limbs, 21); - do_reduction(&mut limbs, 20); - do_reduction(&mut limbs, 19); - do_reduction(&mut limbs, 18); - - for i in (6..18).filter(|x| x % 2 == 0) { - do_carry_centered(&mut limbs, i); - } - for i in (6..16).filter(|x| x % 2 == 1) { - do_carry_centered(&mut limbs, i); - } - - do_reduction(&mut limbs, 17); - do_reduction(&mut limbs, 16); - do_reduction(&mut limbs, 15); - do_reduction(&mut limbs, 14); - do_reduction(&mut limbs, 13); - do_reduction(&mut limbs, 12); - - for i in (0..12).filter(|x| x % 2 == 0) { - do_carry_centered(&mut limbs, i); - } - for i in (0..12).filter(|x| x % 2 == 1) { - do_carry_centered(&mut limbs, i); - } - - do_reduction(&mut limbs, 12); - - for i in 0..12 { - do_carry_uncentered(&mut limbs, i); - } - - do_reduction(&mut limbs, 12); - - for i in 0..11 { - do_carry_uncentered(&mut limbs, i); - } - - UnpackedScalar(*array_ref!(limbs, 0, 12)) + y.from_montgomery() } } #[cfg(test)] mod test { use super::*; + use constants; /// x = 2238329342913194256032495932344128051776374960164957527413114840482143558222 pub static X: Scalar = Scalar( @@ -815,6 +534,12 @@ mod test { 0x59, 0x13, 0xb4, 0x64, 0x1b, 0xc2, 0x7d, 0x52, 0x52, 0xa5, 0x85, 0x10, 0x1b, 0xcc, 0x42, 0x44, 0xd4, 0x49, 0xf4, 0xa8, 0x79, 0xd9, 0xf2, 0x04]); + /// 1/x = 6859937278830797291664592131120606308688036382723378951768035303146619657244 + pub static XINV: Scalar = Scalar( + [0x1c, 0xdc, 0x17, 0xfc, 0xe0, 0xe9, 0xa5, 0xbb, + 0xd9, 0x24, 0x7e, 0x56, 0xbb, 0x01, 0x63, 0x47, + 0xbb, 0xba, 0x31, 0xed, 0xd5, 0xa9, 0xbb, 0x96, + 0xd5, 0x0b, 0xcd, 0x7a, 0x3f, 0x96, 0x2a, 0x0f]); /// y = 2592331292931086675770238855846338635550719849568364935475441891787804997264 pub static Y: Scalar = Scalar( [0x90, 0x76, 0x33, 0xfe, 0x1c, 0x4b, 0x66, 0xa4, @@ -952,6 +677,7 @@ mod test { #[test] fn invert() { let inv_X = X.invert(); + assert_eq!(inv_X, XINV); let should_be_one = &inv_X * &X; assert_eq!(should_be_one, Scalar::one()); } @@ -984,7 +710,7 @@ mod bench { use test::Bencher; use super::*; - use super::test::{X, Y, Z}; + use super::test::{X}; #[bench] fn scalar_random(b: &mut Bencher) { @@ -993,28 +719,9 @@ mod bench { b.iter(|| Scalar::random(&mut csprng)); } - #[bench] - fn scalar_multiply_add(b: &mut Bencher) { - b.iter(|| Scalar::multiply_add(&X, &Y, &Z)); - } - #[bench] fn invert(b: &mut Bencher) { let x = X.unpack(); b.iter(|| x.invert()); } - - #[bench] - fn square(b: &mut Bencher) { - let x = X.unpack(); - b.iter(|| x.square()); - } - - #[bench] - fn scalar_unpacked_multiply_add(b: &mut Bencher) { - let x = X.unpack(); - let y = Y.unpack(); - let z = Z.unpack(); - b.iter(|| UnpackedScalar::multiply_add(&x, &y, &z)); - } } diff --git a/src/scalar_32bit.rs b/src/scalar_32bit.rs new file mode 100644 index 0000000..b0d8b57 --- /dev/null +++ b/src/scalar_32bit.rs @@ -0,0 +1,559 @@ +//! Arithmetic mod 2^252 + 27742317777372353535851937790883648493 +//! with 9 29-bit unsigned limbs +//! +//! To see that this is safe for intermediate results, note that +//! the largest limb in a 9 by 9 product of 29-bit limbs will be +//! (0x1fffffff^2) * 9 = 0x23fffffdc0000009 (62 bits). +//! +//! For a one level Karatsuba decomposition, the specific ranges +//! depend on how the limbs are combined, but will stay within +//! -0x1ffffffe00000008 (62 bits with sign bit) to +//! 0x43fffffbc0000011 (63 bits), which is still safe. +//! +//! (the 9th limb will never exceed 21 bits, so the actual +//! ranges are slightly smaller) + +use core::fmt::Debug; +use core::ops::{Index, IndexMut}; + +use constants; + +/// The `Scalar32` struct represents an element in ℤ/lℤ as 9 29-bit limbs +#[derive(Copy,Clone)] +pub struct Scalar32(pub [u32; 9]); + +impl Debug for Scalar32 { + fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { + write!(f, "Scalar32: {:?}", &self.0[..]) + } +} + +impl Index for Scalar32 { + type Output = u32; + fn index(&self, _index: usize) -> &u32 { + &(self.0[_index]) + } +} + +impl IndexMut for Scalar32 { + fn index_mut(&mut self, _index: usize) -> &mut u32 { + &mut (self.0[_index]) + } +} + +/// u32 * u32 = u64 multiply helper +#[inline(always)] +fn m(x: u32, y: u32) -> u64 { + (x as u64) * (y as u64) +} + +impl Scalar32 { + /// Return the zero scalar. + pub fn zero() -> Scalar32 { + Scalar32([0,0,0,0,0,0,0,0,0]) + } + + /// Unpack a 32 byte / 512 bit scalar into 9 29-bit limbs, ignoring the upper 3 bits. + pub fn from_bytes(bytes: &[u8; 32]) -> Scalar32 { + let mut words = [0u32; 8]; + for i in 0..8 { + for j in 0..4 { + words[i] |= (bytes[(i * 4) + j] as u32) << (j * 8); + } + } + + let mask = (1u32 << 29) - 1; + let top_mask = (1u32 << 21) - 1; + let mut s = Scalar32::zero(); + + s[ 0] = words[0] & mask; + s[ 1] = ((words[0] >> 29) | (words[1] << 3)) & mask; + s[ 2] = ((words[1] >> 26) | (words[2] << 6)) & mask; + s[ 3] = ((words[2] >> 23) | (words[3] << 9)) & mask; + s[ 4] = ((words[3] >> 20) | (words[4] << 12)) & mask; + s[ 5] = ((words[4] >> 17) | (words[5] << 15)) & mask; + s[ 6] = ((words[5] >> 14) | (words[6] << 18)) & mask; + s[ 7] = ((words[6] >> 11) | (words[7] << 21)) & mask; + s[ 8] = (words[7] >> 8) & top_mask; + + s + } + + /// Reduce a 64 byte / 512 bit scalar mod l. + pub fn from_bytes_wide(bytes: &[u8; 64]) -> Scalar32 { + let mut words = [0u32; 16]; + for i in 0..16 { + for j in 0..4 { + words[i] |= (bytes[(i * 4) + j] as u32) << (j * 8); + } + } + + let mask = (1u32 << 29) - 1; + let mut lo = Scalar32::zero(); + let mut hi = Scalar32::zero(); + + lo[0] = words[ 0] & mask; + lo[1] = ((words[ 0] >> 29) | (words[ 1] << 3)) & mask; + lo[2] = ((words[ 1] >> 26) | (words[ 2] << 6)) & mask; + lo[3] = ((words[ 2] >> 23) | (words[ 3] << 9)) & mask; + lo[4] = ((words[ 3] >> 20) | (words[ 4] << 12)) & mask; + lo[5] = ((words[ 4] >> 17) | (words[ 5] << 15)) & mask; + lo[6] = ((words[ 5] >> 14) | (words[ 6] << 18)) & mask; + lo[7] = ((words[ 6] >> 11) | (words[ 7] << 21)) & mask; + lo[8] = ((words[ 7] >> 8) | (words[ 8] << 24)) & mask; + hi[0] = ((words[ 8] >> 5) | (words[ 9] << 27)) & mask; + hi[1] = (words[ 9] >> 2) & mask; + hi[2] = ((words[ 9] >> 31) | (words[10] << 1)) & mask; + hi[3] = ((words[10] >> 28) | (words[11] << 4)) & mask; + hi[4] = ((words[11] >> 25) | (words[12] << 7)) & mask; + hi[5] = ((words[12] >> 22) | (words[13] << 10)) & mask; + hi[6] = ((words[13] >> 19) | (words[14] << 13)) & mask; + hi[7] = ((words[14] >> 16) | (words[15] << 16)) & mask; + hi[8] = (words[15] >> 13) & mask; + + lo = Scalar32::montgomery_mul(&lo, &constants::R); // (lo * R) / R = lo + hi = Scalar32::montgomery_mul(&hi, &constants::RR); // (hi * R^2) / R = hi * R + + Scalar32::add(&hi, &lo) // (hi * R) + lo + } + + /// Pack the limbs of this `Scalar32` into 32 bytes. + pub fn to_bytes(&self) -> [u8; 32] { + let mut s = [0u8; 32]; + + s[0] = (self.0[ 0] >> 0) as u8; + s[1] = (self.0[ 0] >> 8) as u8; + s[2] = (self.0[ 0] >> 16) as u8; + s[3] = ((self.0[ 0] >> 24) | (self.0[ 1] << 5)) as u8; + s[4] = (self.0[ 1] >> 3) as u8; + s[5] = (self.0[ 1] >> 11) as u8; + s[6] = (self.0[ 1] >> 19) as u8; + s[7] = ((self.0[ 1] >> 27) | (self.0[ 2] << 2)) as u8; + s[8] = (self.0[ 2] >> 6) as u8; + s[9] = (self.0[ 2] >> 14) as u8; + s[10] = ((self.0[ 2] >> 22) | (self.0[ 3] << 7)) as u8; + s[11] = (self.0[ 3] >> 1) as u8; + s[12] = (self.0[ 3] >> 9) as u8; + s[13] = (self.0[ 3] >> 17) as u8; + s[14] = ((self.0[ 3] >> 25) | (self.0[ 4] << 4)) as u8; + s[15] = (self.0[ 4] >> 4) as u8; + s[16] = (self.0[ 4] >> 12) as u8; + s[17] = (self.0[ 4] >> 20) as u8; + s[18] = ((self.0[ 4] >> 28) | (self.0[ 5] << 1)) as u8; + s[19] = (self.0[ 5] >> 7) as u8; + s[20] = (self.0[ 5] >> 15) as u8; + s[21] = ((self.0[ 5] >> 23) | (self.0[ 6] << 6)) as u8; + s[22] = (self.0[ 6] >> 2) as u8; + s[23] = (self.0[ 6] >> 10) as u8; + s[24] = (self.0[ 6] >> 18) as u8; + s[25] = ((self.0[ 6] >> 26) | (self.0[ 7] << 3)) as u8; + s[26] = (self.0[ 7] >> 5) as u8; + s[27] = (self.0[ 7] >> 13) as u8; + s[28] = (self.0[ 7] >> 21) as u8; + s[29] = (self.0[ 8] >> 0) as u8; + s[30] = (self.0[ 8] >> 8) as u8; + s[31] = (self.0[ 8] >> 16) as u8; + + s + } + + /// Compute `a + b` (mod l). + pub fn add(a: &Scalar32, b: &Scalar32) -> Scalar32 { + let mut sum = Scalar32::zero(); + let mask = (1u32 << 29) - 1; + + // a + b + let mut carry: u32 = 0; + for i in 0..9 { + carry = a[i] + b[i] + (carry >> 29); + sum[i] = carry & mask; + } + + // subtract l if the sum is >= l + Scalar32::sub(&sum, &constants::L) + } + + /// Compute `a - b` (mod l). + pub fn sub(a: &Scalar32, b: &Scalar32) -> Scalar32 { + let mut difference = Scalar32::zero(); + let mask = (1u32 << 29) - 1; + + // a - b + let mut borrow: u32 = 0; + for i in 0..9 { + borrow = a[i].wrapping_sub(b[i] + (borrow >> 31)); + difference[i] = borrow & mask; + } + + // conditionally add l if the difference is negative + let underflow_mask = ((borrow >> 31) ^ 1).wrapping_sub(1); + let mut carry: u32 = 0; + for i in 0..9 { + carry = (carry >> 29) + difference[i] + (constants::L[i] & underflow_mask); + difference[i] = carry & mask; + } + + difference + } + + /// Compute `a * b`. + /// + /// This is implemented with a one-level refined Karatsuba decomposition + #[inline(always)] + fn mul_internal(a: &Scalar32, b: &Scalar32) -> [u64; 17] { + let mut z = [0u64; 17]; + + z[0] = m(a[0],b[0]); // c00 + z[1] = m(a[0],b[1]) + m(a[1],b[0]); // c01 + z[2] = m(a[0],b[2]) + m(a[1],b[1]) + m(a[2],b[0]); // c02 + z[3] = m(a[0],b[3]) + m(a[1],b[2]) + m(a[2],b[1]) + m(a[3],b[0]); // c03 + z[4] = m(a[0],b[4]) + m(a[1],b[3]) + m(a[2],b[2]) + m(a[3],b[1]) + m(a[4],b[0]); // c04 + z[5] = m(a[1],b[4]) + m(a[2],b[3]) + m(a[3],b[2]) + m(a[4],b[1]); // c05 + z[6] = m(a[2],b[4]) + m(a[3],b[3]) + m(a[4],b[2]); // c06 + z[7] = m(a[3],b[4]) + m(a[4],b[3]); // c07 + z[8] = (m(a[4],b[4])).wrapping_sub(z[3]); // c08 - c03 + + z[10] = z[5].wrapping_sub(m(a[5],b[5])); // c05mc10 + z[11] = z[6].wrapping_sub(m(a[5],b[6]) + m(a[6],b[5])); // c06mc11 + z[12] = z[7].wrapping_sub(m(a[5],b[7]) + m(a[6],b[6]) + m(a[7],b[5])); // c07mc12 + z[13] = m(a[5],b[8]) + m(a[6],b[7]) + m(a[7],b[6]) + m(a[8],b[5]); // c13 + z[14] = m(a[6],b[8]) + m(a[7],b[7]) + m(a[8],b[6]); // c14 + z[15] = m(a[7],b[8]) + m(a[8],b[7]); // c15 + z[16] = m(a[8],b[8]); // c16 + + z[ 5] = z[10].wrapping_sub(z[ 0]); // c05mc10 - c00 + z[ 6] = z[11].wrapping_sub(z[ 1]); // c06mc11 - c01 + z[ 7] = z[12].wrapping_sub(z[ 2]); // c07mc12 - c02 + z[ 8] = z[ 8].wrapping_sub(z[13]); // c08mc13 - c03 + z[ 9] = z[14].wrapping_add(z[ 4]); // c14 + c04 + z[10] = z[15].wrapping_add(z[10]); // c15 + c05mc10 + z[11] = z[16].wrapping_add(z[11]); // c16 + c06mc11 + + let aa = [ + a[0]+a[5], + a[1]+a[6], + a[2]+a[7], + a[3]+a[8] + ]; + + let bb = [ + b[0]+b[5], + b[1]+b[6], + b[2]+b[7], + b[3]+b[8] + ]; + + z[ 5] = (m(aa[0],bb[0])) .wrapping_add(z[ 5]); // c20 + c05mc10 - c00 + z[ 6] = (m(aa[0],bb[1]) + m(aa[1],bb[0])) .wrapping_add(z[ 6]); // c21 + c06mc11 - c01 + z[ 7] = (m(aa[0],bb[2]) + m(aa[1],bb[1]) + m(aa[2],bb[0])) .wrapping_add(z[ 7]); // c22 + c07mc12 - c02 + z[ 8] = (m(aa[0],bb[3]) + m(aa[1],bb[2]) + m(aa[2],bb[1]) + m(aa[3],bb[0])) .wrapping_add(z[ 8]); // c23 + c08mc13 - c03 + z[ 9] = (m(aa[0], b[4]) + m(aa[1],bb[3]) + m(aa[2],bb[2]) + m(aa[3],bb[1]) + m(a[4],bb[0])).wrapping_sub(z[ 9]); // c24 - c14 - c04 + z[10] = ( m(aa[1], b[4]) + m(aa[2],bb[3]) + m(aa[3],bb[2]) + m(a[4],bb[1])).wrapping_sub(z[10]); // c25 - c15 - c05mc10 + z[11] = ( m(aa[2], b[4]) + m(aa[3],bb[3]) + m(a[4],bb[2])).wrapping_sub(z[11]); // c26 - c16 - c06mc11 + z[12] = ( m(aa[3], b[4]) + m(a[4],bb[3])).wrapping_sub(z[12]); // c27 - c07mc12 + + z + } + + /// Compute `a^2`. + #[inline(always)] + fn square_internal(a: &Scalar32) -> [u64; 17] { + let aa = [ + a[0]*2, + a[1]*2, + a[2]*2, + a[3]*2, + a[4]*2, + a[5]*2, + a[6]*2, + a[7]*2 + ]; + + [ + m( a[0],a[0]), + m(aa[0],a[1]), + m(aa[0],a[2]) + m( a[1],a[1]), + m(aa[0],a[3]) + m(aa[1],a[2]), + m(aa[0],a[4]) + m(aa[1],a[3]) + m( a[2],a[2]), + m(aa[0],a[5]) + m(aa[1],a[4]) + m(aa[2],a[3]), + m(aa[0],a[6]) + m(aa[1],a[5]) + m(aa[2],a[4]) + m( a[3],a[3]), + m(aa[0],a[7]) + m(aa[1],a[6]) + m(aa[2],a[5]) + m(aa[3],a[4]), + m(aa[0],a[8]) + m(aa[1],a[7]) + m(aa[2],a[6]) + m(aa[3],a[5]) + m( a[4],a[4]), + m(aa[1],a[8]) + m(aa[2],a[7]) + m(aa[3],a[6]) + m(aa[4],a[5]), + m(aa[2],a[8]) + m(aa[3],a[7]) + m(aa[4],a[6]) + m( a[5],a[5]), + m(aa[3],a[8]) + m(aa[4],a[7]) + m(aa[5],a[6]), + m(aa[4],a[8]) + m(aa[5],a[7]) + m( a[6],a[6]), + m(aa[5],a[8]) + m(aa[6],a[7]), + m(aa[6],a[8]) + m( a[7],a[7]), + m(aa[7],a[8]), + m( a[8],a[8]), + ] + } + + /// Compute `limbs/R` (mod l), where R is the Montgomery modulus 2^261 + #[inline(always)] + fn montgomery_reduce(limbs: &[u64; 17]) -> Scalar32 { + + #[inline(always)] + fn part1(sum: u64) -> (u64, u32) { + let p = (sum as u32).wrapping_mul(constants::LFACTOR) & ((1u32 << 29) - 1); + ((sum + m(p,constants::L[0])) >> 29, p) + } + + #[inline(always)] + fn part2(sum: u64) -> (u64, u32) { + let w = (sum as u32) & ((1u32 << 29) - 1); + (sum >> 29, w) + } + + // note: l5,l6,l7 are zero, so their multiplies can be skipped + let l = &constants::L; + + // the first half computes the Montgomery adjustment factor n, and begins adding n*l to make limbs divisible by R + let (carry, n0) = part1( limbs[ 0]); + let (carry, n1) = part1(carry + limbs[ 1] + m(n0,l[1])); + let (carry, n2) = part1(carry + limbs[ 2] + m(n0,l[2]) + m(n1,l[1])); + let (carry, n3) = part1(carry + limbs[ 3] + m(n0,l[3]) + m(n1,l[2]) + m(n2,l[1])); + let (carry, n4) = part1(carry + limbs[ 4] + m(n0,l[4]) + m(n1,l[3]) + m(n2,l[2]) + m(n3,l[1])); + let (carry, n5) = part1(carry + limbs[ 5] + m(n1,l[4]) + m(n2,l[3]) + m(n3,l[2]) + m(n4,l[1])); + let (carry, n6) = part1(carry + limbs[ 6] + m(n2,l[4]) + m(n3,l[3]) + m(n4,l[2]) + m(n5,l[1])); + let (carry, n7) = part1(carry + limbs[ 7] + m(n3,l[4]) + m(n4,l[3]) + m(n5,l[2]) + m(n6,l[1])); + let (carry, n8) = part1(carry + limbs[ 8] + m(n0,l[8]) + m(n4,l[4]) + m(n5,l[3]) + m(n6,l[2]) + m(n7,l[1])); + + // limbs is divisible by R now, so we can divide by R by simply storing the upper half as the result + let (carry, r0) = part2(carry + limbs[ 9] + m(n1,l[8]) + m(n5,l[4]) + m(n6,l[3]) + m(n7,l[2]) + m(n8,l[1])); + let (carry, r1) = part2(carry + limbs[10] + m(n2,l[8]) + m(n6,l[4]) + m(n7,l[3]) + m(n8,l[2])); + let (carry, r2) = part2(carry + limbs[11] + m(n3,l[8]) + m(n7,l[4]) + m(n8,l[3])); + let (carry, r3) = part2(carry + limbs[12] + m(n4,l[8]) + m(n8,l[4])); + let (carry, r4) = part2(carry + limbs[13] + m(n5,l[8]) ); + let (carry, r5) = part2(carry + limbs[14] + m(n6,l[8]) ); + let (carry, r6) = part2(carry + limbs[15] + m(n7,l[8]) ); + let (carry, r7) = part2(carry + limbs[16] + m(n8,l[8])); + let r8 = carry as u32; + + // result may be >= l, so attempt to subtract l + Scalar32::sub(&Scalar32([r0,r1,r2,r3,r4,r5,r6,r7,r8]), l) + } + + /// Compute `a * b` (mod l). + #[inline(never)] + pub fn mul(a: &Scalar32, b: &Scalar32) -> Scalar32 { + let ab = Scalar32::montgomery_reduce(&Scalar32::mul_internal(a, b)); + Scalar32::montgomery_reduce(&Scalar32::mul_internal(&ab, &constants::RR)) + } + + /// Compute `a^2` (mod l). + #[inline(never)] + pub fn square(&self) -> Scalar32 { + let aa = Scalar32::montgomery_reduce(&Scalar32::square_internal(self)); + Scalar32::montgomery_reduce(&Scalar32::mul_internal(&aa, &constants::RR)) + } + + /// Compute `(a * b) / R` (mod l), where R is the Montgomery modulus 2^261 + #[inline(never)] + pub fn montgomery_mul(a: &Scalar32, b: &Scalar32) -> Scalar32 { + Scalar32::montgomery_reduce(&Scalar32::mul_internal(a, b)) + } + + /// Compute `(a^2) / R` (mod l) in Montgomery form, where R is the Montgomery modulus 2^261 + #[inline(never)] + pub fn montgomery_square(&self) -> Scalar32 { + Scalar32::montgomery_reduce(&Scalar32::square_internal(self)) + } + + /// Puts a Scalar32 in to Montgomery form, i.e. computes `a*R (mod l)` + #[inline(never)] + pub fn to_montgomery(&self) -> Scalar32 { + Scalar32::montgomery_mul(self, &constants::RR) + } + + /// Takes a Scalar32 out of Montgomery form, i.e. computes `a/R (mod l)` + pub fn from_montgomery(&self) -> Scalar32 { + let mut limbs = [0u64; 17]; + for i in 0..9 { + limbs[i] = self[i] as u64; + } + Scalar32::montgomery_reduce(&limbs) + } +} + + +#[cfg(test)] +mod test { + use super::*; + + /// Note: x is 2^253-1 which is slightly larger than the largest scalar produced by + /// this implementation (l-1), and should verify there are no overflows for valid scalars + /// + /// x = 2^253-1 = 14474011154664524427946373126085988481658748083205070504932198000989141204991 + /// x = 7237005577332262213973186563042994240801631723825162898930247062703686954002 mod l + /// x = 5147078182513738803124273553712992179887200054963030844803268920753008712037*R mod l in Montgomery form + pub static X: Scalar32 = Scalar32( + [0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff, + 0x1fffffff, 0x1fffffff, 0x1fffffff, 0x1fffffff, + 0x001fffff]); + + /// x^2 = 3078544782642840487852506753550082162405942681916160040940637093560259278169 mod l + pub static XX: Scalar32 = Scalar32( + [0x00217559, 0x000b3401, 0x103ff43b, 0x1462a62c, + 0x1d6f9f38, 0x18e7a42f, 0x09a3dcee, 0x008dbe18, + 0x0006ce65]); + + /// x^2 = 2912514428060642753613814151688322857484807845836623976981729207238463947987*R mod l in Montgomery form + pub static XX_MONT: Scalar32 = Scalar32( + [0x152b4d2e, 0x0571d53b, 0x1da6d964, 0x188663b6, + 0x1d1b5f92, 0x19d50e3f, 0x12306c29, 0x0c6f26fe, + 0x00030edb]); + + /// y = 6145104759870991071742105800796537629880401874866217824609283457819451087098 + pub static Y: Scalar32 = Scalar32( + [0x1e1458fa, 0x165ba838, 0x1d787b36, 0x0e577f3a, + 0x1d2baf06, 0x1d689a19, 0x1fff3047, 0x117704ab, + 0x000d9601]); + + /// x*y = 36752150652102274958925982391442301741 + pub static XY: Scalar32 = Scalar32( + [0x0ba7632d, 0x017736bb, 0x15c76138, 0x0c69daa1, + 0x000001ba, 0x00000000, 0x00000000, 0x00000000, + 0x00000000]); + + /// x*y = 3783114862749659543382438697751927473898937741870308063443170013240655651591*R mod l in Montgomery form + pub static XY_MONT: Scalar32 = Scalar32( + [0x077b51e1, 0x1c64e119, 0x02a19ef5, 0x18d2129e, + 0x00de0430, 0x045a7bc8, 0x04cfc7c9, 0x1c002681, + 0x000bdc1c]); + + /// a = 2351415481556538453565687241199399922945659411799870114962672658845158063753 + pub static A: Scalar32 = Scalar32( + [0x07b3be89, 0x02291b60, 0x14a99f03, 0x07dc3787, + 0x0a782aae, 0x16262525, 0x0cfdb93f, 0x13f5718d, + 0x000532da]); + + /// b = 4885590095775723760407499321843594317911456947580037491039278279440296187236 + pub static B: Scalar32 = Scalar32( + [0x15421564, 0x1e69fd72, 0x093d9692, 0x161785be, + 0x1587d69f, 0x09d9dada, 0x130246c0, 0x0c0a8e72, + 0x000acd25]); + + /// a+b = 0 + /// a-b = 4702830963113076907131374482398799845891318823599740229925345317690316127506 + pub static AB: Scalar32 = Scalar32( + [0x0f677d12, 0x045236c0, 0x09533e06, 0x0fb86f0f, + 0x14f0555c, 0x0c4c4a4a, 0x19fb727f, 0x07eae31a, + 0x000a65b5]); + + // c = (2^512 - 1) % l = 1627715501170711445284395025044413883736156588369414752970002579683115011840 + pub static C: Scalar32 = Scalar32( + [0x049c0f00, 0x00308f1a, 0x0164d1e9, 0x1c374ed1, + 0x1be65d00, 0x19e90bfa, 0x08f73bb1, 0x036f8613, + 0x00039941]); + + #[test] + fn mul_max() { + let res = Scalar32::mul(&X, &X); + for i in 0..9 { + assert!(res[i] == XX[i]); + } + } + + #[test] + fn square_max() { + let res = X.square(); + for i in 0..9 { + assert!(res[i] == XX[i]); + } + } + + #[test] + fn montgomery_mul_max() { + let res = Scalar32::montgomery_mul(&X, &X); + for i in 0..9 { + assert!(res[i] == XX_MONT[i]); + } + } + + #[test] + fn montgomery_square_max() { + let res = X.montgomery_square(); + for i in 0..9 { + assert!(res[i] == XX_MONT[i]); + } + } + + #[test] + fn mul() { + let res = Scalar32::mul(&X, &Y); + for i in 0..9 { + assert!(res[i] == XY[i]); + } + } + + #[test] + fn montgomery_mul() { + let res = Scalar32::montgomery_mul(&X, &Y); + for i in 0..9 { + assert!(res[i] == XY_MONT[i]); + } + } + + #[test] + fn add() { + let res = Scalar32::add(&A, &B); + let zero = Scalar32::zero(); + for i in 0..9 { + assert!(res[i] == zero[i]); + } + } + + #[test] + fn sub() { + let res = Scalar32::sub(&A, &B); + for i in 0..9 { + assert!(res[i] == AB[i]); + } + } + + #[test] + fn from_bytes_wide() { + let bignum = [255u8; 64]; // 2^512 - 1 + let reduced = Scalar32::from_bytes_wide(&bignum); + for i in 0..9 { + assert!(reduced[i] == C[i]); + } + } +} + + +#[cfg(all(test, feature = "bench"))] +mod bench { + use test::Bencher; + + use super::*; + use super::test::{X, Y}; + + #[bench] + fn square(b: &mut Bencher) { + b.iter(|| X.square()); + } + + #[bench] + fn mul(b: &mut Bencher) { + b.iter(|| Scalar32::mul(&X, &Y)); + } + + #[bench] + fn montgomery_square(b: &mut Bencher) { + b.iter(|| X.montgomery_square()); + } + + #[bench] + fn montgomery_mul(b: &mut Bencher) { + b.iter(|| Scalar32::montgomery_mul(&X, &Y)); + } + + #[bench] + fn from_bytes_wide(b: &mut Bencher) { + let bignum = [255u8; 64]; // 2^512 - 1 + b.iter(|| Scalar32::from_bytes_wide(&bignum)); + } +} diff --git a/src/scalar_64bit.rs b/src/scalar_64bit.rs new file mode 100644 index 0000000..7076b56 --- /dev/null +++ b/src/scalar_64bit.rs @@ -0,0 +1,475 @@ +//! Arithmetic mod 2^252 + 27742317777372353535851937790883648493 +//! with 5 52-bit unsigned limbs. 51-bit limbs would cover the +//! desired bit range (253 bits), but isn't large enough to reduce +//! a 512 bit number with Montgomery multiplication, so 52 bits is +//! used instead +//! +//! To see that this is safe for intermediate results, note that +//! the largest limb in a 5 by 5 product of 52-bit limbs will be +//! (0xfffffffffffff^2) * 5 = 0x4ffffffffffff60000000000005 (107 bits). +//! +//! (the 5th limb will never exceed 45 bits, so the actual +//! ranges are slightly smaller) + + +use core::fmt::Debug; +use core::ops::{Index, IndexMut}; + +use constants; + +/// The `Scalar64` struct represents an element in ℤ/lℤ as 5 52-bit limbs +#[derive(Copy,Clone)] +pub struct Scalar64(pub [u64; 5]); + +impl Debug for Scalar64 { + fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { + write!(f, "Scalar64: {:?}", &self.0[..]) + } +} + +impl Index for Scalar64 { + type Output = u64; + fn index(&self, _index: usize) -> &u64 { + &(self.0[_index]) + } +} + +impl IndexMut for Scalar64 { + fn index_mut(&mut self, _index: usize) -> &mut u64 { + &mut (self.0[_index]) + } +} + +/// u64 * u64 = u128 multiply helper +#[inline(always)] +fn m(x: u64, y: u64) -> u128 { + (x as u128) * (y as u128) +} + +impl Scalar64 { + /// Return the zero scalar + pub fn zero() -> Scalar64 { + Scalar64([0,0,0,0,0]) + } + + /// Unpack a 32 byte / 256 bit scalar into 5 52-bit limbs, ignoring the upper 3 bits + pub fn from_bytes(bytes: &[u8; 32]) -> Scalar64 { + let mut words = [0u64; 8]; + for i in 0..4 { + for j in 0..8 { + words[i] |= (bytes[(i * 8) + j] as u64) << (j * 8); + } + } + + let mask = (1u64 << 52) - 1; + let top_mask = (1u64 << 45) - 1; + let mut s = Scalar64::zero(); + + s[ 0] = words[0] & mask; + s[ 1] = ((words[0] >> 52) | (words[1] << 12)) & mask; + s[ 2] = ((words[1] >> 40) | (words[2] << 24)) & mask; + s[ 3] = ((words[2] >> 28) | (words[3] << 36)) & mask; + s[ 4] = (words[3] >> 16) & top_mask; + + s + } + + /// Reduce a 64 byte / 512 bit scalar mod l + pub fn from_bytes_wide(bytes: &[u8; 64]) -> Scalar64 { + let mut words = [064; 16]; + for i in 0..8 { + for j in 0..8 { + words[i] |= (bytes[(i * 8) + j] as u64) << (j * 8); + } + } + + let mask = (1u64 << 52) - 1; + let mut lo = Scalar64::zero(); + let mut hi = Scalar64::zero(); + + lo[0] = words[ 0] & mask; + lo[1] = ((words[ 0] >> 52) | (words[ 1] << 12)) & mask; + lo[2] = ((words[ 1] >> 40) | (words[ 2] << 24)) & mask; + lo[3] = ((words[ 2] >> 28) | (words[ 3] << 36)) & mask; + lo[4] = ((words[ 3] >> 16) | (words[ 4] << 48)) & mask; + hi[0] = (words[ 4] >> 4) & mask; + hi[1] = ((words[ 4] >> 56) | (words[ 5] << 8)) & mask; + hi[2] = ((words[ 5] >> 44) | (words[ 6] << 20)) & mask; + hi[3] = ((words[ 6] >> 32) | (words[ 7] << 32)) & mask; + hi[4] = words[ 7] >> 20 ; + + lo = Scalar64::montgomery_mul(&lo, &constants::R); // (lo * R) / R = lo + hi = Scalar64::montgomery_mul(&hi, &constants::RR); // (hi * R^2) / R = hi * R + + Scalar64::add(&hi, &lo) + } + + /// Pack the limbs of this `Scalar64` into 32 bytes + pub fn to_bytes(&self) -> [u8; 32] { + let mut s = [0u8; 32]; + + s[0] = (self.0[ 0] >> 0) as u8; + s[1] = (self.0[ 0] >> 8) as u8; + s[2] = (self.0[ 0] >> 16) as u8; + s[3] = (self.0[ 0] >> 24) as u8; + s[4] = (self.0[ 0] >> 32) as u8; + s[5] = (self.0[ 0] >> 40) as u8; + s[6] = ((self.0[ 0] >> 48) | (self.0[ 1] << 4)) as u8; + s[7] = (self.0[ 1] >> 4) as u8; + s[8] = (self.0[ 1] >> 12) as u8; + s[9] = (self.0[ 1] >> 20) as u8; + s[10] = (self.0[ 1] >> 28) as u8; + s[11] = (self.0[ 1] >> 36) as u8; + s[12] = (self.0[ 1] >> 44) as u8; + s[13] = (self.0[ 2] >> 0) as u8; + s[14] = (self.0[ 2] >> 8) as u8; + s[15] = (self.0[ 2] >> 16) as u8; + s[16] = (self.0[ 2] >> 24) as u8; + s[17] = (self.0[ 2] >> 32) as u8; + s[18] = (self.0[ 2] >> 40) as u8; + s[19] = ((self.0[ 2] >> 48) | (self.0[ 3] << 4)) as u8; + s[20] = (self.0[ 3] >> 4) as u8; + s[21] = (self.0[ 3] >> 12) as u8; + s[22] = (self.0[ 3] >> 20) as u8; + s[23] = (self.0[ 3] >> 28) as u8; + s[24] = (self.0[ 3] >> 36) as u8; + s[25] = (self.0[ 3] >> 44) as u8; + s[26] = (self.0[ 4] >> 0) as u8; + s[27] = (self.0[ 4] >> 8) as u8; + s[28] = (self.0[ 4] >> 16) as u8; + s[29] = (self.0[ 4] >> 24) as u8; + s[30] = (self.0[ 4] >> 32) as u8; + s[31] = (self.0[ 4] >> 40) as u8; + + s + } + + /// Compute `a + b` (mod l) + pub fn add(a: &Scalar64, b: &Scalar64) -> Scalar64 { + let mut sum = Scalar64::zero(); + let mask = (1u64 << 52) - 1; + + // a + b + let mut carry: u64 = 0; + for i in 0..5 { + carry = a[i] + b[i] + (carry >> 52); + sum[i] = carry & mask; + } + + // subtract l if the sum is >= l + Scalar64::sub(&sum, &constants::L) + } + + /// Compute `a - b` (mod l) + pub fn sub(a: &Scalar64, b: &Scalar64) -> Scalar64 { + let mut difference = Scalar64::zero(); + let mask = (1u64 << 52) - 1; + + // a - b + let mut borrow: u64 = 0; + for i in 0..5 { + borrow = a[i].wrapping_sub(b[i] + (borrow >> 63)); + difference[i] = borrow & mask; + } + + // conditionally add l if the difference is negative + let underflow_mask = ((borrow >> 63) ^ 1).wrapping_sub(1); + let mut carry: u64 = 0; + for i in 0..5 { + carry = (carry >> 52) + difference[i] + (constants::L[i] & underflow_mask); + difference[i] = carry & mask; + } + + difference + } + + /// Compute `a * b` + #[inline(always)] + fn mul_internal(a: &Scalar64, b: &Scalar64) -> [u128; 9] { + [ + m(a[0],b[0]), + m(a[0],b[1]) + m(a[1],b[0]), + m(a[0],b[2]) + m(a[1],b[1]) + m(a[2],b[0]), + m(a[0],b[3]) + m(a[1],b[2]) + m(a[2],b[1]) + m(a[3],b[0]), + m(a[0],b[4]) + m(a[1],b[3]) + m(a[2],b[2]) + m(a[3],b[1]) + m(a[4],b[0]), + m(a[1],b[4]) + m(a[2],b[3]) + m(a[3],b[2]) + m(a[4],b[1]), + m(a[2],b[4]) + m(a[3],b[3]) + m(a[4],b[2]), + m(a[3],b[4]) + m(a[4],b[3]), + m(a[4],b[4]) + ] + } + + /// Compute `a^2` + #[inline(always)] + fn square_internal(a: &Scalar64) -> [u128; 9] { + let aa = [ + a[0]*2, + a[1]*2, + a[2]*2, + a[3]*2, + ]; + + [ + m( a[0],a[0]), + m(aa[0],a[1]), + m(aa[0],a[2]) + m( a[1],a[1]), + m(aa[0],a[3]) + m(aa[1],a[2]), + m(aa[0],a[4]) + m(aa[1],a[3]) + m( a[2],a[2]), + m(aa[1],a[4]) + m(aa[2],a[3]), + m(aa[2],a[4]) + m( a[3],a[3]), + m(aa[3],a[4]), + m(a[4],a[4]) + ] + } + + /// Compute `limbs/R` (mod l), where R is the Montgomery modulus 2^260 + #[inline(always)] + fn montgomery_reduce(limbs: &[u128; 9]) -> Scalar64 { + + #[inline(always)] + fn part1(sum: u128) -> (u128, u64) { + let p = (sum as u64).wrapping_mul(constants::LFACTOR) & ((1u64 << 52) - 1); + ((sum + m(p,constants::L[0])) >> 52, p) + } + + #[inline(always)] + fn part2(sum: u128) -> (u128, u64) { + let w = (sum as u64) & ((1u64 << 52) - 1); + (sum >> 52, w) + } + + // note: l3 is zero, so its multiplies can be skipped + let l = &constants::L; + + // the first half computes the Montgomery adjustment factor n, and begins adding n*l to make limbs divisible by R + let (carry, n0) = part1( limbs[0]); + let (carry, n1) = part1(carry + limbs[1] + m(n0,l[1])); + let (carry, n2) = part1(carry + limbs[2] + m(n0,l[2]) + m(n1,l[1])); + let (carry, n3) = part1(carry + limbs[3] + m(n1,l[2]) + m(n2,l[1])); + let (carry, n4) = part1(carry + limbs[4] + m(n0,l[4]) + m(n2,l[2]) + m(n3,l[1])); + + // limbs is divisible by R now, so we can divide by R by simply storing the upper half as the result + let (carry, r0) = part2(carry + limbs[5] + m(n1,l[4]) + m(n3,l[2]) + m(n4,l[1])); + let (carry, r1) = part2(carry + limbs[6] + m(n2,l[4]) + m(n4,l[2])); + let (carry, r2) = part2(carry + limbs[7] + m(n3,l[4]) ); + let (carry, r3) = part2(carry + limbs[8] + m(n4,l[4])); + let r4 = carry as u64; + + // result may be >= l, so attempt to subtract l + Scalar64::sub(&Scalar64([r0,r1,r2,r3,r4]), l) + } + + /// Compute `a * b` (mod l) + #[inline(never)] + pub fn mul(a: &Scalar64, b: &Scalar64) -> Scalar64 { + let ab = Scalar64::montgomery_reduce(&Scalar64::mul_internal(a, b)); + Scalar64::montgomery_reduce(&Scalar64::mul_internal(&ab, &constants::RR)) + } + + /// Compute `a^2` (mod l) + #[inline(never)] + pub fn square(&self) -> Scalar64 { + let aa = Scalar64::montgomery_reduce(&Scalar64::square_internal(self)); + Scalar64::montgomery_reduce(&Scalar64::mul_internal(&aa, &constants::RR)) + } + + /// Compute `(a * b) / R` (mod l), where R is the Montgomery modulus 2^260 + #[inline(never)] + pub fn montgomery_mul(a: &Scalar64, b: &Scalar64) -> Scalar64 { + Scalar64::montgomery_reduce(&Scalar64::mul_internal(a, b)) + } + + /// Compute `(a^2) / R` (mod l) in Montgomery form, where R is the Montgomery modulus 2^260 + #[inline(never)] + pub fn montgomery_square(&self) -> Scalar64 { + Scalar64::montgomery_reduce(&Scalar64::square_internal(self)) + } + + /// Puts a Scalar64 in to Montgomery form, i.e. computes `a*R (mod l)` + #[inline(never)] + pub fn to_montgomery(&self) -> Scalar64 { + Scalar64::montgomery_mul(self, &constants::RR) + } + + /// Takes a Scalar64 out of Montgomery form, i.e. computes `a/R (mod l)` + #[inline(never)] + pub fn from_montgomery(&self) -> Scalar64 { + let mut limbs = [0u128; 9]; + for i in 0..5 { + limbs[i] = self[i] as u128; + } + Scalar64::montgomery_reduce(&limbs) + } +} + + +#[cfg(test)] +mod test { + use super::*; + + /// Note: x is 2^253-1 which is slightly larger than the largest scalar produced by + /// this implementation (l-1), and should show there are no overflows for valid scalars + /// + /// x = 14474011154664524427946373126085988481658748083205070504932198000989141204991 + /// x = 7237005577332262213973186563042994240801631723825162898930247062703686954002 mod l + /// x = 3057150787695215392275360544382990118917283750546154083604586903220563173085*R mod l in Montgomery form + pub static X: Scalar64 = Scalar64( + [0x000fffffffffffff, 0x000fffffffffffff, 0x000fffffffffffff, 0x000fffffffffffff, + 0x00001fffffffffff]); + + /// x^2 = 3078544782642840487852506753550082162405942681916160040940637093560259278169 mod l + pub static XX: Scalar64 = Scalar64( + [0x0001668020217559, 0x000531640ffd0ec0, 0x00085fd6f9f38a31, 0x000c268f73bb1cf4, + 0x000006ce65046df0]); + + /// x^2 = 4413052134910308800482070043710297189082115023966588301924965890668401540959*R mod l in Montgomery form + pub static XX_MONT: Scalar64 = Scalar64( + [0x000c754eea569a5c, 0x00063b6ed36cb215, 0x0008ffa36bf25886, 0x000e9183614e7543, + 0x0000061db6c6f26f]); + + /// y = 6145104759870991071742105800796537629880401874866217824609283457819451087098 + pub static Y: Scalar64 = Scalar64( + [0x000b75071e1458fa, 0x000bf9d75e1ecdac, 0x000433d2baf0672b, 0x0005fffcc11fad13, + 0x00000d96018bb825]); + + /// x*y = 36752150652102274958925982391442301741 mod l + pub static XY: Scalar64 = Scalar64( + [0x000ee6d76ba7632d, 0x000ed50d71d84e02, 0x00000000001ba634, 0x0000000000000000, + 0x0000000000000000]); + + /// x*y = 658448296334113745583381664921721413881518248721417041768778176391714104386*R mod l in Montgomery form + pub static XY_MONT: Scalar64 = Scalar64( + [0x0006d52bf200cfd5, 0x00033fb1d7021570, 0x000f201bc07139d8, 0x0001267e3e49169e, + 0x000007b839c00268]); + + /// a = 2351415481556538453565687241199399922945659411799870114962672658845158063753 + pub static A: Scalar64 = Scalar64( + [0x0005236c07b3be89, 0x0001bc3d2a67c0c4, 0x000a4aa782aae3ee, 0x0006b3f6e4fec4c4, + 0x00000532da9fab8c]); + + /// b = 4885590095775723760407499321843594317911456947580037491039278279440296187236 + pub static B: Scalar64 = Scalar64( + [0x000d3fae55421564, 0x000c2df24f65a4bc, 0x0005b5587d69fb0b, 0x00094c091b013b3b, + 0x00000acd25605473]); + + /// a+b = 0 + /// a-b = 4702830963113076907131374482398799845891318823599740229925345317690316127506 + pub static AB: Scalar64 = Scalar64( + [0x000a46d80f677d12, 0x0003787a54cf8188, 0x0004954f0555c7dc, 0x000d67edc9fd8989, + 0x00000a65b53f5718]); + + // c = (2^512 - 1) % l = 1627715501170711445284395025044413883736156588369414752970002579683115011840 + pub static C: Scalar64 = Scalar64( + [0x000611e3449c0f00, 0x000a768859347a40, 0x0007f5be65d00e1b, 0x0009a3dceec73d21, + 0x00000399411b7c30]); + + #[test] + fn mul_max() { + let res = Scalar64::mul(&X, &X); + for i in 0..5 { + assert!(res[i] == XX[i]); + } + } + + #[test] + fn square_max() { + let res = X.square(); + for i in 0..5 { + assert!(res[i] == XX[i]); + } + } + + #[test] + fn montgomery_mul_max() { + let res = Scalar64::montgomery_mul(&X, &X); + for i in 0..5 { + assert!(res[i] == XX_MONT[i]); + } + } + + #[test] + fn montgomery_square_max() { + let res = X.montgomery_square(); + for i in 0..5 { + assert!(res[i] == XX_MONT[i]); + } + } + + #[test] + fn mul() { + let res = Scalar64::mul(&X, &Y); + for i in 0..5 { + assert!(res[i] == XY[i]); + } + } + + #[test] + fn montgomery_mul() { + let res = Scalar64::montgomery_mul(&X, &Y); + for i in 0..5 { + assert!(res[i] == XY_MONT[i]); + } + } + + #[test] + fn add() { + let res = Scalar64::add(&A, &B); + let zero = Scalar64::zero(); + for i in 0..5 { + assert!(res[i] == zero[i]); + } + } + + #[test] + fn sub() { + let res = Scalar64::sub(&A, &B); + for i in 0..5 { + assert!(res[i] == AB[i]); + } + } + + #[test] + fn from_bytes_wide() { + let bignum = [255u8; 64]; // 2^512 - 1 + let reduced = Scalar64::from_bytes_wide(&bignum); + println!("{:?}", reduced); + for i in 0..5 { + assert!(reduced[i] == C[i]); + } + } +} + + +#[cfg(all(test, feature = "bench"))] +mod bench { + use test::Bencher; + + use super::*; + use super::test::{X, Y}; + + #[bench] + fn square(b: &mut Bencher) { + b.iter(|| X.square()); + } + + #[bench] + fn mul(b: &mut Bencher) { + b.iter(|| Scalar64::mul(&X, &Y)); + } + + #[bench] + fn montgomery_square(b: &mut Bencher) { + b.iter(|| X.montgomery_square()); + } + + #[bench] + fn montgomery_mul(b: &mut Bencher) { + b.iter(|| Scalar64::montgomery_mul(&X, &Y)); + } + + #[bench] + fn from_bytes_wide(b: &mut Bencher) { + let bignum = [255u8; 64]; // 2^512 - 1 + b.iter(|| Scalar64::from_bytes_wide(&bignum)); + } +} From ecef4d836ea7f71cfbe3e39b91eb3286dabf8eb7 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Sun, 1 Oct 2017 23:36:57 +0000 Subject: [PATCH 10/23] Make FieldElement limbs private to the curve25519-dalek crate. Limbs are no longer accessible outside of the curve25519-dalek crate. If you were relying on this behaviour, first you probably shouldn't be doing that, second please contact us so we can determine the best way forward for your use case. --- src/field_32bit.rs | 2 +- src/field_64bit.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/field_32bit.rs b/src/field_32bit.rs index ad85599..8a5cc12 100644 --- a/src/field_32bit.rs +++ b/src/field_32bit.rs @@ -56,7 +56,7 @@ use utils::{load3, load4}; /// faster. However, the `FieldElement64` implementation requires Rust's /// `u128`, which is not yet stable. #[derive(Copy, Clone)] -pub struct FieldElement32(pub [i32; 10]); +pub struct FieldElement32(pub (crate) [i32; 10]); impl Debug for FieldElement32 { fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { diff --git a/src/field_64bit.rs b/src/field_64bit.rs index a4da3a0..9eb8f6e 100644 --- a/src/field_64bit.rs +++ b/src/field_64bit.rs @@ -51,7 +51,7 @@ pub type Limb = u64; /// which gives even better performance. This implementation requires Rust's /// `u128`, which is not yet stable. #[derive(Copy, Clone)] -pub struct FieldElement64(pub [u64; 5]); +pub struct FieldElement64(pub (crate) [u64; 5]); impl Debug for FieldElement64 { fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { From c08591a7d81da1f5471e399f571621811351b723 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Wed, 4 Oct 2017 04:43:04 +0000 Subject: [PATCH 11/23] Improve documentation for Mongomery code. --- src/montgomery.rs | 39 ++++++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/src/montgomery.rs b/src/montgomery.rs index 1f4f52d..1020a24 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -35,7 +35,9 @@ use field::FieldElement; use edwards::{ExtendedPoint, CompressedEdwardsY}; use scalar::Scalar; -// XXX move these to a common "traits" or "group" module? —isis +// 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::slices_equal; @@ -54,8 +56,6 @@ use subtle::Mask; /// 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]); @@ -211,7 +211,7 @@ impl CompressedMontgomeryU { /// 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 +/// v-coordinate is superfluous for the purposes of scalar multiplication, we merely /// use `(U:W)`. #[derive(Copy, Clone, Debug)] #[allow(missing_docs)] @@ -224,9 +224,9 @@ pub struct MontgomeryPoint{ /// /// 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) . +///     ⎧ (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 @@ -325,11 +325,12 @@ impl MontgomeryPoint { /// 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? + /// The doubling case is degenerate, in that using this method to accomplish + /// point doubling is less efficient than using `differential_double()`. fn differential_add(&self, that: &MontgomeryPoint, difference: &MontgomeryPoint) -> MontgomeryPoint { - // debug_assert!(self.ct_eq(that) != 1); // The doubling case is degenerate + // 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()); @@ -342,13 +343,20 @@ impl MontgomeryPoint { } } - /// Differential doubling for single-coordinate Montgomery points. + /// Pseudo-doubling for single-coordinate Montgomery points. /// - /// DOCDOC + /// 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. + /// 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; @@ -370,8 +378,9 @@ impl MontgomeryPoint { /// Multiply this `MontgomeryPoint` by a `Scalar`. /// -/// DOCDOC -/// explain montgomery laddering +/// 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; From 49f8781681016b0ac626dcbe8e007e15509030d5 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 00:49:49 +0000 Subject: [PATCH 12/23] Revert "Use subtle branch with ConditionallySwappable trait for now." This reverts commit d86fcd4e17e8a6cdb3913d3a2fa82a91d8403801. --- Cargo.toml | 2 -- 1 file changed, 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index c1cd666..8c22506 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -35,8 +35,6 @@ version = "0.6" [dependencies.subtle] version = "^0.2" default-features = false -git = "https://github.com/isislovecruft/subtle" -branch = "feature/conditional-swap_r1" [dependencies.generic-array] # same version that digest depends on From cb98a23230b94ff96693b4589a9a00c2c042ccbd Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 00:51:42 +0000 Subject: [PATCH 13/23] Change Cargo.toml to use subtle-0.3.0. --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 8c22506..f8a879e 100644 --- a/Cargo.toml +++ b/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] From d39cb275c5a8d8c315d916760753dc9604a015fd Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 00:56:40 +0000 Subject: [PATCH 14/23] Remove DecafPoint.to_edwards() method. --- src/decaf.rs | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/decaf.rs b/src/decaf.rs index 23511dd..9bc2c8a 100644 --- a/src/decaf.rs +++ b/src/decaf.rs @@ -189,16 +189,11 @@ 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); impl DecafPoint { - /// Convert this `DecafPoint` to its underlying `ExtendedPoint`. - pub fn to_edwards(&self) -> ExtendedPoint { - self.0 - } - /// Compress in Decaf format. pub fn compress(&self) -> CompressedDecaf { // Q: Do we want to encode twisted or untwisted? @@ -757,7 +752,7 @@ mod test { fn decaf_decompress_id() { let compressed_id = CompressedDecaf::identity(); let id = compressed_id.decompress().unwrap(); - assert_eq!(id.to_edwards().compress(), CompressedEdwardsY::identity()); + assert_eq!(id.0.compress(), CompressedEdwardsY::identity()); } #[test] @@ -769,7 +764,7 @@ mod test { #[test] fn decaf_basepoint_roundtrip() { let bp_compressed_decaf = constants::DECAF_ED25519_BASEPOINT_POINT.compress(); - let bp_recaf = bp_compressed_decaf.decompress().unwrap().to_edwards(); + let bp_recaf = bp_compressed_decaf.decompress().unwrap().0; // 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 @@ -843,7 +838,7 @@ mod test { for _ in 0..100 { let P = DecafPoint::random(&mut rng); // Check that P is on the curve - assert!(P.to_edwards().is_valid()); + assert!(P.0.is_valid()); // Check that P is in the image of the decaf map P.compress(); } From 7e4fd5677cde54837e2583e3662a1bccfbc49d97 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 01:26:15 +0000 Subject: [PATCH 15/23] Add tests and benchmark for MontgomeryPoint.ct_eq(). --- src/montgomery.rs | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/montgomery.rs b/src/montgomery.rs index 1020a24..6c32d81 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -425,6 +425,7 @@ impl<'a, 'b> Mul<&'b MontgomeryPoint> for &'a Scalar { #[cfg(test)] mod test { + use constants::ED25519_BASEPOINT_TABLE; use constants::BASE_COMPRESSED_MONTGOMERY; use edwards::Identity; use super::*; @@ -481,6 +482,26 @@ mod test { 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(); @@ -532,6 +553,17 @@ mod bench { 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()); From 5e6e6c3fa824ea0a9b7678113e25325324df0b04 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 01:32:44 +0000 Subject: [PATCH 16/23] Eliminate extra inversions in MontgomeryPoint.ct_eq(). --- src/montgomery.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/montgomery.rs b/src/montgomery.rs index 6c32d81..ba3d556 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -257,8 +257,9 @@ impl Identity for MontgomeryPoint { /// `1` if the points are equal, and `0` otherwise. impl Equal for MontgomeryPoint { fn ct_eq(&self, that: &MontgomeryPoint) -> u8 { - slices_equal(self.compress().as_bytes(), - that.compress().as_bytes()) + // (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)) } } From 7b378ada6bfacb237cfae2d24a0df4046a171959 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 01:46:42 +0000 Subject: [PATCH 17/23] Rephrase doc note on exceptional projective Montgomery points. --- src/montgomery.rs | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/montgomery.rs b/src/montgomery.rs index ba3d556..c42fb42 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -267,9 +267,10 @@ impl Equal for MontgomeryPoint { /// /// # 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). +/// 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 /// From ca5b58c2b75d48656d142bb8c686369448f33945 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 01:54:26 +0000 Subject: [PATCH 18/23] Clarify doc note on degenerate cases for differential addition. --- src/montgomery.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/montgomery.rs b/src/montgomery.rs index c42fb42..bb6a9ad 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -327,8 +327,8 @@ impl MontgomeryPoint { /// 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 using this method to accomplish - /// point doubling is less efficient than using `differential_double()`. + /// 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 From 9da24d8afaf0b9443664604d3272dd336cf99c88 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 02:15:28 +0000 Subject: [PATCH 19/23] Add test for Montgomery ladder with a scalar with high bit set. --- src/montgomery.rs | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/montgomery.rs b/src/montgomery.rs index bb6a9ad..37f447a 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -545,6 +545,19 @@ mod test { 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"))] From 29f9090411e2e9d3d06e8a916a38ddcba67e7b2d Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 02:23:34 +0000 Subject: [PATCH 20/23] Fix two typos in docstrings for constants. --- src/constants_32bit.rs | 2 +- src/constants_64bit.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/constants_32bit.rs b/src/constants_32bit.rs index 02b629b..bf46fe5 100644 --- a/src/constants_32bit.rs +++ b/src/constants_32bit.rs @@ -75,7 +75,7 @@ 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.) +/// `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) diff --git a/src/constants_64bit.rs b/src/constants_64bit.rs index 8beee81..ba7b2d9 100644 --- a/src/constants_64bit.rs +++ b/src/constants_64bit.rs @@ -54,7 +54,7 @@ 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.) +/// `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) From 4965238b5a4e3a16ec0009a52f5af0e0943e72dd Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 02:23:55 +0000 Subject: [PATCH 21/23] Removed now unused subtle import from montgomery module. --- src/montgomery.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/montgomery.rs b/src/montgomery.rs index 37f447a..5c34114 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -40,7 +40,6 @@ use scalar::Scalar; // XXX Rust. —isis use edwards::{Identity, ValidityCheck}; -use subtle::slices_equal; use subtle::ConditionallyAssignable; use subtle::ConditionallySwappable; use subtle::Equal; From d39e47ff116cf1c6c47c84b7354225772e5f2d97 Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 02:27:27 +0000 Subject: [PATCH 22/23] Remove comment on non-canonical encodings in CompressedMontgomeryU.decompress(). --- src/montgomery.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/montgomery.rs b/src/montgomery.rs index 5c34114..9ba71ed 100644 --- a/src/montgomery.rs +++ b/src/montgomery.rs @@ -120,7 +120,6 @@ impl CompressedMontgomeryU { /// A projective `MontgomeryPoint` corresponding to this compressed point. pub fn decompress(&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(), } From bf25b2767b6a4afbbd3d19595996e1b9069ff0bc Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Thu, 5 Oct 2017 06:29:39 +0000 Subject: [PATCH 23/23] Bump curve25519-dalek version to 0.12.0. --- Cargo.toml | 2 +- README.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index f8a879e..4b477e2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "curve25519-dalek" -version = "0.11.0" +version = "0.12.0" authors = ["Isis Lovecruft ", "Henry de Valence "] readme = "README.md" diff --git a/README.md b/README.md index 3df4040..77420cd 100644 --- a/README.md +++ b/README.md @@ -44,7 +44,7 @@ Extensive documentation is available [here](https://docs.rs/curve25519-dalek). To install, add the following to the dependencies section of your project's `Cargo.toml`: - curve25519-dalek = "^0.11" + curve25519-dalek = "^0.12" Then, in your library or executable source, add: