// -*- mode: rust; coding: utf-8; -*- // // This file is part of curve25519-dalek. // Copyright (c) 2016-2018 Isis Lovecruft, Henry de Valence // See LICENSE for licensing information. // // Authors: // - Isis Agora Lovecruft // - Henry de Valence //! Field arithmetic modulo \\(p = 2\^{255} - 19\\), using \\(64\\)-bit //! limbs with \\(128\\)-bit products. use core::fmt::Debug; use core::ops::{Add, AddAssign}; use core::ops::{Sub, SubAssign}; use core::ops::{Mul, MulAssign}; use core::ops::Neg; use subtle::ConditionallyAssignable; use subtle::Choice; /// A `FieldElement64` represents an element of the field /// \\( \mathbb Z / (2\^{255} - 19)\\). /// /// In the 64-bit implementation, a `FieldElement` is represented in /// radix \\(2\^{51}\\) as five `u64`s; the coefficients are allowed to /// grow up to \\(2\^{54}\\) between reductions modulo \\(p\\). /// /// # Note /// /// The `curve25519_dalek::field` module provides a type alias /// `curve25519_dalek::field::FieldElement` to either `FieldElement64` /// or `FieldElement32`. /// /// The backend-specific type `FieldElement64` should not be used /// outside of the `curve25519_dalek::field` module. #[derive(Copy, Clone)] pub struct FieldElement64(pub (crate) [u64; 5]); impl Debug for FieldElement64 { fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result { write!(f, "FieldElement64({:?})", &self.0[..]) } } impl<'b> AddAssign<&'b FieldElement64> for FieldElement64 { fn add_assign(&mut self, _rhs: &'b FieldElement64) { for i in 0..5 { self.0[i] += _rhs.0[i]; } } } impl<'a, 'b> Add<&'b FieldElement64> for &'a FieldElement64 { type Output = FieldElement64; fn add(self, _rhs: &'b FieldElement64) -> FieldElement64 { let mut output = *self; output += _rhs; output } } impl<'b> SubAssign<&'b FieldElement64> for FieldElement64 { fn sub_assign(&mut self, _rhs: &'b FieldElement64) { let result = (self as &FieldElement64) - _rhs; self.0 = result.0; } } impl<'a, 'b> Sub<&'b FieldElement64> for &'a FieldElement64 { type Output = FieldElement64; fn sub(self, _rhs: &'b FieldElement64) -> FieldElement64 { // To avoid underflow, first add a multiple of p. // Choose 16*p = p << 4 to be larger than 54-bit _rhs. // // If we could statically track the bitlengths of the limbs // of every FieldElement64, we could choose a multiple of p // just bigger than _rhs and avoid having to do a reduction. // // Since we don't yet have type-level integers to do this, we // have to add an explicit reduction call here. FieldElement64::reduce([ (self.0[0] + 36028797018963664u64) - _rhs.0[0], (self.0[1] + 36028797018963952u64) - _rhs.0[1], (self.0[2] + 36028797018963952u64) - _rhs.0[2], (self.0[3] + 36028797018963952u64) - _rhs.0[3], (self.0[4] + 36028797018963952u64) - _rhs.0[4], ]) } } impl<'b> MulAssign<&'b FieldElement64> for FieldElement64 { fn mul_assign(&mut self, _rhs: &'b FieldElement64) { let result = (self as &FieldElement64) * _rhs; self.0 = result.0; } } impl<'a, 'b> Mul<&'b FieldElement64> for &'a FieldElement64 { type Output = FieldElement64; fn mul(self, _rhs: &'b FieldElement64) -> FieldElement64 { /// Helper function to multiply two 64-bit integers with 128 /// bits of output. #[inline(always)] fn m(x: u64, y: u64) -> u128 { (x as u128) * (y as u128) } // Alias self, _rhs for more readable formulas let a: &[u64; 5] = &self.0; let b: &[u64; 5] = &_rhs.0; // Precondition: assume input limbs a[i], b[i] are bounded as // // a[i], b[i] < 2^(51 + b) // // where b is a real parameter measuring the "bit excess" of the limbs. // 64-bit precomputations to avoid 128-bit multiplications. // // This fits into a u64 whenever 51 + b + lg(19) < 64. // // Since 51 + b + lg(19) < 51 + 4.25 + b // = 55.25 + b, // this fits if b < 8.75. let b1_19 = b[1] * 19; let b2_19 = b[2] * 19; let b3_19 = b[3] * 19; let b4_19 = b[4] * 19; // Multiply to get 128-bit coefficients of output let c0: u128 = m(a[0],b[0]) + m(a[4],b1_19) + m(a[3],b2_19) + m(a[2],b3_19) + m(a[1],b4_19); let mut c1: u128 = m(a[1],b[0]) + m(a[0],b[1]) + m(a[4],b2_19) + m(a[3],b3_19) + m(a[2],b4_19); let mut c2: u128 = m(a[2],b[0]) + m(a[1],b[1]) + m(a[0],b[2]) + m(a[4],b3_19) + m(a[3],b4_19); let mut c3: u128 = m(a[3],b[0]) + m(a[2],b[1]) + m(a[1],b[2]) + m(a[0],b[3]) + m(a[4],b4_19); let mut c4: u128 = m(a[4],b[0]) + m(a[3],b[1]) + m(a[2],b[2]) + m(a[1],b[3]) + m(a[0],b[4]); // How big are the c[i]? We have // // c[i] < 2^(102 + 2*b) * (1+i + (4-i)*19) // < 2^(102 + lg(1 + 4*19) + 2*b) // < 2^(108.27 + 2*b) // // The carry (c[i] >> 51) fits into a u64 when // 108.27 + 2*b - 51 < 64 // 2*b < 6.73 // b < 3.365. // // So we require b < 3 to ensure this fits. debug_assert!(a[0] < (1 << 54)); debug_assert!(b[0] < (1 << 54)); debug_assert!(a[1] < (1 << 54)); debug_assert!(b[1] < (1 << 54)); debug_assert!(a[2] < (1 << 54)); debug_assert!(b[2] < (1 << 54)); debug_assert!(a[3] < (1 << 54)); debug_assert!(b[3] < (1 << 54)); debug_assert!(a[4] < (1 << 54)); debug_assert!(b[4] < (1 << 54)); // Casting to u64 and back tells the compiler that the carry is // bounded by 2^64, so that the addition is a u128 + u64 rather // than u128 + u128. const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1; let mut out = [0u64; 5]; c1 += ((c0 >> 51) as u64) as u128; out[0] = (c0 as u64) & LOW_51_BIT_MASK; c2 += ((c1 >> 51) as u64) as u128; out[1] = (c1 as u64) & LOW_51_BIT_MASK; c3 += ((c2 >> 51) as u64) as u128; out[2] = (c2 as u64) & LOW_51_BIT_MASK; c4 += ((c3 >> 51) as u64) as u128; out[3] = (c3 as u64) & LOW_51_BIT_MASK; let carry: u64 = (c4 >> 51) as u64; out[4] = (c4 as u64) & LOW_51_BIT_MASK; // To see that this does not overflow, we need out[0] + carry * 19 < 2^64. // // c4 < a0*b4 + a1*b3 + a2*b2 + a3*b1 + a4*b0 + (carry from c3) // < 5*(2^(51 + b) * 2^(51 + b)) + (carry from c3) // < 2^(102 + 2*b + lg(5)) + 2^64. // // When b < 3 we get // // c4 < 2^110.33 so that carry < 2^59.33 // // so that // // out[0] + carry * 19 < 2^51 + 19 * 2^59.33 < 2^63.58 // // and there is no overflow. out[0] = out[0] + carry * 19; // Now out[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + epsilon). out[1] += out[0] >> 51; out[0] &= LOW_51_BIT_MASK; // Now out[i] < 2^(51 + epsilon) for all i. FieldElement64(out) } } impl<'a> Neg for &'a FieldElement64 { type Output = FieldElement64; fn neg(self) -> FieldElement64 { let mut output = *self; output.negate(); output } } impl ConditionallyAssignable for FieldElement64 { fn conditional_assign(&mut self, other: &FieldElement64, choice: Choice) { for i in 0..5 { self.0[i].conditional_assign(&other.0[i], choice); } } } impl FieldElement64 { /// Invert the sign of this field element pub fn negate(&mut self) { // See commentary in the Sub impl let neg = FieldElement64::reduce([ 36028797018963664u64 - self.0[0], 36028797018963952u64 - self.0[1], 36028797018963952u64 - self.0[2], 36028797018963952u64 - self.0[3], 36028797018963952u64 - self.0[4], ]); self.0 = neg.0; } /// Construct zero. pub fn zero() -> FieldElement64 { FieldElement64([ 0, 0, 0, 0, 0 ]) } /// Construct one. pub fn one() -> FieldElement64 { FieldElement64([ 1, 0, 0, 0, 0 ]) } /// Construct -1. pub fn minus_one() -> FieldElement64 { FieldElement64([2251799813685228, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247]) } /// Given 64-bit input limbs, reduce to enforce the bound 2^(51 + epsilon). #[inline(always)] fn reduce(mut limbs: [u64; 5]) -> FieldElement64 { const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1; // Since the input limbs are bounded by 2^64, the biggest // carry-out is bounded by 2^13. // // The biggest carry-in is c4 * 19, resulting in // // 2^51 + 19*2^13 < 2^51.0000000001 // // Because we don't need to canonicalize, only to reduce the // limb sizes, it's OK to do a "weak reduction", where we // compute the carry-outs in parallel. let c0 = limbs[0] >> 51; let c1 = limbs[1] >> 51; let c2 = limbs[2] >> 51; let c3 = limbs[3] >> 51; let c4 = limbs[4] >> 51; limbs[0] &= LOW_51_BIT_MASK; limbs[1] &= LOW_51_BIT_MASK; limbs[2] &= LOW_51_BIT_MASK; limbs[3] &= LOW_51_BIT_MASK; limbs[4] &= LOW_51_BIT_MASK; limbs[0] += c4 * 19; limbs[1] += c0; limbs[2] += c1; limbs[3] += c2; limbs[4] += c3; FieldElement64(limbs) } /// Load a `FieldElement64` from the low 255 bits of a 256-bit /// input. /// /// # Warning /// /// This function does not check that the input used the canonical /// representative. It masks the high bit, but it will happily /// decode 2^255 - 18 to 1. Applications that require a canonical /// encoding of every field element should decode, re-encode to /// the canonical encoding, and check that the input was /// canonical. /// pub fn from_bytes(bytes: &[u8; 32]) -> FieldElement64 { let load8 = |input: &[u8]| -> u64 { (input[0] as u64) | ((input[1] as u64) << 8) | ((input[2] as u64) << 16) | ((input[3] as u64) << 24) | ((input[4] as u64) << 32) | ((input[5] as u64) << 40) | ((input[6] as u64) << 48) | ((input[7] as u64) << 56) }; let low_51_bit_mask = (1u64 << 51) - 1; FieldElement64( // load bits [ 0, 64), no shift [ load8(&bytes[ 0..]) & low_51_bit_mask // load bits [ 48,112), shift to [ 51,112) , (load8(&bytes[ 6..]) >> 3) & low_51_bit_mask // load bits [ 96,160), shift to [102,160) , (load8(&bytes[12..]) >> 6) & low_51_bit_mask // load bits [152,216), shift to [153,216) , (load8(&bytes[19..]) >> 1) & low_51_bit_mask // load bits [192,256), shift to [204,112) , (load8(&bytes[24..]) >> 12) & low_51_bit_mask ]) } /// Serialize this `FieldElement64` to a 32-byte array. The /// encoding is canonical. pub fn to_bytes(&self) -> [u8; 32] { // Let h = limbs[0] + limbs[1]*2^51 + ... + limbs[4]*2^204. // // Write h = pq + r with 0 <= r < p. // // We want to compute r = h mod p. // // If h < 2*p = 2^256 - 38, // then 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. // First, reduce the limbs to ensure h < 2*p. let mut limbs = FieldElement64::reduce(self.0).0; let mut q = (limbs[0] + 19) >> 51; q = (limbs[1] + q) >> 51; q = (limbs[2] + q) >> 51; q = (limbs[3] + q) >> 51; q = (limbs[4] + q) >> 51; // Now we can compute r as r = h - pq = r - (2^255-19)q = r + 19q - 2^255q limbs[0] += 19*q; // Now carry the result to compute r + 19q ... let low_51_bit_mask = (1u64 << 51) - 1; limbs[1] += limbs[0] >> 51; limbs[0] = limbs[0] & low_51_bit_mask; limbs[2] += limbs[1] >> 51; limbs[1] = limbs[1] & low_51_bit_mask; limbs[3] += limbs[2] >> 51; limbs[2] = limbs[2] & low_51_bit_mask; limbs[4] += limbs[3] >> 51; limbs[3] = limbs[3] & low_51_bit_mask; // ... but instead of carrying (limbs[4] >> 51) = 2^255q // into another limb, discard it, subtracting the value limbs[4] = limbs[4] & low_51_bit_mask; // Now arrange the bits of the limbs. let mut s = [0u8;32]; s[ 0] = limbs[0] as u8; s[ 1] = (limbs[0] >> 8) as u8; s[ 2] = (limbs[0] >> 16) as u8; s[ 3] = (limbs[0] >> 24) as u8; s[ 4] = (limbs[0] >> 32) as u8; s[ 5] = (limbs[0] >> 40) as u8; s[ 6] = ((limbs[0] >> 48) | (limbs[1] << 3)) as u8; s[ 7] = (limbs[1] >> 5) as u8; s[ 8] = (limbs[1] >> 13) as u8; s[ 9] = (limbs[1] >> 21) as u8; s[10] = (limbs[1] >> 29) as u8; s[11] = (limbs[1] >> 37) as u8; s[12] = ((limbs[1] >> 45) | (limbs[2] << 6)) as u8; s[13] = (limbs[2] >> 2) as u8; s[14] = (limbs[2] >> 10) as u8; s[15] = (limbs[2] >> 18) as u8; s[16] = (limbs[2] >> 26) as u8; s[17] = (limbs[2] >> 34) as u8; s[18] = (limbs[2] >> 42) as u8; s[19] = ((limbs[2] >> 50) | (limbs[3] << 1)) as u8; s[20] = (limbs[3] >> 7) as u8; s[21] = (limbs[3] >> 15) as u8; s[22] = (limbs[3] >> 23) as u8; s[23] = (limbs[3] >> 31) as u8; s[24] = (limbs[3] >> 39) as u8; s[25] = ((limbs[3] >> 47) | (limbs[4] << 4)) as u8; s[26] = (limbs[4] >> 4) as u8; s[27] = (limbs[4] >> 12) as u8; s[28] = (limbs[4] >> 20) as u8; s[29] = (limbs[4] >> 28) as u8; s[30] = (limbs[4] >> 36) as u8; s[31] = (limbs[4] >> 44) as u8; // High bit should be zero. debug_assert!((s[31] & 0b1000_0000u8) == 0u8); s } /// Given `k > 0`, return `self^(2^k)`. pub fn pow2k(&self, mut k: u32) -> FieldElement64 { debug_assert!( k > 0 ); /// Multiply two 64-bit integers with 128 bits of output. #[inline(always)] fn m(x: u64, y: u64) -> u128 { (x as u128) * (y as u128) } let mut a: [u64; 5] = self.0; loop { // Precondition: assume input limbs a[i] are bounded as // // a[i] < 2^(51 + b) // // where b is a real parameter measuring the "bit excess" of the limbs. // Precomputation: 64-bit multiply by 19. // // This fits into a u64 whenever 51 + b + lg(19) < 64. // // Since 51 + b + lg(19) < 51 + 4.25 + b // = 55.25 + b, // this fits if b < 8.75. let a3_19 = 19 * a[3]; let a4_19 = 19 * a[4]; // Multiply to get 128-bit coefficients of output. // // The 128-bit multiplications by 2 turn into 1 slr + 1 slrd each, // which doesn't seem any better or worse than doing them as precomputations // on the 64-bit inputs. let c0: u128 = m(a[0], a[0]) + 2*( m(a[1], a4_19) + m(a[2], a3_19) ); let mut c1: u128 = m(a[3], a3_19) + 2*( m(a[0], a[1]) + m(a[2], a4_19) ); let mut c2: u128 = m(a[1], a[1]) + 2*( m(a[0], a[2]) + m(a[4], a3_19) ); let mut c3: u128 = m(a[4], a4_19) + 2*( m(a[0], a[3]) + m(a[1], a[2]) ); let mut c4: u128 = m(a[2], a[2]) + 2*( m(a[0], a[4]) + m(a[1], a[3]) ); // Same bound as in multiply: // c[i] < 2^(102 + 2*b) * (1+i + (4-i)*19) // < 2^(102 + lg(1 + 4*19) + 2*b) // < 2^(108.27 + 2*b) // // The carry (c[i] >> 51) fits into a u64 when // 108.27 + 2*b - 51 < 64 // 2*b < 6.73 // b < 3.365. // // So we require b < 3 to ensure this fits. debug_assert!(a[0] < (1 << 54)); debug_assert!(a[1] < (1 << 54)); debug_assert!(a[2] < (1 << 54)); debug_assert!(a[3] < (1 << 54)); debug_assert!(a[4] < (1 << 54)); const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1; // Casting to u64 and back tells the compiler that the carry is bounded by 2^64, so // that the addition is a u128 + u64 rather than u128 + u128. c1 += ((c0 >> 51) as u64) as u128; a[0] = (c0 as u64) & LOW_51_BIT_MASK; c2 += ((c1 >> 51) as u64) as u128; a[1] = (c1 as u64) & LOW_51_BIT_MASK; c3 += ((c2 >> 51) as u64) as u128; a[2] = (c2 as u64) & LOW_51_BIT_MASK; c4 += ((c3 >> 51) as u64) as u128; a[3] = (c3 as u64) & LOW_51_BIT_MASK; let carry: u64 = (c4 >> 51) as u64; a[4] = (c4 as u64) & LOW_51_BIT_MASK; // To see that this does not overflow, we need a[0] + carry * 19 < 2^64. // // c4 < a2^2 + 2*a0*a4 + 2*a1*a3 + (carry from c3) // < 2^(102 + 2*b + lg(5)) + 2^64. // // When b < 3 we get // // c4 < 2^110.33 so that carry < 2^59.33 // // so that // // a[0] + carry * 19 < 2^51 + 19 * 2^59.33 < 2^63.58 // // and there is no overflow. a[0] = a[0] + carry * 19; // Now a[1] < 2^51 + 2^(64 -51) = 2^51 + 2^13 < 2^(51 + epsilon). a[1] += a[0] >> 51; a[0] &= LOW_51_BIT_MASK; // Now all a[i] < 2^(51 + epsilon) and a = self^(2^k). k = k - 1; if k == 0 { break; } } FieldElement64(a) } /// Returns the square of this field element. pub fn square(&self) -> FieldElement64 { self.pow2k(1) } /// Returns 2 times the square of this field element. pub fn square2(&self) -> FieldElement64 { let mut square = self.pow2k(1); for i in 0..5 { square.0[i] *= 2; } square } }