From 6079b0269f221a63ccef7a60d8bcca6b25b579ca Mon Sep 17 00:00:00 2001 From: Isis Lovecruft Date: Mon, 30 Oct 2017 19:27:44 +0000 Subject: [PATCH] Revert "Revert "Merge remote-tracking branch 'floodyberry/optimized_scalar' into develop"" This reverts commit 90b69c13eed20c8edf734756b8cb24e3ff721523. Signed-off-by: Isis Lovecruft --- 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 bf46fe5..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; @@ -96,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 ba7b2d9..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; @@ -69,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/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)); + } +}