saymrwulf/curve25519-dalek-source
1
0
Fork 0
mirror of https://github.com/saymrwulf/curve25519-dalek-source.git synced 2026-09-10 21:21:13 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions
curve25519-dalek-source/src/backend/serial/u64/field.rs
Henry de Valence f1d2b5182b Restructure source tree into serial and vector backends. 
This begins to attempt to restructure the source tree so that the common parts
are common and the different parts are different.

The backend is now split into two parts:
- serial (containing the implementation using serial formulas and mixed-model arithmetic).
- vector (containing the implementation using parallel formulas and single-model arithmetic).

The serial scalar_mul tree is now under backend::serial::scalar_mul.
The avx2 scalar_mul tree is now under backend::avx2::scalar_mul.
2019-01-18 01:49:40 -08:00

555 lines
19 KiB
Rust
Raw Blame History

// -*- 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 <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca>
//! Field arithmetic modulo \\(p = 2\^{255} - 19\\), using \\(64\\)-bit
//! limbs with \\(128\\)-bit products.
use core::fmt::Debug;
use core::ops::Neg;
use core::ops::{Add, AddAssign};
use core::ops::{Mul, MulAssign};
use core::ops::{Sub, SubAssign};
use subtle::Choice;
use subtle::ConditionallySelectable;
/// A `FieldElement51` 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 `FieldElement51`
/// or `FieldElement2625`.
///
/// The backend-specific type `FieldElement51` should not be used
/// outside of the `curve25519_dalek::field` module.
#[derive(Copy, Clone)]
pub struct FieldElement51(pub (crate) [u64; 5]);
impl Debug for FieldElement51 {
fn fmt(&self, f: &mut ::core::fmt::Formatter) -> ::core::fmt::Result {
write!(f, "FieldElement51({:?})", &self.0[..])
}
}
impl<'b> AddAssign<&'b FieldElement51> for FieldElement51 {
fn add_assign(&mut self, _rhs: &'b FieldElement51) {
for i in 0..5 {
self.0[i] += _rhs.0[i];
}
}
}
impl<'a, 'b> Add<&'b FieldElement51> for &'a FieldElement51 {
type Output = FieldElement51;
fn add(self, _rhs: &'b FieldElement51) -> FieldElement51 {
let mut output = *self;
output += _rhs;
output
}
}
impl<'b> SubAssign<&'b FieldElement51> for FieldElement51 {
fn sub_assign(&mut self, _rhs: &'b FieldElement51) {
let result = (self as &FieldElement51) - _rhs;
self.0 = result.0;
}
}
impl<'a, 'b> Sub<&'b FieldElement51> for &'a FieldElement51 {
type Output = FieldElement51;
fn sub(self, _rhs: &'b FieldElement51) -> FieldElement51 {
// 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 FieldElement51, 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.
FieldElement51::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 FieldElement51> for FieldElement51 {
fn mul_assign(&mut self, _rhs: &'b FieldElement51) {
let result = (self as &FieldElement51) * _rhs;
self.0 = result.0;
}
}
impl<'a, 'b> Mul<&'b FieldElement51> for &'a FieldElement51 {
type Output = FieldElement51;
fn mul(self, _rhs: &'b FieldElement51) -> FieldElement51 {
/// 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.
FieldElement51(out)
}
}
impl<'a> Neg for &'a FieldElement51 {
type Output = FieldElement51;
fn neg(self) -> FieldElement51 {
let mut output = *self;
output.negate();
output
}
}
impl ConditionallySelectable for FieldElement51 {
fn conditional_select(
a: &FieldElement51,
b: &FieldElement51,
choice: Choice,
) -> FieldElement51 {
FieldElement51([
u64::conditional_select(&a.0[0], &b.0[0], choice),
u64::conditional_select(&a.0[1], &b.0[1], choice),
u64::conditional_select(&a.0[2], &b.0[2], choice),
u64::conditional_select(&a.0[3], &b.0[3], choice),
u64::conditional_select(&a.0[4], &b.0[4], choice),
])
}
fn conditional_swap(a: &mut FieldElement51, b: &mut FieldElement51, choice: Choice) {
u64::conditional_swap(&mut a.0[0], &mut b.0[0], choice);
u64::conditional_swap(&mut a.0[1], &mut b.0[1], choice);
u64::conditional_swap(&mut a.0[2], &mut b.0[2], choice);
u64::conditional_swap(&mut a.0[3], &mut b.0[3], choice);
u64::conditional_swap(&mut a.0[4], &mut b.0[4], choice);
}
fn conditional_assign(&mut self, other: &FieldElement51, choice: Choice) {
self.0[0].conditional_assign(&other.0[0], choice);
self.0[1].conditional_assign(&other.0[1], choice);
self.0[2].conditional_assign(&other.0[2], choice);
self.0[3].conditional_assign(&other.0[3], choice);
self.0[4].conditional_assign(&other.0[4], choice);
}
}
impl FieldElement51 {
/// Invert the sign of this field element
pub fn negate(&mut self) {
// See commentary in the Sub impl
let neg = FieldElement51::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() -> FieldElement51 {
FieldElement51([ 0, 0, 0, 0, 0 ])
}
/// Construct one.
pub fn one() -> FieldElement51 {
FieldElement51([ 1, 0, 0, 0, 0 ])
}
/// Construct -1.
pub fn minus_one() -> FieldElement51 {
FieldElement51([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]) -> FieldElement51 {
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;
FieldElement51(limbs)
}
/// Load a `FieldElement51` 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]) -> FieldElement51 {
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;
FieldElement51(
// 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 `FieldElement51` 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 = FieldElement51::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) -> FieldElement51 {
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;
}
}
FieldElement51(a)
}
/// Returns the square of this field element.
pub fn square(&self) -> FieldElement51 {
self.pow2k(1)
}
/// Returns 2 times the square of this field element.
pub fn square2(&self) -> FieldElement51 {
let mut square = self.pow2k(1);
for i in 0..5 {
square.0[i] *= 2;
}
square
}
}
Reference in a new issue View git blame Copy permalink