diff --git a/CHANGELOG.md b/CHANGELOG.md index c2fa0ad..fc7a3ab 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,15 @@ Entries are listed in reverse chronological order. +## 1.2.0 + +* New multiscalar multiplication algorithm with better performance for + large problem sizes. The backend algorithm is selected + transparently using the size hints of the input iterators, so no + changes are required for client crates to start using it. +* Equality of Edwards points is now checked in projective coordinates. +* Serde can now be used with `no_std`. + ## 1.1.4 * Fix typos in documentation comments. diff --git a/Cargo.toml b/Cargo.toml index 6f3ed15..b2026d3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "curve25519-dalek" -version = "1.1.4" +version = "1.2.0" authors = ["Isis Lovecruft ", "Henry de Valence "] readme = "README.md" @@ -49,7 +49,7 @@ byteorder = { version = "^1.2.3", default-features = false, features = ["i128"] digest = { version = "0.8", default-features = false } clear_on_drop = "=0.2.3" subtle = { version = "2", default-features = false } -serde = { version = "1.0", optional = true } +serde = { version = "1.0", default-features = false, optional = true } packed_simd = { version = "0.3.0", features = ["into_bits"], optional = true } [build-dependencies] @@ -58,7 +58,7 @@ byteorder = { version = "^1.2.3", default-features = false, features = ["i128"] digest = { version = "0.8", default-features = false } clear_on_drop = "=0.2.3" subtle = { version = "2", default-features = false } -serde = { version = "1.0", optional = true } +serde = { version = "1.0", default-features = false, optional = true } packed_simd = { version = "0.3.0", features = ["into_bits"], optional = true } [features] diff --git a/src/backend/serial/scalar_mul/mod.rs b/src/backend/serial/scalar_mul/mod.rs index bec874b..8d859eb 100644 --- a/src/backend/serial/scalar_mul/mod.rs +++ b/src/backend/serial/scalar_mul/mod.rs @@ -26,3 +26,6 @@ pub mod straus; #[cfg(feature = "alloc")] pub mod precomputed_straus; + +#[cfg(feature = "alloc")] +pub mod pippenger; diff --git a/src/backend/serial/scalar_mul/pippenger.rs b/src/backend/serial/scalar_mul/pippenger.rs new file mode 100644 index 0000000..89ea723 --- /dev/null +++ b/src/backend/serial/scalar_mul/pippenger.rs @@ -0,0 +1,205 @@ +// -*- mode: rust; -*- +// +// This file is part of curve25519-dalek. +// Copyright (c) 2019 Oleg Andreev +// See LICENSE for licensing information. +// +// Authors: +// - Oleg Andreev + +//! Implementation of a variant of Pippenger's algorithm. + +#![allow(non_snake_case)] + +use core::borrow::Borrow; + +use edwards::EdwardsPoint; +use scalar::Scalar; +use traits::VartimeMultiscalarMul; + +#[allow(unused_imports)] +use prelude::*; + +/// Implements a version of Pippenger's algorithm. +/// +/// The algorithm works as follows: +/// +/// Let `n` be a number of point-scalar pairs. +/// Let `w` be a window of bits (6..8, chosen based on `n`, see cost factor). +/// +/// 1. Prepare `2^(w-1) - 1` buckets with indices `[1..2^(w-1))` initialized with identity points. +/// Bucket 0 is not needed as it would contain points multiplied by 0. +/// 2. Convert scalars to a radix-`2^w` representation with signed digits in `[-2^w/2, 2^w/2]`. +/// Note: only the last digit may equal `2^w/2`. +/// 3. Starting with the last window, for each point `i=[0..n)` add it to a a bucket indexed by +/// the point's scalar's value in the window. +/// 4. Once all points in a window are sorted into buckets, add buckets by multiplying each +/// by their index. Efficient way of doing it is to start with the last bucket and compute two sums: +/// intermediate sum from the last to the first, and the full sum made of all intermediate sums. +/// 5. Shift the resulting sum of buckets by `w` bits by using `w` doublings. +/// 6. Add to the return value. +/// 7. Repeat the loop. +/// +/// Approximate cost w/o wNAF optimizations (A = addition, D = doubling): +/// +/// ```ascii +/// cost = (n*A + 2*(2^w/2)*A + w*D + A)*256/w +/// | | | | | +/// | | | | looping over 256/w windows +/// | | | adding to the result +/// sorting points | shifting the sum by w bits (to the next window, starting from last window) +/// one by one | +/// into buckets adding/subtracting all buckets +/// multiplied by their indexes +/// using a sum of intermediate sums +/// ``` +/// +/// For large `n`, dominant factor is (n*256/w) additions. +/// However, if `w` is too big and `n` is not too big, then `(2^w/2)*A` could dominate. +/// Therefore, the optimal choice of `w` grows slowly as `n` grows. +/// +/// This algorithm is adapted from section 4 of https://eprint.iacr.org/2012/549.pdf. +pub struct Pippenger; + +#[cfg(any(feature = "alloc", feature = "std"))] +impl VartimeMultiscalarMul for Pippenger { + type Point = EdwardsPoint; + + fn optional_multiscalar_mul(scalars: I, points: J) -> Option + where + I: IntoIterator, + I::Item: Borrow, + J: IntoIterator>, + { + use traits::Identity; + + let mut scalars = scalars.into_iter(); + let size = scalars.by_ref().size_hint().0; + + // Digit width in bits. As digit width grows, + // number of point additions goes down, but amount of + // buckets and bucket additions grows exponentially. + let w = if size < 500 { + 6 + } else if size < 800 { + 7 + } else { + 8 + }; + + let max_digit: usize = 1 << w; + let digits_count: usize = (256 + w - 1) / w; // == ceil(256/w) + let buckets_count: usize = max_digit / 2; // digits are signed+centered hence 2^w/2, excluding 0-th bucket + + // Collect optimized scalars and points in buffers for repeated access + // (scanning the whole set per digit position). + let scalars = scalars + .into_iter() + .map(|s| s.borrow().to_radix_2w(w).0); + + let points = points + .into_iter() + .map(|p| p.map(|P| P.to_projective_niels())); + + let scalars_points = scalars.zip(points).map(|(s,maybe_p)| maybe_p.map(|p| (s,p) ) ) + .collect::>>(); + let scalars_points = match scalars_points { + Some(sp) => sp, + None => return None, + }; + + // Prepare 2^w/2 buckets. + // buckets[i] corresponds to a multiplication factor (i+1). + let mut buckets: Vec<_> = (0..buckets_count) + .map(|_| EdwardsPoint::identity()) + .collect(); + + let mut columns = (0..digits_count).rev().map(|digit_index| { + // Clear the buckets when processing another digit. + for i in 0..buckets_count { + buckets[i] = EdwardsPoint::identity(); + } + + // Iterate over pairs of (point, scalar) + // and add/sub the point to the corresponding bucket. + // Note: if we add support for precomputed lookup tables, + // we'll be adding/subtracting point premultiplied by `digits[i]` to buckets[0]. + for (digits, pt) in scalars_points.iter() { + let digit = digits[digit_index]; + if digit > 0 { + let b = (digit - 1) as usize; + buckets[b] = (&buckets[b] + pt).to_extended(); + } else if digit < 0 { + let b = (-digit - 1) as usize; + buckets[b] = (&buckets[b] - pt).to_extended(); + } + } + + // Add the buckets applying the multiplication factor to each bucket. + // The most efficient way to do that is to have a single sum with two running sums: + // an intermediate sum from last bucket to the first, and a sum of intermediate sums. + // + // For example, to add buckets 1*A, 2*B, 3*C we need to add these points: + // C + // C B + // C B A Sum = C + (C+B) + (C+B+A) + let mut buckets_intermediate_sum = buckets[buckets_count - 1]; + let mut buckets_sum = buckets[buckets_count - 1]; + for i in (0..(buckets_count - 1)).rev() { + buckets_intermediate_sum += buckets[i]; + buckets_sum += buckets_intermediate_sum; + } + + buckets_sum + }); + + // Take the high column as an initial value to avoid wasting time doubling the identity element in `fold()`. + // `unwrap()` always succeeds because we know we have more than zero digits. + let hi_column = columns.next().unwrap(); + + Some( + columns + .fold(hi_column, |total, p| total.mul_by_pow_2(w as u32) + p) + .into(), + ) + } +} + +#[cfg(test)] +mod test { + use super::*; + use constants; + use scalar::Scalar; + + #[test] + fn test_vartime_pippenger() { + // Reuse points across different tests + let mut n = 512; + let x = Scalar::from(2128506u64).invert(); + let y = Scalar::from(4443282u64).invert(); + let points: Vec<_> = (0..n) + .map(|i| constants::ED25519_BASEPOINT_POINT * Scalar::from(1 + i as u64)) + .collect(); + let scalars: Vec<_> = (0..n) + .map(|i| x + (Scalar::from(i as u64) * y)) // fast way to make ~random but deterministic scalars + .collect(); + + let premultiplied: Vec = scalars + .iter() + .zip(points.iter()) + .map(|(sc, pt)| sc * pt) + .collect(); + + while n > 0 { + let scalars = &scalars[0..n].to_vec(); + let points = &points[0..n].to_vec(); + let control: EdwardsPoint = premultiplied[0..n].iter().sum(); + + let subject = Pippenger::vartime_multiscalar_mul(scalars.clone(), points.clone()); + + assert_eq!(subject.compress(), control.compress()); + + n = n / 2; + } + } +} diff --git a/src/backend/serial/scalar_mul/variable_base.rs b/src/backend/serial/scalar_mul/variable_base.rs index 588a84e..2471375 100644 --- a/src/backend/serial/scalar_mul/variable_base.rs +++ b/src/backend/serial/scalar_mul/variable_base.rs @@ -3,7 +3,7 @@ use traits::Identity; use scalar::Scalar; use edwards::EdwardsPoint; -use backend::serial::curve_models::ProjectiveNielsPoint; +use backend::serial::curve_models::{ProjectiveNielsPoint, ProjectivePoint}; use window::LookupTable; /// Perform constant-time, variable-base scalar multiplication. @@ -23,10 +23,24 @@ pub(crate) fn mul(point: &EdwardsPoint, scalar: &Scalar) -> EdwardsPoint { // s*P = P*s_0 + 16*(P*s_1 + 16*(P*s_2 + 16*( ... + P*s_63)...)) // // We sum right-to-left. - let mut Q = EdwardsPoint::identity(); - for i in (0..64).rev() { - Q = Q.mul_by_pow_2(4); - Q = (&Q + &lookup_table.select(scalar_digits[i])).to_extended(); + + // Unwrap first loop iteration to save computing 16*identity + let mut tmp2 = ProjectivePoint::identity(); + let mut tmp3 = EdwardsPoint::identity(); + let mut tmp1 = &tmp3 + &lookup_table.select(scalar_digits[63]); + // Now tmp1 = s_63*P in P1xP1 coords + for i in (0..63).rev() { + tmp2 = tmp1.to_projective(); // tmp2 = (prev) in P2 coords + tmp1 = tmp2.double(); // tmp1 = 2*(prev) in P1xP1 coords + tmp2 = tmp1.to_projective(); // tmp2 = 2*(prev) in P2 coords + tmp1 = tmp2.double(); // tmp1 = 4*(prev) in P1xP1 coords + tmp2 = tmp1.to_projective(); // tmp2 = 4*(prev) in P2 coords + tmp1 = tmp2.double(); // tmp1 = 8*(prev) in P1xP1 coords + tmp2 = tmp1.to_projective(); // tmp2 = 8*(prev) in P2 coords + tmp1 = tmp2.double(); // tmp1 = 16*(prev) in P1xP1 coords + tmp3 = tmp1.to_extended(); // tmp3 = 16*(prev) in P3 coords + tmp1 = &tmp3 + &lookup_table.select(scalar_digits[i]); + // Now tmp1 = s_i*P + 16*(prev) in P1xP1 coords } - Q + tmp1.to_extended() } diff --git a/src/backend/vector/scalar_mul/mod.rs b/src/backend/vector/scalar_mul/mod.rs index 5c8734d..bdd6c4a 100644 --- a/src/backend/vector/scalar_mul/mod.rs +++ b/src/backend/vector/scalar_mul/mod.rs @@ -17,3 +17,6 @@ pub mod straus; #[cfg(feature = "alloc")] pub mod precomputed_straus; + +#[cfg(feature = "alloc")] +pub mod pippenger; diff --git a/src/backend/vector/scalar_mul/pippenger.rs b/src/backend/vector/scalar_mul/pippenger.rs new file mode 100644 index 0000000..0053e67 --- /dev/null +++ b/src/backend/vector/scalar_mul/pippenger.rs @@ -0,0 +1,165 @@ +// -*- mode: rust; -*- +// +// This file is part of curve25519-dalek. +// Copyright (c) 2019 Oleg Andreev +// See LICENSE for licensing information. +// +// Authors: +// - Oleg Andreev + +#![allow(non_snake_case)] + +use core::borrow::Borrow; + +use backend::vector::{CachedPoint, ExtendedPoint}; +use edwards::EdwardsPoint; +use scalar::Scalar; +use traits::{Identity, VartimeMultiscalarMul}; + +#[allow(unused_imports)] +use prelude::*; + +/// Implements a version of Pippenger's algorithm. +/// +/// See the documentation in the serial `scalar_mul::pippenger` module for details. +pub struct Pippenger; + +#[cfg(any(feature = "alloc", feature = "std"))] +impl VartimeMultiscalarMul for Pippenger { + type Point = EdwardsPoint; + + fn optional_multiscalar_mul(scalars: I, points: J) -> Option + where + I: IntoIterator, + I::Item: Borrow, + J: IntoIterator>, + { + let mut scalars = scalars.into_iter(); + let size = scalars.by_ref().size_hint().0; + let w = if size < 500 { + 6 + } else if size < 800 { + 7 + } else { + 8 + }; + + let max_digit: usize = 1 << w; + let digits_count: usize = (256 + w - 1) / w; // == ceil(256/w) + let buckets_count: usize = max_digit / 2; // digits are signed+centered hence 2^w/2, excluding 0-th bucket + + // Collect optimized scalars and points in a buffer for repeated access + // (scanning the whole collection per each digit position). + let scalars = scalars + .into_iter() + .map(|s| s.borrow().to_radix_2w(w).0); + + let points = points + .into_iter() + .map(|p| p.map(|P| CachedPoint::from(ExtendedPoint::from(P)))); + + let scalars_points = scalars.zip(points).map(|(s,maybe_p)| maybe_p.map(|p| (s,p) ) ) + .collect::>>(); + let scalars_points = match scalars_points { + Some(sp) => sp, + None => return None, + }; + + // Prepare 2^w/2 buckets. + // buckets[i] corresponds to a multiplication factor (i+1). + let mut buckets: Vec = (0..buckets_count) + .map(|_| ExtendedPoint::identity()) + .collect(); + + let mut columns = (0..digits_count).rev().map(|digit_index| { + // Clear the buckets when processing another digit. + for i in 0..buckets_count { + buckets[i] = ExtendedPoint::identity(); + } + + // Iterate over pairs of (point, scalar) + // and add/sub the point to the corresponding bucket. + // Note: if we add support for precomputed lookup tables, + // we'll be adding/subtractiong point premultiplied by `digits[i]` to buckets[0]. + for (digits, pt) in scalars_points.iter() { + let digit = digits[digit_index]; + if digit > 0 { + let b = (digit - 1) as usize; + buckets[b] = &buckets[b] + pt; + } else if digit < 0 { + let b = (-digit - 1) as usize; + buckets[b] = &buckets[b] - pt; + } + } + + // Add the buckets applying the multiplication factor to each bucket. + // The most efficient way to do that is to have a single sum with two running sums: + // an intermediate sum from last bucket to the first, and a sum of intermediate sums. + // + // For example, to add buckets 1*A, 2*B, 3*C we need to add these points: + // C + // C B + // C B A Sum = C + (C+B) + (C+B+A) + let mut buckets_intermediate_sum = buckets[buckets_count - 1]; + let mut buckets_sum = buckets[buckets_count - 1]; + for i in (0..(buckets_count - 1)).rev() { + buckets_intermediate_sum = + &buckets_intermediate_sum + &CachedPoint::from(buckets[i]); + buckets_sum = &buckets_sum + &CachedPoint::from(buckets_intermediate_sum); + } + + buckets_sum + }); + + // Take the high column as an initial value to avoid wasting time doubling the identity element in `fold()`. + // `unwrap()` always succeeds because we know we have more than zero digits. + let hi_column = columns.next().unwrap(); + + Some( + columns + .fold(hi_column, |total, p| { + &total.mul_by_pow_2(w as u32) + &CachedPoint::from(p) + }) + .into(), + ) + } +} + +#[cfg(test)] +mod test { + use super::*; + use constants; + use scalar::Scalar; + + #[test] + fn test_vartime_pippenger() { + // Reuse points across different tests + let mut n = 512; + let x = Scalar::from(2128506u64).invert(); + let y = Scalar::from(4443282u64).invert(); + let points: Vec<_> = (0..n) + .map(|i| constants::ED25519_BASEPOINT_POINT * Scalar::from(1 + i as u64)) + .collect(); + let scalars: Vec<_> = (0..n) + .map(|i| x + (Scalar::from(i as u64) * y)) // fast way to make ~random but deterministic scalars + .collect(); + + let premultiplied: Vec = scalars + .iter() + .zip(points.iter()) + .map(|(sc, pt)| sc * pt) + .collect(); + + while n > 0 { + let scalars = &scalars[0..n].to_vec(); + let points = &points[0..n].to_vec(); + let control: EdwardsPoint = premultiplied[0..n].iter().sum(); + + let subject = Pippenger::vartime_multiscalar_mul(scalars.clone(), points.clone()); + + assert_eq!(subject.compress(), control.compress()); + + n = n / 2; + } + } +} diff --git a/src/edwards.rs b/src/edwards.rs index 10b4614..e43adb4 100644 --- a/src/edwards.rs +++ b/src/edwards.rs @@ -394,7 +394,14 @@ impl ConditionallySelectable for EdwardsPoint { impl ConstantTimeEq for EdwardsPoint { fn ct_eq(&self, other: &EdwardsPoint) -> Choice { - self.compress().ct_eq(&other.compress()) + // We would like to check that the point (X/Z, Y/Z) is equal to + // the point (X'/Z', Y'/Z') without converting into affine + // coordinates (x, y) and (x', y'), which requires two inversions. + // We have that X = xZ and X' = x'Z'. Thus, x = x' is equivalent to + // (xZ)Z' = (x'Z')Z, and similarly for the y-coordinate. + + (&self.X * &other.Z).ct_eq(&(&other.X * &self.Z)) + & (&self.Y * &other.Z).ct_eq(&(&other.Y * &self.Z)) } } @@ -669,11 +676,15 @@ impl VartimeMultiscalarMul for EdwardsPoint { assert_eq!(s_hi, Some(s_lo)); assert_eq!(p_hi, Some(p_lo)); - // Now we know there's a single size. When we do - // size-dependent algorithm dispatch, use this as the hint. - let _size = s_lo; + // Now we know there's a single size. + // Use this as the hint to decide which algorithm to use. + let size = s_lo; - scalar_mul::straus::Straus::optional_multiscalar_mul(scalars, points) + if size < 190 { + scalar_mul::straus::Straus::optional_multiscalar_mul(scalars, points) + } else { + scalar_mul::pippenger::Pippenger::optional_multiscalar_mul(scalars, points) + } } } diff --git a/src/scalar.rs b/src/scalar.rs index a7a8310..00ebad8 100644 --- a/src/scalar.rs +++ b/src/scalar.rs @@ -961,6 +961,80 @@ impl Scalar { output } + /// Creates a representation of a Scalar in radix 64, 128 or 256 for use with the Pippenger algorithm. + /// For lower radix, use `to_radix_16`, which is used by the Straus multi-scalar multiplication. + /// Higher radixes are not supported to save cache space. Radix 256 is near-optimal even for very + /// large inputs. + /// + /// Radix below 64 or above 256 is prohibited. + /// This method returns digits in a fixed-sized array, excess digits are zeroes. + /// The second returned value is the number of digits. + /// + /// ## Scalar representation + /// + /// Radix \\(2\^w\\), with \\(n = ceil(256/w)\\) coefficients in \\([-(2\^w)/2,(2\^w)/2)\\), + /// i.e., scalar is represented using digits \\(a\_i\\) such that + /// $$ + /// a = a\_0 + a\_1 2\^1w + \cdots + a_{n-1} 2\^{w*(n-1)}, + /// $$ + /// with \\(-2\^w/2 \leq a_i < 2\^w/2\\) for \\(0 \leq i < (n-1)\\) and \\(-2\^w/2 \leq a_{n-1} \leq 2\^w/2\\). + /// + pub(crate) fn to_radix_2w(&self, w: usize) -> ([i8; 43], usize) { + debug_assert!(w >= 6); + debug_assert!(w <= 8); + + let digits_count = (256 + w - 1)/w as usize; + debug_assert!(digits_count <= 43); + + use byteorder::{ByteOrder, LittleEndian}; + + // Scalar formatted as four `u64`s with carry bit packed into the highest bit. + let mut scalar64x4 = [0u64; 4]; + LittleEndian::read_u64_into(&self.bytes, &mut scalar64x4[0..4]); + + let radix: u64 = 1 << w; + let window_mask: u64 = radix - 1; + + let mut carry = 0u64; + let mut digits = [0i8; 43]; + for i in 0..digits_count { + // Construct a buffer of bits of the scalar, starting at `bit_offset`. + let bit_offset = i*w; + let u64_idx = bit_offset / 64; + let bit_idx = bit_offset % 64; + + // Read the bits from the scalar + let bit_buf: u64; + if bit_idx < 64 - w || u64_idx == 3 { + // This window's bits are contained in a single u64, + // or it's the last u64 anyway. + bit_buf = scalar64x4[u64_idx] >> bit_idx; + } else { + // Combine the current u64's bits with the bits from the next u64 + bit_buf = (scalar64x4[u64_idx] >> bit_idx) | (scalar64x4[1+u64_idx] << (64 - bit_idx)); + } + + // Read the actual coefficient value from the window + let coef = carry + (bit_buf & window_mask); // coef = [0, 2^r) + + // Recenter coefficients from [0,2^r) to [-2^r/2, 2^r/2) + carry = (coef + (radix/2) as u64) >> w; + digits[i] = ((coef as i64) - (carry << w) as i64) as i8; + } + + // Apply the resulting carry to the last digit + // Since the highest bit of the 256-bit integer is 0, + // the last coefficient would always be in the lower half _inclusive_, + // so the carry in the end can be 1 iff the word equals 2^r/2. + // Since ±2^r/2 values are valid, to avoid adding an extra word, + // we allow the last word to touch the value 2^r/2. + // XXX: make sure tests cover this case, so the carry is non-zero and this line matters. + // Maybe it never happens to be non-zero for r=6/7/8?... + digits[digits_count-1] += (carry << w) as i8; + + (digits, digits_count) + } + /// Unpack this `Scalar` to an `UnpackedScalar` for faster arithmetic. pub(crate) fn unpack(&self) -> UnpackedScalar { UnpackedScalar::from_bytes(&self.bytes) @@ -1437,4 +1511,33 @@ mod test { assert_eq!(a * b, Scalar::one()); } } + + #[test] + fn test_pippenger_radix() { + use core::iter; + // For each valid radix it tests that 1000 random-ish scalars can be restored + // from the produced representation precisely. + for w in 6..9 { + for scalar in (2..100).map(|s| Scalar::from(s as u64).invert() ).chain(iter::once(-Scalar::one())) { + let (digits, digits_count) = scalar.to_radix_2w(w); + + let radix = Scalar::from((1<(scalars: I, points: J) -> Self::Point where I: IntoIterator,