diff --git a/Cargo.toml b/Cargo.toml index 8c22506..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" @@ -33,7 +33,7 @@ version = "0.3" version = "0.6" [dependencies.subtle] -version = "^0.2" +version = "^0.3" default-features = false [dependencies.generic-array] 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: 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..004add7 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; @@ -75,6 +76,9 @@ pub const HALF: FieldElement32 = FieldElement32([ pub const A: FieldElement32 = FieldElement32([ 486662, 0, 0, 0, 0, 0, 0, 0, 0, 0, ]); +/// `APLUS2_OVER_FOUR` is (A+2)/4. (This is used internally within the Montgomery ladder.) +pub const APLUS2_OVER_FOUR: FieldElement32 = FieldElement32([121666, 0, 0, 0, 0, 0, 0, 0, 0, 0]); + /// `SQRT_MINUS_A` is sqrt(-486662) // XXX I think that this was used in Adam's code for his elligator // implementation, but that should maybe be using sqrt(-486664) @@ -93,6 +97,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..e650757 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; @@ -54,6 +55,9 @@ pub const HALF: FieldElement64 = FieldElement64([2251799813685239, 2251799813685 /// In Montgomery form y² = x³+Ax²+x, Curve25519 has A=486662. pub const A: FieldElement64 = FieldElement64([486662, 0, 0, 0, 0]); +/// `APLUS2_OVER_FOUR` is (A+2)/4. (This is used internally within the Montgomery ladder.) +pub const APLUS2_OVER_FOUR: FieldElement64 = FieldElement64([121666, 0, 0, 0, 0]); + /// `SQRT_MINUS_A` is sqrt(-486662) // XXX I think that this was used in Adam's code for his elligator // implementation, but that should maybe be using sqrt(-486664) @@ -66,6 +70,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/decaf.rs b/src/decaf.rs index f8eb0c4..9bc2c8a 100644 --- a/src/decaf.rs +++ b/src/decaf.rs @@ -189,7 +189,7 @@ impl<'de> Deserialize<'de> for DecafPoint { /// A point in a prime-order group. /// -/// XXX think about how this API should work +// XXX think about how this API should work #[derive(Copy, Clone)] pub struct DecafPoint(pub ExtendedPoint); @@ -752,7 +752,7 @@ mod test { fn decaf_decompress_id() { let compressed_id = CompressedDecaf::identity(); let id = compressed_id.decompress().unwrap(); - assert_eq!(id.0.compress_edwards(), CompressedEdwardsY::identity()); + assert_eq!(id.0.compress(), CompressedEdwardsY::identity()); } #[test] @@ -768,7 +768,7 @@ mod test { // Check that bp_recaf differs from bp by a point of order 4 let diff = &constants::ED25519_BASEPOINT_POINT - &bp_recaf; let diff4 = diff.mult_by_pow_2(4); // XXX this is wrong - assert_eq!(diff4.compress_edwards(), CompressedEdwardsY::identity()); + assert_eq!(diff4.compress(), CompressedEdwardsY::identity()); } #[test] diff --git a/src/edwards.rs b/src/edwards.rs index 2456954..0b8b333 100644 --- a/src/edwards.rs +++ b/src/edwards.rs @@ -89,7 +89,7 @@ use core::ops::Index; use constants; use field::FieldElement; use scalar::Scalar; -use montgomery::CompressedMontgomeryU; +use montgomery::MontgomeryPoint; use subtle::slices_equal; use subtle::bytes_equal; @@ -123,7 +123,6 @@ impl CompressedEdwardsY { } /// Copy this `CompressedEdwardsY` to an array of bytes. - /// XXX is this useful? pub fn to_bytes(&self) -> [u8; 32] { self.0 } @@ -169,7 +168,7 @@ impl Serialize for ExtendedPoint { fn serialize(&self, serializer: S) -> Result where S: Serializer { - serializer.serialize_bytes(self.compress_edwards().as_bytes()) + serializer.serialize_bytes(self.compress().as_bytes()) } } @@ -393,8 +392,8 @@ impl ConditionallyAssignable for ExtendedPoint { impl Equal for ExtendedPoint { fn ct_eq(&self, other: &ExtendedPoint) -> u8 { - slices_equal(self.compress_edwards().as_bytes(), - other.compress_edwards().as_bytes()) + slices_equal(self.compress().as_bytes(), + other.compress().as_bytes()) } } @@ -436,7 +435,7 @@ impl ProjectivePoint { } /// Convert this point to a `CompressedEdwardsY` - pub fn compress_edwards(&self) -> CompressedEdwardsY { + pub fn compress(&self) -> CompressedEdwardsY { let recip = self.Z.invert(); let x = &self.X * &recip; let y = &self.Y * &recip; @@ -447,32 +446,67 @@ impl ProjectivePoint { CompressedEdwardsY(s) } - /// Convert this point to a `CompressedMontgomeryU`. - /// 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(CompressedMontgomeryU)` otherwise. + /// Taking the Montgomery curve equation in affine coordinates: /// - pub fn compress_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(CompressedMontgomeryU(u.to_bytes())) - } else { - None + ///     E_(A,B) = Bv² = u³ + Au² + u   (1) + /// + /// and given its relations to the coordinates of the Edwards model: + /// + ///     u = (1+y)/(1-y)        (2) + ///     v = (λu)/(x) + /// + /// 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, } } } @@ -515,20 +549,15 @@ impl ExtendedPoint { } } - /// Compress this point to `CompressedEdwardsY` format. - pub fn compress_edwards(&self) -> CompressedEdwardsY { - self.to_projective().compress_edwards() + /// Convert this point to its equivalent on the Montgomery form of the + /// curve. + pub fn to_montgomery(&self) -> MontgomeryPoint { + self.to_projective().to_montgomery() } - /// 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() + /// Compress this point to `CompressedEdwardsY` format. + pub fn compress(&self) -> CompressedEdwardsY { + self.to_projective().compress() } } @@ -1305,7 +1334,7 @@ mod test { assert!(bp.is_valid()); // Check that decompression actually gives the correct X coordinate assert_eq!(base_X, bp.X); - assert_eq!(bp.compress_edwards(), constants::BASE_CMPRSSD); + assert_eq!(bp.compress(), constants::BASE_CMPRSSD); } /// Test sign handling in decompression @@ -1328,7 +1357,7 @@ mod test { #[test] fn basepoint_mult_one_vs_basepoint() { let bp = &constants::ED25519_BASEPOINT_TABLE * &Scalar::one(); - let compressed = bp.compress_edwards(); + let compressed = bp.compress(); assert_eq!(compressed, constants::BASE_CMPRSSD); } @@ -1336,7 +1365,7 @@ mod test { #[test] fn basepoint_table_basepoint_function_correct() { let bp = constants::ED25519_BASEPOINT_TABLE.basepoint(); - assert_eq!(bp.compress_edwards(), constants::BASE_CMPRSSD); + assert_eq!(bp.compress(), constants::BASE_CMPRSSD); } /// Test `impl Add for ExtendedPoint` @@ -1345,7 +1374,7 @@ mod test { fn basepoint_plus_basepoint_vs_basepoint2() { let bp = constants::ED25519_BASEPOINT_POINT; let bp_added = &bp + &bp; - assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD); + assert_eq!(bp_added.compress(), BASE2_CMPRSSD); } /// Test `impl Add for ExtendedPoint` @@ -1354,7 +1383,7 @@ mod test { fn basepoint_plus_basepoint_projective_niels_vs_basepoint2() { let bp = constants::ED25519_BASEPOINT_POINT; let bp_added = (&bp + &bp.to_projective_niels()).to_extended(); - assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD); + assert_eq!(bp_added.compress(), BASE2_CMPRSSD); } /// Test `impl Add for ExtendedPoint` @@ -1364,7 +1393,7 @@ mod test { let bp = constants::ED25519_BASEPOINT_POINT; let bp_affine_niels = bp.to_affine_niels(); let bp_added = (&bp + &bp_affine_niels).to_extended(); - assert_eq!(bp_added.compress_edwards(), BASE2_CMPRSSD); + assert_eq!(bp_added.compress(), BASE2_CMPRSSD); } /// Check that equality of `ExtendedPoints` handles projective @@ -1389,15 +1418,15 @@ mod test { let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR; let aB_affine_niels = aB.to_affine_niels(); let also_aB = (&ExtendedPoint::identity() + &aB_affine_niels).to_extended(); - assert_eq!( aB.compress_edwards(), - 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 @@ -1415,20 +1444,20 @@ mod test { let table = EdwardsBasepointTable::create(&constants::ED25519_BASEPOINT_POINT); let aB_1 = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR; let aB_2 = &table * &A_SCALAR; - assert_eq!(aB_1.compress_edwards(), aB_2.compress_edwards()); + assert_eq!(aB_1.compress(), aB_2.compress()); } /// Test scalar_mult versus a known scalar multiple from ed25519.py #[test] fn scalar_mult_vs_ed25519py() { let aB = &constants::ED25519_BASEPOINT_POINT * &A_SCALAR; - assert_eq!(aB.compress_edwards(), A_TIMES_BASEPOINT); + assert_eq!(aB.compress(), A_TIMES_BASEPOINT); } /// Test basepoint.double() versus the 2*basepoint constant. #[test] fn basepoint_double_vs_basepoint2() { - assert_eq!(constants::ED25519_BASEPOINT_POINT.double().compress_edwards(), + assert_eq!(constants::ED25519_BASEPOINT_POINT.double().compress(), BASE2_CMPRSSD); } @@ -1437,14 +1466,14 @@ mod test { fn basepoint_mult_two_vs_basepoint2() { let mut two_bytes = [0u8; 32]; two_bytes[0] = 2; let bp2 = &constants::ED25519_BASEPOINT_TABLE * &Scalar(two_bytes); - assert_eq!(bp2.compress_edwards(), BASE2_CMPRSSD); + assert_eq!(bp2.compress(), BASE2_CMPRSSD); } /// Check that converting to projective and then back to extended round-trips. #[test] fn basepoint_projective_extended_round_trip() { assert_eq!(constants::ED25519_BASEPOINT_POINT - .to_projective().to_extended().compress_edwards(), + .to_projective().to_extended().compress(), constants::BASE_CMPRSSD); } @@ -1452,7 +1481,7 @@ mod test { #[test] fn basepoint16_vs_mult_by_pow_2_4() { let bp16 = constants::ED25519_BASEPOINT_POINT.mult_by_pow_2(4); - assert_eq!(bp16.compress_edwards(), BASE16_CMPRSSD); + assert_eq!(bp16.compress(), BASE16_CMPRSSD); } /// Test that the conditional assignment trait works for AffineNielsPoints. @@ -1480,7 +1509,7 @@ mod test { #[test] fn compressed_identity() { - assert_eq!(ExtendedPoint::identity().compress_edwards(), + assert_eq!(ExtendedPoint::identity().compress(), CompressedEdwardsY::identity()); } @@ -1518,7 +1547,7 @@ mod test { let P1 = &G * &s; let P2 = &s * &G; - assert!(P1.compress_edwards().to_bytes() == P2.compress_edwards().to_bytes()); + assert!(P1.compress().to_bytes() == P2.compress().to_bytes()); } #[test] @@ -1542,7 +1571,7 @@ mod test { fn double_scalar_mult_basepoint_vs_ed25519py() { let A = A_TIMES_BASEPOINT.decompress().unwrap(); let result = vartime::double_scalar_mult_basepoint(&A_SCALAR, &A, &B_SCALAR); - assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT); + assert_eq!(result.compress(), DOUBLE_SCALAR_MULT_RESULT); } #[test] @@ -1552,7 +1581,7 @@ mod test { &[A_SCALAR, B_SCALAR], &[A, constants::ED25519_BASEPOINT_POINT] ); - assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT); + assert_eq!(result.compress(), DOUBLE_SCALAR_MULT_RESULT); } #[test] @@ -1567,7 +1596,7 @@ mod test { &[A, constants::ED25519_BASEPOINT_POINT] ); - assert_eq!(result_vartime.compress_edwards(), result_consttime.compress_edwards()); + assert_eq!(result_vartime.compress(), result_consttime.compress()); } } @@ -1579,7 +1608,7 @@ mod test { fn serde_cbor_basepoint_roundtrip() { let output = serde_cbor::to_vec(&constants::ED25519_BASEPOINT_POINT).unwrap(); let parsed: ExtendedPoint = serde_cbor::from_slice(&output).unwrap(); - assert_eq!(parsed.compress_edwards(), constants::BASE_CMPRSSD); + assert_eq!(parsed.compress(), constants::BASE_CMPRSSD); } #[test] @@ -1615,7 +1644,7 @@ mod bench { #[bench] fn edwards_compress(b: &mut Bencher) { let B = &constants::ED25519_BASEPOINT_POINT; - b.iter(|| B.compress_edwards()); + b.iter(|| B.compress()); } #[bench] 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/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 { 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/montgomery.rs b/src/montgomery.rs index 29e1b20..9ba71ed 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,22 @@ // affine and projective cakes and eat both of them too. #![allow(non_snake_case)] +use core::ops::{Mul, MulAssign}; use constants; use field::FieldElement; use edwards::{ExtendedPoint, CompressedEdwardsY}; +use scalar::Scalar; + +// XXX Move these to a common "group" module? At the same time, we should +// XXX probably make a `trait Group` once const generics are implemented in +// XXX Rust. —isis +use edwards::{Identity, ValidityCheck}; use subtle::ConditionallyAssignable; +use subtle::ConditionallySwappable; +use subtle::Equal; +use subtle::Mask; /// In "Montgomery u" format, as used in X25519, a point `(u,v)` on /// the Montgomery curve @@ -33,13 +55,16 @@ use subtle::ConditionallyAssignable; /// coordinates. For Montgomery curves, it is possible to compute the /// `u`-coordinate of `n(u,v)` just from `n` and `u`, so it is not /// necessary to use `v` for a Diffie-Hellman key exchange. -/// -/// XXX add note on monty, twist security, edwards impl of x25519, rfc7748 #[derive(Copy, Clone, Debug, PartialEq, Eq)] pub struct CompressedMontgomeryU(pub [u8; 32]); impl CompressedMontgomeryU { /// View this `CompressedMontgomeryU` as an array of bytes. + pub fn as_bytes<'a>(&'a self) -> &'a [u8; 32] { + &self.0 + } + + /// Convert this `CompressedMontgomeryU` to an array of bytes. pub fn to_bytes(&self) -> [u8; 32] { self.0 } @@ -65,7 +90,7 @@ impl CompressedMontgomeryU { /// * `v` is not square. // // XXX any other exceptional points for the birational map? - pub fn decompress(&self) -> Option { + 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 +109,22 @@ impl CompressedMontgomeryU { CompressedEdwardsY(y.to_bytes()).decompress() } + /// Decompress this `CompressedMontgomeryU` to a `MontgomeryPoint`. + /// + /// Going from affine to projective coordinates, we have: + /// + ///     u → U/W + /// + /// # Returns + /// + /// A projective `MontgomeryPoint` corresponding to this compressed point. + pub fn decompress(&self) -> MontgomeryPoint { + MontgomeryPoint{ + U: FieldElement::from_bytes(&self.0), + W: FieldElement::one(), + } + } + /// Given a Montgomery `u` coordinate, compute an Edwards `y` via /// `y = (u-1)/(u+1)`. /// @@ -150,34 +191,259 @@ impl CompressedMontgomeryU { } } +/// A point on the Montgomery form of the curve, in projective 𝗣^2 coordinates. +/// +/// The transition between affine and projective is given by +/// +///     u → U/W +///     v → V/W +/// +/// thus the Montgomery curve equation +/// +///     E_(A,B) : Bv² = u(u² + Au + 1) +/// +/// becomes +/// +///     E_(A,B) : BV²W = U(U² + AUW + W²) ⊆ 𝗣^2 +/// +/// Here, again, to differentiate from points in the twisted Edwards model, we +/// call the point `(x,y)` in affine coordinates `(u,v)` and similarly in projective +/// space we use `(U:V:W)`. However, since (as per Montgomery's original work) the +/// v-coordinate is superfluous for the purposes of scalar multiplication, we merely +/// use `(U:W)`. +#[derive(Copy, Clone, Debug)] +#[allow(missing_docs)] +pub struct MontgomeryPoint{ + pub U: FieldElement, + pub W: FieldElement, +} + +/// The identity point is a unique point (the only where `W = 0`) on the curve. +/// +/// In projective coordinates, the quotient map `x : E (A,B) → E/<⦵> = 𝗣¹` is +/// +///     ⎧ (x_P:1) if P = (x_P:y_P:1) , +///     x : P ↦ ⎨ +///     ⎩ (1:0) if P = O = (0:1:0) . +/// +/// We emphasize that the formula `x((U: V : W)) = (U : W)` only holds on the +/// open subset of `E_(A,B)` where `W ≠ 0`; it does not extend to the point +/// `O = (0:1:0)` at infinity, because `(0:0)` is not a projective point. +/// +/// # Returns +/// +/// The (exceptional) point at infinity in the Montgomery model. +impl Identity for MontgomeryPoint { + fn identity() -> MontgomeryPoint { + MontgomeryPoint { + U: FieldElement::one(), + W: FieldElement::zero(), + } + } +} + +/// Determine if two `MontgomeryPoint`s are equal, in constant time. +/// +/// # Note +/// +/// Because a compressed point on the Montgomery form of the curve doesn't +/// include the sign bit, there's two points here (if translated from the +/// Edwards form) which will equate. +/// +/// # Returns +/// +/// `1` if the points are equal, and `0` otherwise. +impl Equal for MontgomeryPoint { + fn ct_eq(&self, that: &MontgomeryPoint) -> u8 { + // (U_P:W_P) = (U_Q:W_Q) iff U_P * W_Q == U_Q * W_P, + // since U_P/W_P == U_Q/W_Q. + (&self.U * &that.W).ct_eq(&(&self.W * &that.U)) + } +} + +/// Determine if this `MontgomeryPoint` is valid. +/// +/// # Note +/// +/// All projective points, except for `(X:W) = (0:0)`, are valid, since the +/// projective model is linear through the origin and is comprised by all `X` in +/// ℤ/(2²⁵⁵-19), thus `(0:0)` is the only element in Fₚ² which is not a +/// projective point. +/// +/// # Returns +/// +/// `true` if it is valid, and `false` otherwise. +impl ValidityCheck for MontgomeryPoint { + fn is_valid(&self) -> bool { + let zero = FieldElement::zero(); + + if (self.U.ct_eq(&zero) & self.W.ct_eq(&zero)) == 1 { + return true; + } + false + } +} + +/// Conditionally assign another `MontgomeryPoint` to this point, in constant time. +/// +/// If `choice == 1`, assign `that` to `self`. Otherwise, leave `self` +/// unchanged. +impl ConditionallyAssignable for MontgomeryPoint { + fn conditional_assign(&mut self, that: &MontgomeryPoint, choice: Mask) { + self.U.conditional_assign(&that.U, choice); + self.W.conditional_assign(&that.W, choice); + } +} + +impl MontgomeryPoint { + /// Compress this point to only its u-coordinate (note: affine). + /// + /// # Returns + /// + /// A `CompressedMontgomeryU`. + pub fn compress(&self) -> CompressedMontgomeryU { + let u_affine: FieldElement = &self.U * &self.W.invert(); + + CompressedMontgomeryU(u_affine.to_bytes()) + } + + /// Differential addition for single-coordinate Montgomery points. + /// + /// Montgomery coordinates in projective 𝗣¹ space are odd in that 𝗣¹ + /// inherits none of the group structure from E_(A,B). Hence, the mapping + /// of the group operation, `⊕`, is undefined for the pair `(x(P), x(Q))`; + /// that is, given `x(P)` and `x(Q)`, we cannot derive `x(P ⊕ Q)`. This is + /// due to the fact that, in Montgomery coordinates, `x(P)` determines `P` + /// only up to a sign, and thus we cannot differentiate `x(P ⊕ Q)` from + /// `x(P ⊖ Q)`. However, via differential addition, any three of the values + /// `{x(P), x(Q), x(P ⊕ Q), x(P ⊖ Q)}` determines the forth, so we can + /// define *pseudo-addition* for a singular coordinate. + /// + /// # Warning + /// + /// If the `difference` is the identity point, or a two torsion point, the + /// results of this method are not correct, but instead result in `(0:0)` + /// (an invalid projective point in the Montgomery model). + /// + /// The doubling case is degenerate, in that `P ⦵ Q ∉ {O,T}`, where `T` is + /// the two torsion point. + fn differential_add(&self, that: &MontgomeryPoint, + difference: &MontgomeryPoint) -> MontgomeryPoint { + // XXX Do we want these debug assertions? We would need to implement + // XXX is_two_torsion_point(). —isis + // debug_assert!(!difference.is_identity()); // P ⦵ Q ∉ {O,T} + // debug_assert!(!difference.is_two_torsion_point()); + + let v1: FieldElement = &(&self.U + &self.W) * &(&that.U - &that.W); + let v2: FieldElement = &(&self.U - &self.W) * &(&that.U + &that.W); + + MontgomeryPoint { + U: &difference.W * &(&v1 + &v2).square(), // does reduction on square() + W: &difference.U * &(&v1 - &v2).square(), // does reduction on square() + } + } + + /// Pseudo-doubling for single-coordinate Montgomery points. + /// + /// Given a Montgomery U-coordinate of a point `P`, compute the + /// U-coordinate given by + /// + ///     differential_double: x(P) ⟼ x([2]P) + /// + /// # Returns + /// + /// A Montgomery point equal to doubling this one. + /// + // XXX It seems possible that combining the differential_add() and + // XXX differential_double() methods would save a non-trivial amount of + // XXX computation in the ladder. —isis + fn differential_double(&self) -> MontgomeryPoint { + let mut v1: FieldElement; + let v2: FieldElement; + let v3: FieldElement; + + v1 = (&self.U + &self.W).square(); + v2 = (&self.U - &self.W).square(); + + let U: FieldElement = &v1 * &v2; + + v1 -= &v2; + v3 = &(&constants::APLUS2_OVER_FOUR * &v1) + &v2; + + let W: FieldElement = &v1 * &v3; + + MontgomeryPoint{ U: U, W: W } + } +} + +/// Multiply this `MontgomeryPoint` by a `Scalar`. +/// +/// The reader is refered to §5.3 of ["Montgomery Curves and Their Arithmetic" +/// by Craig Costello and Benjamin Smith](https://eprint.iacr.org/2017/212.pdf) +/// for an overview of side-channel-free Montgomery laddering algorithms. +impl<'a, 'b> Mul<&'b Scalar> for &'a MontgomeryPoint { + type Output = MontgomeryPoint; + + fn mul(self, scalar: &'b Scalar) -> MontgomeryPoint { + let mut x0: MontgomeryPoint = MontgomeryPoint::identity(); + let mut x1: MontgomeryPoint = *self; + + let bits: [i8; 256] = scalar.bits(); + + for i in (0..255).rev() { + let mask: u8 = (bits[i+1] ^ bits[i]) as u8; + + debug_assert!(mask == 0 || mask == 1); + + x0.conditional_swap(&mut x1, mask); + x1 = x0.differential_add(&x1, &self); + x0 = x0.differential_double(); + } + x0.conditional_swap(&mut x1, bits[0] as u8); + x0 + } +} + +impl<'b> MulAssign<&'b Scalar> for MontgomeryPoint { + fn mul_assign(&mut self, scalar: &'b Scalar) { + let result = (self as &MontgomeryPoint) * scalar; + *self = result; + } +} + +impl<'a, 'b> Mul<&'b MontgomeryPoint> for &'a Scalar { + type Output = MontgomeryPoint; + + fn mul(self, point: &'b MontgomeryPoint) -> MontgomeryPoint { + point * &self + } +} + // ------------------------------------------------------------------------ // Tests // ------------------------------------------------------------------------ #[cfg(test)] mod test { + use constants::ED25519_BASEPOINT_TABLE; + use constants::BASE_COMPRESSED_MONTGOMERY; use edwards::Identity; use super::*; - /// The X25519 basepoint, in compressed Montgomery form. - 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); + 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_CMPRSSD_MONTY.decompress().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. @@ -189,14 +455,146 @@ mod test { let minus_one = FieldElement::minus_one(); let minus_one_bytes = minus_one.to_bytes(); let div_by_zero_u = CompressedMontgomeryU(minus_one_bytes); - assert!(div_by_zero_u.decompress().is_none()); + assert!(div_by_zero_u.decompress_edwards().is_none()); } - /// Montgomery compression of the identity point should - /// 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() { + assert_eq!(BASE_COMPRESSED_MONTGOMERY.decompress().compress(), + BASE_COMPRESSED_MONTGOMERY); + + } + + #[test] + fn differential_double_matches_double() { + let p: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double(); + let q: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress().differential_double(); + + assert_eq!(p.to_montgomery().compress(), q.compress()); + } + + #[test] + fn montgomery_ct_eq_ne() { + let mut csprng: OsRng = OsRng::new().unwrap(); + let s1: Scalar = Scalar::random(&mut csprng); + let s2: Scalar = Scalar::random(&mut csprng); + let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery(); + let p2: MontgomeryPoint = (&s2 * &ED25519_BASEPOINT_TABLE).to_montgomery(); + + assert_eq!(p1.ct_eq(&p2), 0); + } + + #[test] + fn montgomery_ct_eq_eq() { + let mut csprng: OsRng = OsRng::new().unwrap(); + let s1: Scalar = Scalar::random(&mut csprng); + let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery(); + + assert_eq!(p1.ct_eq(&p1), 1); + } + + #[test] + fn differential_add_matches_edwards_model() { + let mut csprng: OsRng = OsRng::new().unwrap(); + + let s1: Scalar = Scalar::random(&mut csprng); + let s2: Scalar = Scalar::random(&mut csprng); + let p1: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s1; + let p2: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s2; + let diff: ExtendedPoint = &p1 - &p2; + + let p1m: MontgomeryPoint = p1.to_montgomery(); + let p2m: MontgomeryPoint = p2.to_montgomery(); + let diffm: MontgomeryPoint = diff.to_montgomery(); + + let result = p1m.differential_add(&p2m, &diffm); + + assert_eq!(result.compress(), (&p1 + &p2).to_montgomery().compress()); + } + + #[test] + fn ladder_matches_scalarmult() { + let mut csprng: OsRng = OsRng::new().unwrap(); + + let s: Scalar = Scalar::random(&mut csprng); + let p_edwards: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s; + let p_montgomery: MontgomeryPoint = p_edwards.to_montgomery(); + + let expected = &s * &p_edwards; + let result = &s * &p_montgomery; + + assert_eq!(result.compress(), expected.to_montgomery().compress()) + } + + #[test] + fn ladder_basepoint_times_two_matches_double() { + let two: Scalar = Scalar::from_u64(2u64); + let result: MontgomeryPoint = &BASE_COMPRESSED_MONTGOMERY.decompress() * &two; + let expected: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double(); + + assert_eq!(result.compress(), expected.to_montgomery().compress()); + } + + #[test] + #[should_panic(expected = "assertion failed: self[31] <= 127")] + fn ladder_matches_scalarmult_with_scalar_high_bit_set() { + let mut s: Scalar = Scalar::one(); + + s[31] = 255; + + let result: MontgomeryPoint = &BASE_COMPRESSED_MONTGOMERY.decompress() * &s; + let expected: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s; + + assert_eq!(result.compress(), expected.to_montgomery().compress()) + } +} + +#[cfg(all(test, feature = "bench"))] +mod bench { + use rand::OsRng; + use constants::ED25519_BASEPOINT_TABLE; + use constants::BASE_COMPRESSED_MONTGOMERY; + use test::Bencher; + use super::*; + + #[bench] + fn montgomery_ct_eq(b: &mut Bencher) { + let mut csprng: OsRng = OsRng::new().unwrap(); + let s1: Scalar = Scalar::random(&mut csprng); + let s2: Scalar = Scalar::random(&mut csprng); + let p1: MontgomeryPoint = (&s1 * &ED25519_BASEPOINT_TABLE).to_montgomery(); + let p2: MontgomeryPoint = (&s2 * &ED25519_BASEPOINT_TABLE).to_montgomery(); + + b.iter(| | p1.ct_eq(&p2)) + } + + #[bench] + fn montgomery_decompress(b: &mut Bencher) { + b.iter(| | BASE_COMPRESSED_MONTGOMERY.decompress()); + } + + #[bench] + fn montgomery_compress(b: &mut Bencher) { + let p: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress(); + + b.iter(| | p.compress()); + } + + #[bench] + fn montgomery_ladder(b: &mut Bencher) { + let mut csprng: OsRng = OsRng::new().unwrap(); + let s: Scalar = Scalar::random(&mut csprng); + let p: MontgomeryPoint = (&Scalar::random(&mut csprng) * &ED25519_BASEPOINT_TABLE).to_montgomery(); + + b.iter(| | &s * &p); } } diff --git a/src/scalar.rs b/src/scalar.rs index 876061d..399d00d 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. //! @@ -20,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; @@ -41,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; @@ -106,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) } } @@ -232,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. /// @@ -382,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 /// @@ -440,293 +427,106 @@ 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. 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 = UnpackedScalar::multiply_add(&y, &y, &UnpackedScalar::zero()); - 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.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::montgomery_mul(&_1111, &_1); + + #[inline] + fn square_multiply(y: &mut UnpackedScalar, squarings: usize, x: &UnpackedScalar) { + for _ in 0..squarings { + *y = y.montgomery_square(); } - } - y - } - - /// 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; + *y = UnpackedScalar::montgomery_mul(y, x); } - 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); - } + 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, 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); + 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); - 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( @@ -734,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, @@ -838,6 +644,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]; @@ -862,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()); } @@ -894,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) { @@ -903,22 +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 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)); + } +}