diff --git a/src/constants.rs b/src/constants.rs index 814f84e..4c159cd 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -37,8 +37,8 @@ use scalar::Scalar; #[cfg(feature="radix_51")] pub use constants_64bit::*; -#[cfg(not(feature="radix_51"))] -pub use constants_32bit::*; +//#[cfg(not(feature="radix_51"))] +//pub use constants_32bit::*; /// Basepoint has y = 4/5. /// diff --git a/src/field_32bit.rs b/src/field_32bit.rs index 25a352d..10c17e1 100644 --- a/src/field_32bit.rs +++ b/src/field_32bit.rs @@ -11,8 +11,8 @@ //! Field arithmetic for ℤ/(2²⁵⁵-19), using 32-bit arithmetic with //! 64-bit products. //! -//! Based on Adam Langley's curve25519-donna and (Golang) ed25519 -//! implementations. +//! This code was originally derived from Adam Langley's +//! curve25519-donna and (Golang) ed25519 implementations. //! //! This implementation is intended for platforms that can multiply //! 32-bit inputs to produce 64-bit outputs. @@ -30,14 +30,11 @@ use core::ops::Neg; use subtle::ConditionallyAssignable; -use utils::{load3, load4}; - /// A `FieldElement32` represents an element of the field GF(2^255 - 19). /// /// In the 32-bit implementation, a `FieldElement32` is represented in -/// radix 2^25.5 as ten `i32`s, so that an element t, entries -/// t[0],...,t[9], represents the integer t[0]+2^26 t[1]+2^51 -/// t[2]+2^77 t[3]+2^102 t[4]+...+2^230 t[9]. +/// radix 2^25.5 as ten `u32`s, so that an element t, entries +/// t[0],...,t[9], represents `sum(t[i]*2^ceil(i*51/2))`. /// /// The coefficients t[i] are allowed to grow between multiplications. /// @@ -56,7 +53,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 (crate) [i32; 10]); +pub struct FieldElement32(pub (crate) [u32; 10]); impl Debug for FieldElement32 { fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { @@ -83,9 +80,22 @@ impl<'a, 'b> Add<&'b FieldElement32> for &'a FieldElement32 { impl<'b> SubAssign<&'b FieldElement32> for FieldElement32 { fn sub_assign(&mut self, _rhs: &'b FieldElement32) { - for i in 0..10 { - self.0[i] -= _rhs.0[i]; - } + // See comment in FieldElement64::Sub + // + // Compute a - b as ((a + 2^4 * p) - b) to avoid underflow. + let b = &_rhs.0; + self.0 = FieldElement32::reduce([ + ((self.0[0] + (0x3ffffed << 4)) - b[0]) as u64, + ((self.0[1] + (0x1ffffff << 4)) - b[1]) as u64, + ((self.0[2] + (0x3ffffff << 4)) - b[2]) as u64, + ((self.0[3] + (0x1ffffff << 4)) - b[3]) as u64, + ((self.0[4] + (0x3ffffff << 4)) - b[4]) as u64, + ((self.0[5] + (0x1ffffff << 4)) - b[5]) as u64, + ((self.0[6] + (0x3ffffff << 4)) - b[6]) as u64, + ((self.0[7] + (0x1ffffff << 4)) - b[7]) as u64, + ((self.0[8] + (0x3ffffff << 4)) - b[8]) as u64, + ((self.0[9] + (0x1ffffff << 4)) - b[9]) as u64, + ]).0; } } @@ -108,86 +118,101 @@ impl<'b> MulAssign<&'b FieldElement32> for FieldElement32 { impl<'a, 'b> Mul<&'b FieldElement32> for &'a FieldElement32 { type Output = FieldElement32; fn mul(self, _rhs: &'b FieldElement32) -> FieldElement32 { - // Notes preserved from ed25519.go (presumably originally from ref10): - // - // Calculates h = f * g. Can overlap h with f or g. - // - // # Preconditions - // - // * |f[i]| bounded by 1.1*2^26, 1.1*2^25, 1.1*2^26, 1.1*2^25, etc. - // * |g[i]| bounded by 1.1*2^26, 1.1*2^25, 1.1*2^26, 1.1*2^25, etc. - // - // # Postconditions - // - // * |h| bounded by 1.1*2^25, 1.1*2^24, 1.1*2^25, 1.1*2^24, etc. - // - // ## Notes on implementation strategy - // - // * Using schoolbook multiplication. - // * Karatsuba would save a little in some cost models. - // - // * Most multiplications by 2 and 19 are 32-bit precomputations; - // cheaper than 64-bit postcomputations. - // - // * There is one remaining multiplication by 19 in the carry chain; - // one *19 precomputation can be merged into this, - // but the resulting data flow is considerably less clean. - // - // * There are 12 carries below. - // 10 of them are 2-way parallelizable and vectorizable. - // Can get away with 11 carries, but then data flow is much deeper. - // - // * With tighter constraints on inputs can squeeze carries into int32. - let f0 = self.0[0] as i64; - let f1 = self.0[1] as i64; - let f2 = self.0[2] as i64; - let f3 = self.0[3] as i64; - let f4 = self.0[4] as i64; - let f5 = self.0[5] as i64; - let f6 = self.0[6] as i64; - let f7 = self.0[7] as i64; - let f8 = self.0[8] as i64; - let f9 = self.0[9] as i64; + /// Helper function to multiply two 32-bit integers with 64 bits + /// of output. + #[inline(always)] + fn m(x: u32, y: u32) -> u64 { (x as u64) * (y as u64) } - let f1_2 = (2 * self.0[1]) as i64; - let f3_2 = (2 * self.0[3]) as i64; - let f5_2 = (2 * self.0[5]) as i64; - let f7_2 = (2 * self.0[7]) as i64; - let f9_2 = (2 * self.0[9]) as i64; + // Alias self, _rhs for more readable formulas + let x: &[u32;10] = &self.0; let y: &[u32;10] = &_rhs.0; - let g0 = _rhs.0[0] as i64; - let g1 = _rhs.0[1] as i64; - let g2 = _rhs.0[2] as i64; - let g3 = _rhs.0[3] as i64; - let g4 = _rhs.0[4] as i64; - let g5 = _rhs.0[5] as i64; - let g6 = _rhs.0[6] as i64; - let g7 = _rhs.0[7] as i64; - let g8 = _rhs.0[8] as i64; - let g9 = _rhs.0[9] as i64; + // We assume that the input limbs x[i], y[i] are bounded by: + // + // x[i], y[i] < 2^(26 + b) if i even + // x[i], y[i] < 2^(25 + b) if i odd + // + // where b is a (real) parameter representing the excess bits of + // the limbs. We track the bitsizes of all variables through + // the computation and solve at the end for the allowable + // headroom bitsize b (which determines how many additions we + // can perform between reductions or multiplications). - let g1_19 = (19 * _rhs.0[1]) as i64; /* 1.4*2^29 */ - let g2_19 = (19 * _rhs.0[2]) as i64; /* 1.4*2^30; still ok */ - let g3_19 = (19 * _rhs.0[3]) as i64; - let g4_19 = (19 * _rhs.0[4]) as i64; - let g5_19 = (19 * _rhs.0[5]) as i64; - let g6_19 = (19 * _rhs.0[6]) as i64; - let g7_19 = (19 * _rhs.0[7]) as i64; - let g8_19 = (19 * _rhs.0[8]) as i64; - let g9_19 = (19 * _rhs.0[9]) as i64; + let y1_19 = 19 * y[1]; // This fits in a u32 + let y2_19 = 19 * y[2]; // iff 26 + b + lg(19) < 32 + let y3_19 = 19 * y[3]; // if b < 32 - 26 - 4.248 = 1.752 + let y4_19 = 19 * y[4]; + let y5_19 = 19 * y[5]; // below, b<2.5: this is a bottleneck, + let y6_19 = 19 * y[6]; // could be avoided by promoting to + let y7_19 = 19 * y[7]; // u64 here instead of in m() + let y8_19 = 19 * y[8]; + let y9_19 = 19 * y[9]; - let h0 = f0*g0 + f1_2*g9_19 + f2*g8_19 + f3_2*g7_19 + f4*g6_19 + f5_2*g5_19 + f6*g4_19 + f7_2*g3_19 + f8*g2_19 + f9_2*g1_19; - let h1 = f0*g1 + f1*g0 + f2*g9_19 + f3*g8_19 + f4*g7_19 + f5*g6_19 + f6*g5_19 + f7*g4_19 + f8*g3_19 + f9*g2_19; - let h2 = f0*g2 + f1_2*g1 + f2*g0 + f3_2*g9_19 + f4*g8_19 + f5_2*g7_19 + f6*g6_19 + f7_2*g5_19 + f8*g4_19 + f9_2*g3_19; - let h3 = f0*g3 + f1*g2 + f2*g1 + f3*g0 + f4*g9_19 + f5*g8_19 + f6*g7_19 + f7*g6_19 + f8*g5_19 + f9*g4_19; - let h4 = f0*g4 + f1_2*g3 + f2*g2 + f3_2*g1 + f4*g0 + f5_2*g9_19 + f6*g8_19 + f7_2*g7_19 + f8*g6_19 + f9_2*g5_19; - let h5 = f0*g5 + f1*g4 + f2*g3 + f3*g2 + f4*g1 + f5*g0 + f6*g9_19 + f7*g8_19 + f8*g7_19 + f9*g6_19; - let h6 = f0*g6 + f1_2*g5 + f2*g4 + f3_2*g3 + f4*g2 + f5_2*g1 + f6*g0 + f7_2*g9_19 + f8*g8_19 + f9_2*g7_19; - let h7 = f0*g7 + f1*g6 + f2*g5 + f3*g4 + f4*g3 + f5*g2 + f6*g1 + f7*g0 + f8*g9_19 + f9*g8_19; - let h8 = f0*g8 + f1_2*g7 + f2*g6 + f3_2*g5 + f4*g4 + f5_2*g3 + f6*g2 + f7_2*g1 + f8*g0 + f9_2*g9_19; - let h9 = f0*g9 + f1*g8 + f2*g7 + f3*g6 + f4*g5 + f5*g4 + f6*g3 + f7*g2 + f8*g1 + f9*g0; + // What happens when we multiply x[i] with y[j] and place the + // result into the (i+j)-th limb? + // + // x[i] represents the value x[i]*2^ceil(i*51/2) + // y[j] represents the value y[j]*2^ceil(j*51/2) + // z[i+j] represents the value z[i+j]*2^ceil((i+j)*51/2) + // x[i]*y[j] represents the value x[i]*y[i]*2^(ceil(i*51/2)+ceil(j*51/2)) + // + // Since the radix is already accounted for, the result placed + // into the (i+j)-th limb should be + // + // x[i]*y[i]*2^(ceil(i*51/2)+ceil(j*51/2) - ceil((i+j)*51/2)). + // + // The value of ceil(i*51/2)+ceil(j*51/2) - ceil((i+j)*51/2) is + // 1 when both i and j are odd, and 0 otherwise. So we add + // + // x[i]*y[j] if either i or j is even + // 2*x[i]*y[j] if i and j are both odd + // + // by using precomputed multiples of x[i] for odd i: - FieldElement32::reduce([h0, h1, h2, h3, h4, h5, h6, h7, h8, h9]) + let x1_2 = 2 * x[1]; // This fits in a u32 iff 25 + b + 1 < 32 + let x3_2 = 2 * x[3]; // iff b < 6 + let x5_2 = 2 * x[5]; + let x7_2 = 2 * x[7]; + let x9_2 = 2 * x[9]; + + let z0 = m(x[0],y[0]) + m(x1_2,y9_19) + m(x[2],y8_19) + m(x3_2,y7_19) + m(x[4],y6_19) + m(x5_2,y5_19) + m(x[6],y4_19) + m(x7_2,y3_19) + m(x[8],y2_19) + m(x9_2,y1_19); + let z1 = m(x[0],y[1]) + m(x[1],y[0]) + m(x[2],y9_19) + m(x[3],y8_19) + m(x[4],y7_19) + m(x[5],y6_19) + m(x[6],y5_19) + m(x[7],y4_19) + m(x[8],y3_19) + m(x[9],y2_19); + let z2 = m(x[0],y[2]) + m(x1_2,y[1]) + m(x[2],y[0]) + m(x3_2,y9_19) + m(x[4],y8_19) + m(x5_2,y7_19) + m(x[6],y6_19) + m(x7_2,y5_19) + m(x[8],y4_19) + m(x9_2,y3_19); + let z3 = m(x[0],y[3]) + m(x[1],y[2]) + m(x[2],y[1]) + m(x[3],y[0]) + m(x[4],y9_19) + m(x[5],y8_19) + m(x[6],y7_19) + m(x[7],y6_19) + m(x[8],y5_19) + m(x[9],y4_19); + let z4 = m(x[0],y[4]) + m(x1_2,y[3]) + m(x[2],y[2]) + m(x3_2,y[1]) + m(x[4],y[0]) + m(x5_2,y9_19) + m(x[6],y8_19) + m(x7_2,y7_19) + m(x[8],y6_19) + m(x9_2,y5_19); + let z5 = m(x[0],y[5]) + m(x[1],y[4]) + m(x[2],y[3]) + m(x[3],y[2]) + m(x[4],y[1]) + m(x[5],y[0]) + m(x[6],y9_19) + m(x[7],y8_19) + m(x[8],y7_19) + m(x[9],y6_19); + let z6 = m(x[0],y[6]) + m(x1_2,y[5]) + m(x[2],y[4]) + m(x3_2,y[3]) + m(x[4],y[2]) + m(x5_2,y[1]) + m(x[6],y[0]) + m(x7_2,y9_19) + m(x[8],y8_19) + m(x9_2,y7_19); + let z7 = m(x[0],y[7]) + m(x[1],y[6]) + m(x[2],y[5]) + m(x[3],y[4]) + m(x[4],y[3]) + m(x[5],y[2]) + m(x[6],y[1]) + m(x[7],y[0]) + m(x[8],y9_19) + m(x[9],y8_19); + let z8 = m(x[0],y[8]) + m(x1_2,y[7]) + m(x[2],y[6]) + m(x3_2,y[5]) + m(x[4],y[4]) + m(x5_2,y[3]) + m(x[6],y[2]) + m(x7_2,y[1]) + m(x[8],y[0]) + m(x9_2,y9_19); + let z9 = m(x[0],y[9]) + m(x[1],y[8]) + m(x[2],y[7]) + m(x[3],y[6]) + m(x[4],y[5]) + m(x[5],y[4]) + m(x[6],y[3]) + m(x[7],y[2]) + m(x[8],y[1]) + m(x[9],y[0]); + + // How big is the contribution to z[i+j] from x[i], y[j]? + // + // Using the bounds above, we get: + // + // i even, j even: x[i]*y[j] < 2^(26+b)*2^(26+b) = 2*2^(51+2*b) + // i odd, j even: x[i]*y[j] < 2^(25+b)*2^(26+b) = 1*2^(51+2*b) + // i even, j odd: x[i]*y[j] < 2^(26+b)*2^(25+b) = 1*2^(51+2*b) + // i odd, j odd: 2*x[i]*y[j] < 2*2^(25+b)*2^(25+b) = 1*2^(51+2*b) + // + // We perform inline reduction mod p by replacing 2^255 by 19 + // (since 2^255 - 19 = 0 mod p). This adds a factor of 19, so + // we get the bounds (z0 is the biggest one, but calculated for + // posterity here in case finer estimation is needed later): + // + // z0 < ( 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 249*2^(51 + 2*b) + // z1 < ( 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 154*2^(51 + 2*b) + // z2 < ( 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 195*2^(51 + 2*b) + // z3 < ( 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 118*2^(51 + 2*b) + // z4 < ( 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 141*2^(51 + 2*b) + // z5 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 82*2^(51 + 2*b) + // z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 87*2^(51 + 2*b) + // z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 )*2^(51 + 2b) = 46*2^(51 + 2*b) + // z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 )*2^(51 + 2b) = 33*2^(51 + 2*b) + // z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 )*2^(51 + 2b) = 10*2^(51 + 2*b) + // + // So z[0] fits into a u64 if 51 + 2*b + lg(249) < 64 + // if b < 2.5. + FieldElement32::reduce([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9]) } } @@ -202,7 +227,7 @@ impl<'a> Neg for &'a FieldElement32 { impl ConditionallyAssignable for FieldElement32 { fn conditional_assign(&mut self, f: &FieldElement32, choice: u8) { - let mask = -(choice as i32); + let mask = (-(choice as i32)) as u32; for i in 0..10 { self.0[i] ^= mask & (self.0[i] ^ f.0[i]); } @@ -212,9 +237,20 @@ impl ConditionallyAssignable for FieldElement32 { impl FieldElement32 { /// Invert the sign of this field element pub fn negate(&mut self) { - for i in 0..10 { - self.0[i] = -self.0[i]; - } + // Compute -b as ((2^4 * p) - b) to avoid underflow. + let neg = FieldElement32::reduce([ + ((0x3ffffed << 4) - self.0[0]) as u64, + ((0x1ffffff << 4) - self.0[1]) as u64, + ((0x3ffffff << 4) - self.0[2]) as u64, + ((0x1ffffff << 4) - self.0[3]) as u64, + ((0x3ffffff << 4) - self.0[4]) as u64, + ((0x1ffffff << 4) - self.0[5]) as u64, + ((0x3ffffff << 4) - self.0[6]) as u64, + ((0x1ffffff << 4) - self.0[7]) as u64, + ((0x3ffffff << 4) - self.0[8]) as u64, + ((0x1ffffff << 4) - self.0[9]) as u64, + ]); + self.0 = neg.0; } /// Construct zero. @@ -229,98 +265,64 @@ impl FieldElement32 { /// Construct -1. pub fn minus_one() -> FieldElement32 { - FieldElement32([-1, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]) + FieldElement32([ + 0x3ffffec, 0x1ffffff, 0x3ffffff, 0x1ffffff, 0x3ffffff, + 0x1ffffff, 0x3ffffff, 0x1ffffff, 0x3ffffff, 0x1ffffff, + ]) } - fn reduce(mut h: [i64; 10]) -> FieldElement32 { //FeCombine - let mut c = [0i64; 10]; + /// Given unreduced coefficients `z[0], ..., z[9]` of any size, + /// carry and reduce them mod p to obtain a `FieldElement32` + /// whose coefficients have excess `b < 0.007`. + /// + /// In other words, each coefficient of the result is bounded by + /// either `2^(25 + 0.007)` or `2^(26 + 0.007)`, as appropriate. + fn reduce(mut z: [u64; 10]) -> FieldElement32 { - /* - |h[0]| <= (1.1*1.1*2^52*(1+19+19+19+19)+1.1*1.1*2^50*(38+38+38+38+38)) - i.e. |h[0]| <= 1.2*2^59; narrower ranges for h[2], h[4], h[6], h[8] - |h[1]| <= (1.1*1.1*2^51*(1+1+19+19+19+19+19+19+19+19)) - i.e. |h[1]| <= 1.5*2^58; narrower ranges for h[3], h[5], h[7], h[9] - */ + const LOW_25_BITS: u64 = (1 << 25) - 1; + const LOW_26_BITS: u64 = (1 << 26) - 1; - c[0] = (h[0] + (1 << 25)) >> 26; - h[1] += c[0]; - h[0] -= c[0] << 26; - c[4] = (h[4] + (1 << 25)) >> 26; - h[5] += c[4]; - h[4] -= c[4] << 26; - /* |h[0]| <= 2^25 */ - /* |h[4]| <= 2^25 */ - /* |h[1]| <= 1.51*2^58 */ - /* |h[5]| <= 1.51*2^58 */ + /// Carry the value from limb i = 0..8 to limb i+1 + #[inline(always)] + fn carry(z: &mut [u64; 10], i: usize) { + debug_assert!(i < 9); + if i % 2 == 0 { + // Even limbs have 26 bits + z[i+1] += z[i] >> 26; + z[i] &= LOW_26_BITS; + } else { + // Odd limbs have 25 bits + z[i+1] += z[i] >> 25; + z[i] &= LOW_25_BITS; + } + } - c[1] = (h[1] + (1 << 24)) >> 25; - h[2] += c[1]; - h[1] -= c[1] << 25; - c[5] = (h[5] + (1 << 24)) >> 25; - h[6] += c[5]; - h[5] -= c[5] << 25; - /* |h[1]| <= 2^24; from now on fits into int32 */ - /* |h[5]| <= 2^24; from now on fits into int32 */ - /* |h[2]| <= 1.21*2^59 */ - /* |h[6]| <= 1.21*2^59 */ + // Perform two halves of the carry chain in parallel. + carry(&mut z, 0); carry(&mut z, 4); + carry(&mut z, 1); carry(&mut z, 5); + carry(&mut z, 2); carry(&mut z, 6); + carry(&mut z, 3); carry(&mut z, 7); + // Since z[3] < 2^64, c < 2^(64-25) = 2^39, + // so z[4] < 2^26 + 2^39 < 2^39.0002 + carry(&mut z, 4); carry(&mut z, 8); + // Now z[4] < 2^26 + // and z[5] < 2^25 + 2^13.0002 < 2^25.0004 (good enough) - c[2] = (h[2] + (1 << 25)) >> 26; - h[3] += c[2]; - h[2] -= c[2] << 26; - c[6] = (h[6] + (1 << 25)) >> 26; - h[7] += c[6]; - h[6] -= c[6] << 26; - /* |h[2]| <= 2^25; from now on fits into int32 unchanged */ - /* |h[6]| <= 2^25; from now on fits into int32 unchanged */ - /* |h[3]| <= 1.51*2^58 */ - /* |h[7]| <= 1.51*2^58 */ + // Last carry has a multiplication by 19: + z[0] += 19*(z[9] >> 25); + z[9] &= LOW_25_BITS; - c[3] = (h[3] + (1 << 24)) >> 25; - h[4] += c[3]; - h[3] -= c[3] << 25; - c[7] = (h[7] + (1 << 24)) >> 25; - h[8] += c[7]; - h[7] -= c[7] << 25; - /* |h[3]| <= 2^24; from now on fits into int32 unchanged */ - /* |h[7]| <= 2^24; from now on fits into int32 unchanged */ - /* |h[4]| <= 1.52*2^33 */ - /* |h[8]| <= 1.52*2^33 */ + // Since z[9] < 2^64, c < 2^(64-25) = 2^39, + // so z[0] + 19*c < 2^26 + 2^43.248 < 2^43.249. + carry(&mut z, 0); + // Now z[1] < 2^25 - 2^(43.249 - 26) + // < 2^25.007 (good enough) + // and we're done. - c[4] = (h[4] + (1 << 25)) >> 26; - h[5] += c[4]; - h[4] -= c[4] << 26; - c[8] = (h[8] + (1 << 25)) >> 26; - h[9] += c[8]; - h[8] -= c[8] << 26; - /* |h[4]| <= 2^25; from now on fits into int32 unchanged */ - /* |h[8]| <= 2^25; from now on fits into int32 unchanged */ - /* |h[5]| <= 1.01*2^24 */ - /* |h[9]| <= 1.51*2^58 */ - - c[9] = (h[9] + (1 << 24)) >> 25; - h[0] += c[9] * 19; - h[9] -= c[9] << 25; - /* |h[9]| <= 2^24; from now on fits into int32 unchanged */ - /* |h[0]| <= 1.8*2^37 */ - - c[0] = (h[0] + (1 << 25)) >> 26; - h[1] += c[0]; - h[0] -= c[0] << 26; - /* |h[0]| <= 2^25; from now on fits into int32 unchanged */ - /* |h[1]| <= 1.01*2^24 */ - - let mut output = FieldElement32([0i32; 10]); - output.0[0] = h[0] as i32; - output.0[1] = h[1] as i32; - output.0[2] = h[2] as i32; - output.0[3] = h[3] as i32; - output.0[4] = h[4] as i32; - output.0[5] = h[5] as i32; - output.0[6] = h[6] as i32; - output.0[7] = h[7] as i32; - output.0[8] = h[8] as i32; - output.0[9] = h[9] as i32; - output + FieldElement32([ + z[0] as u32, z[1] as u32, z[2] as u32, z[3] as u32, z[4] as u32, + z[5] as u32, z[6] as u32, z[7] as u32, z[8] as u32, z[9] as u32, + ]) } /// Load a `FieldElement64` from the low 255 bits of a 256-bit @@ -334,11 +336,19 @@ impl FieldElement32 { /// encoding of every field element should decode, re-encode to /// the canonical encoding, and check that the input was /// canonical. - /// - /// XXX the above applies to the 64-bit implementation; check that - /// it applies here too. pub fn from_bytes(data: &[u8; 32]) -> FieldElement32 { //FeFromBytes - let mut h = [0i64;10]; + #[inline] + fn load3(b: &[u8]) -> u64 { + (b[0] as u64) | ((b[1] as u64) << 8) | ((b[2] as u64) << 16) + } + + #[inline] + fn load4(b: &[u8]) -> u64 { + (b[0] as u64) | ((b[1] as u64) << 8) | ((b[2] as u64) << 16) | ((b[3] as u64) << 24) + } + + let mut h = [0u64;10]; + const LOW_23_BITS: u64 = (1 << 23) - 1; h[0] = load4(&data[ 0..]); h[1] = load3(&data[ 4..]) << 6; h[2] = load3(&data[ 7..]) << 5; @@ -348,52 +358,33 @@ impl FieldElement32 { h[6] = load3(&data[20..]) << 7; h[7] = load3(&data[23..]) << 5; h[8] = load3(&data[26..]) << 4; - h[9] = (load3(&data[29..]) & 8388607) << 2; + h[9] = (load3(&data[29..]) & LOW_23_BITS) << 2; FieldElement32::reduce(h) } /// Serialize this `FieldElement64` to a 32-byte array. The /// encoding is canonical. - pub fn to_bytes(&self) -> [u8; 32] { //FeToBytes - // Comment preserved from ed25519.go (presumably originally from ref10): - // - // # Preconditions - // - // * `|h[i]|` bounded by 1.1*2^25, 1.1*2^24, 1.1*2^25, 1.1*2^24, etc. - // - // # Lemma - // - // Write p = 2^255 - 19 and q = floor(h/p). - // - // Basic claim: q = floor(2^(-255)(h + 19 * 2^-25 h9 + 2^-1)). - // - // # Proof - // - // Have |h|<=p so |q|<=1 so |19^2 * 2^-255 * q| < 1/4. - // - // Also have |h-2^230 * h9| < 2^230 so |19 * 2^-255 * (h-2^230 * h9)| < 1/4. - // - // Write y=2^(-1)-19^2 2^(-255)q-19 2^(-255)(h-2^230 h9), then 0 [u8; 32] { - let mut q:i32 = (19*h[9] + (1 << 24)) >> 25; - q = (h[0] + q) >> 26; + let inp = &self.0; + // Reduce the value represented by `in` to the range [0,2*p) + let mut h: [u32; 10] = FieldElement32::reduce([ + // XXX this cast is annoying + inp[0] as u64, inp[1] as u64, inp[2] as u64, inp[3] as u64, inp[4] as u64, + inp[5] as u64, inp[6] as u64, inp[7] as u64, inp[8] as u64, inp[9] as u64, + ]).0; + + // Let h be the value to encode. + // + // Write h = pq + r with 0 <= r < p. We want to compute r = h mod p. + // + // Since h < 2*p, q = 0 or 1, with q = 0 when h < p and q = 1 when h >= p. + // + // Notice that h >= p <==> h + 19 >= p + 19 <==> h + 19 >= 2^255. + // Therefore q can be computed as the carry bit of h + 19. + + let mut q: u32 = (h[0] + 19) >> 26; q = (h[1] + q) >> 25; q = (h[2] + q) >> 26; q = (h[3] + q) >> 25; @@ -404,45 +395,40 @@ impl FieldElement32 { q = (h[8] + q) >> 26; q = (h[9] + q) >> 25; - // Goal: Output h-(2^255-19)q, which is between 0 and 2^255-20. - h[0] += 19 * q; - // Goal: Output h-2^255 q, which is between 0 and 2^255-20. + debug_assert!( q == 0 || q == 1 ); - carry[0] = h[0] >> 26; - h[1] += carry[0]; - h[0] -= carry[0] << 26; - carry[1] = h[1] >> 25; - h[2] += carry[1]; - h[1] -= carry[1] << 25; - carry[2] = h[2] >> 26; - h[3] += carry[2]; - h[2] -= carry[2] << 26; - carry[3] = h[3] >> 25; - h[4] += carry[3]; - h[3] -= carry[3] << 25; - carry[4] = h[4] >> 26; - h[5] += carry[4]; - h[4] -= carry[4] << 26; - carry[5] = h[5] >> 25; - h[6] += carry[5]; - h[5] -= carry[5] << 25; - carry[6] = h[6] >> 26; - h[7] += carry[6]; - h[6] -= carry[6] << 26; - carry[7] = h[7] >> 25; - h[8] += carry[7]; - h[7] -= carry[7] << 25; - carry[8] = h[8] >> 26; - h[9] += carry[8]; - h[8] -= carry[8] << 26; - carry[9] = h[9] >> 25; - h[9] -= carry[9] << 25; - // h10 = carry9 + // Now we can compute r as r = h - pq = r - (2^255-19)q = r + 19q - 2^255q - // Goal: Output h[0]+...+2^255 h10-2^255 q, which is between 0 and 2^255-20. - // Have h[0]+...+2^230 h[9] between 0 and 2^255-1; - // evidently 2^255 h10-2^255 q = 0. - // Goal: Output h[0]+...+2^230 h[9]. + const LOW_25_BITS: u32 = (1 << 25) - 1; + const LOW_26_BITS: u32 = (1 << 26) - 1; + + h[0] += 19*q; + + // Now carry the result to compute r + 19q... + h[1] += h[0] >> 26; + h[0] = h[0] & LOW_26_BITS; + h[2] += h[1] >> 25; + h[1] = h[1] & LOW_25_BITS; + h[3] += h[2] >> 26; + h[2] = h[2] & LOW_26_BITS; + h[4] += h[3] >> 25; + h[3] = h[3] & LOW_25_BITS; + h[5] += h[4] >> 26; + h[4] = h[4] & LOW_26_BITS; + h[6] += h[5] >> 25; + h[5] = h[5] & LOW_25_BITS; + h[7] += h[6] >> 26; + h[6] = h[6] & LOW_26_BITS; + h[8] += h[7] >> 25; + h[7] = h[7] & LOW_25_BITS; + h[9] += h[8] >> 26; + h[8] = h[8] & LOW_26_BITS; + + // ... but instead of carrying the value + // (h[9] >> 25) = q*2^255 into another limb, + // discard it, subtracting the value from h. + debug_assert!( (h[9] >> 25) == 0 || (h[9] >> 25) == 1); + h[9] = h[9] & LOW_25_BITS; let mut s = [0u8; 32]; s[0] = (h[0] >> 0) as u8; @@ -484,77 +470,53 @@ impl FieldElement32 { s } - fn square_inner(&self) -> [i64; 10] { - let f0 = self.0[0] as i64; - let f1 = self.0[1] as i64; - let f2 = self.0[2] as i64; - let f3 = self.0[3] as i64; - let f4 = self.0[4] as i64; - let f5 = self.0[5] as i64; - let f6 = self.0[6] as i64; - let f7 = self.0[7] as i64; - let f8 = self.0[8] as i64; - let f9 = self.0[9] as i64; - let f0_2 = (2 * self.0[0]) as i64; - let f1_2 = (2 * self.0[1]) as i64; - let f2_2 = (2 * self.0[2]) as i64; - let f3_2 = (2 * self.0[3]) as i64; - let f4_2 = (2 * self.0[4]) as i64; - let f5_2 = (2 * self.0[5]) as i64; - let f6_2 = (2 * self.0[6]) as i64; - let f7_2 = (2 * self.0[7]) as i64; - let f5_38 = 38 * f5; // 1.31*2^30 - let f6_19 = 19 * f6; // 1.31*2^30 - let f7_38 = 38 * f7; // 1.31*2^30 - let f8_19 = 19 * f8; // 1.31*2^30 - let f9_38 = 38 * f9; // 1.31*2^30 + fn square_inner(&self) -> [u64; 10] { + // Optimized version of multiplication for the case of squaring. + // Pre- and post- conditions identical to multiplication function. + let x = &self.0; + let x0_2 = 2 * x[0]; + let x1_2 = 2 * x[1]; + let x2_2 = 2 * x[2]; + let x3_2 = 2 * x[3]; + let x4_2 = 2 * x[4]; + let x5_2 = 2 * x[5]; + let x6_2 = 2 * x[6]; + let x7_2 = 2 * x[7]; + let x5_19 = 19 * x[5]; + let x6_19 = 19 * x[6]; + let x7_19 = 19 * x[7]; + let x8_19 = 19 * x[8]; + let x9_19 = 19 * x[9]; - let mut h = [0i64;10]; - h[0] = f0*f0 + f1_2*f9_38 + f2_2*f8_19 + f3_2*f7_38 + f4_2*f6_19 + f5*f5_38; - h[1] = f0_2*f1 + f2*f9_38 + f3_2*f8_19 + f4*f7_38 + f5_2*f6_19; - h[2] = f0_2*f2 + f1_2*f1 + f3_2*f9_38 + f4_2*f8_19 + f5_2*f7_38 + f6*f6_19; - h[3] = f0_2*f3 + f1_2*f2 + f4*f9_38 + f5_2*f8_19 + f6*f7_38; - h[4] = f0_2*f4 + f1_2*f3_2 + f2*f2 + f5_2*f9_38 + f6_2*f8_19 + f7*f7_38; - h[5] = f0_2*f5 + f1_2*f4 + f2_2*f3 + f6*f9_38 + f7_2*f8_19; - h[6] = f0_2*f6 + f1_2*f5_2 + f2_2*f4 + f3_2*f3 + f7_2*f9_38 + f8*f8_19; - h[7] = f0_2*f7 + f1_2*f6 + f2_2*f5 + f3_2*f4 + f8*f9_38; - h[8] = f0_2*f8 + f1_2*f7_2 + f2_2*f6 + f3_2*f5_2 + f4*f4 + f9*f9_38; - h[9] = f0_2*f9 + f1_2*f8 + f2_2*f7 + f3_2*f6 + f4_2*f5; + /// Helper function to multiply two 32-bit integers with 64 bits + /// of output. + #[inline(always)] + fn m(x: u32, y: u32) -> u64 { (x as u64) * (y as u64) } - h + // This block is rearranged so that instead of doing a 32-bit multiplication by 38, we do a + // 64-bit multiplication by 2 on the results. This is because lg(38) is too big: we would + // have less than 1 bit of headroom left, which is too little. + let mut z = [0u64;10]; + z[0] = m(x[0],x[0]) + m(x2_2,x8_19) + m(x4_2,x6_19) + (m(x1_2,x9_19) + m(x3_2,x7_19) + m(x[5],x5_19))*2; + z[1] = m(x0_2,x[1]) + m(x3_2,x8_19) + m(x5_2,x6_19) + (m(x[2],x9_19) + m(x[4],x7_19))*2; + z[2] = m(x0_2,x[2]) + m(x1_2,x[1]) + m(x4_2,x8_19) + m(x[6],x6_19) + (m(x3_2,x9_19) + m(x5_2,x7_19))*2; + z[3] = m(x0_2,x[3]) + m(x1_2,x[2]) + m(x5_2,x8_19) + (m(x[4],x9_19) + m(x[6],x7_19))*2; + z[4] = m(x0_2,x[4]) + m(x1_2,x3_2) + m(x[2],x[2]) + m(x6_2,x8_19) + (m(x5_2,x9_19) + m(x[7],x7_19))*2; + z[5] = m(x0_2,x[5]) + m(x1_2,x[4]) + m(x2_2,x[3]) + m(x7_2,x8_19) + m(x[6],x9_19)*2; + z[6] = m(x0_2,x[6]) + m(x1_2,x5_2) + m(x2_2,x[4]) + m(x3_2,x[3]) + m(x[8],x8_19) + m(x7_2,x9_19)*2; + z[7] = m(x0_2,x[7]) + m(x1_2,x[6]) + m(x2_2,x[5]) + m(x3_2,x[4]) + m(x[8],x9_19)*2; + z[8] = m(x0_2,x[8]) + m(x1_2,x7_2) + m(x2_2,x[6]) + m(x3_2,x5_2) + m(x[4],x[4]) + m(x[9],x9_19)*2; + z[9] = m(x0_2,x[9]) + m(x1_2,x[8]) + m(x2_2,x[7]) + m(x3_2,x[6]) + m(x4_2,x[5]) ; + + z } - /// Calculates h = f*f. Can overlap h with f. - /// - /// XXX limbs: better to talk about headroom? - /// - /// # Preconditions - /// - /// * |f[i]| bounded by 1.1*2^26, 1.1*2^25, 1.1*2^26, 1.1*2^25, etc. - /// - /// # Postconditions - /// - /// * |h[i]| bounded by 1.1*2^25, 1.1*2^24, 1.1*2^25, 1.1*2^24, etc. + /// Compute `self^2`. pub fn square(&self) -> FieldElement32 { FieldElement32::reduce(self.square_inner()) } - /// Square this field element and multiply the result by 2. - /// - /// XXX explain why square2 exists vs square (overflow) - /// - /// # Preconditions - /// - /// * |f[i]| bounded by 1.65*2^26, 1.65*2^25, 1.65*2^26, 1.65*2^25, etc. - /// - /// # Postconditions - /// - /// * |h[i]| bounded by 1.01*2^25, 1.01*2^24, 1.01*2^25, 1.01*2^24, etc. - /// - /// # Notes - /// - /// See fe_mul.c in ref10 implementation for discussion of implementation - /// strategy. + /// Compute `2*self^2`. pub fn square2(&self) -> FieldElement32 { let mut coeffs = self.square_inner(); for i in 0..self.0.len() { diff --git a/src/lib.rs b/src/lib.rs index ad62fdf..2c86b46 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -82,7 +82,7 @@ pub mod ristretto; // Other miscelaneous utilities. -pub mod utils; +mod utils; // Low-level curve and point constants, as well as pre-computed curve group elements. @@ -91,4 +91,3 @@ pub mod constants; mod constants_32bit; #[cfg(feature="radix_51")] mod constants_64bit; -