diff --git a/.travis.yml b/.travis.yml index cffc1c2..4a28c7e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -6,10 +6,11 @@ rust: - nightly env: - - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='yolocrypto' - - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='yolocrypto serde' + - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='' + - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='serde' + - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='nightly' - TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='yolocrypto nightly' - - TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='yolocrypto bench' + - TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='nightly bench' - TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='yolocrypto nightly bench' - TEST_COMMAND=build EXTRA_FLAGS=--no-default-features FEATURES='' @@ -20,14 +21,18 @@ matrix: # run benchmarks, which causes dalek not to build on stable. See # https://github.com/isislovecruft/curve25519-dalek/pull/38#issuecomment-286027562 - rust: stable - env: TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='yolocrypto bench' + env: TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='nightly bench' - rust: beta - env: TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='yolocrypto bench' + env: TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='nightly bench' - rust: stable env: TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='yolocrypto nightly bench' - rust: beta env: TEST_COMMAND=bench EXTRA_FLAGS='' FEATURES='yolocrypto nightly bench' # Test nightly features, such as radix_51, only on nightly. + - rust: stable + env: TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='nightly' + - rust: beta + env: TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='nightly' - rust: stable env: TEST_COMMAND=test EXTRA_FLAGS='' FEATURES='yolocrypto nightly' - rust: beta diff --git a/Cargo.toml b/Cargo.toml index aa3fd59..dcb4747 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -24,6 +24,10 @@ rustdoc-args = ["--html-in-header", ".cargo/registry/src/github.com-1ecc6299db9e [badges] travis-ci = { repository = "isislovecruft/curve25519-dalek", branch = "master"} +[dependencies.stdsimd] +git = "https://github.com/rust-lang-nursery/stdsimd" +optional = true + [dependencies.serde] version = "1.0" optional = true @@ -60,6 +64,10 @@ digest = "0.7" arrayref = "0.3.4" clear_on_drop = "=0.2.3" +[build-dependencies.stdsimd] +git = "https://github.com/rust-lang-nursery/stdsimd" +optional = true + [build-dependencies.serde] version = "1.0" optional = true @@ -69,11 +77,11 @@ nightly = ["radix_51", "subtle/nightly", "clear_on_drop/nightly"] default = ["std"] std = ["rand", "subtle/std"] alloc = [] -# This isn't used at the moment, but keep it around for future yolocrypto features. -yolocrypto = [] +yolocrypto = ["avx2_backend"] bench = [] # Radix-51 arithmetic using u128 radix_51 = [] # Include precomputed basepoint tables. This is off by default so that build.rs can generate the tables, and then re-enabled by build.rs in the main-stage compilation. precomputed_tables = [] - +# experimental avx2 support +avx2_backend = ["nightly", "stdsimd"] diff --git a/build.rs b/build.rs index 28044c3..2c30435 100644 --- a/build.rs +++ b/build.rs @@ -1,4 +1,5 @@ #![cfg_attr(feature = "nightly", feature(i128_type))] +#![cfg_attr(feature = "nightly", feature(cfg_target_feature))] #![allow(unused_variables)] #![allow(non_snake_case)] #![allow(dead_code)] @@ -22,6 +23,8 @@ use std::path::Path; // For instance, this shouldn't exist here at all, but it does. #[cfg(feature = "serde")] extern crate serde; +#[cfg(feature = "yolocrypto")] +extern crate stdsimd; // Public modules diff --git a/src/backend/avx2/constants.rs b/src/backend/avx2/constants.rs new file mode 100644 index 0000000..087d96c --- /dev/null +++ b/src/backend/avx2/constants.rs @@ -0,0 +1,112 @@ +// -*- mode: rust; -*- +// +// This file is part of curve25519-dalek. +// Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence +// See LICENSE for licensing information. +// +// Authors: +// - Isis Agora Lovecruft +// - Henry de Valence + +//! This module contains constants used by the AVX2 backend. + +use stdsimd::simd::u32x8; + +use backend::avx2::field::FieldElement32x4; +use backend::avx2::edwards::ExtendedPoint; + +/// The low limbs of (2p, 2p, 2p, 2p), so that +/// ```no_run +/// (2p, 2p, 2p, 2p) = [P_TIMES_2_LO, P_TIMES_2_HI, P_TIMES_2_HI, P_TIMES_2_HI, P_TIMES_2_HI] +/// ``` +pub(crate) static P_TIMES_2_LO: u32x8 = + u32x8::new(67108845 << 1, 67108845 << 1, 33554431 << 1, 33554431 << 1, 67108845 << 1, 67108845 << 1, 33554431 << 1, 33554431 << 1); + +/// The high limbs of (2p, 2p, 2p, 2p), so that +/// ```no_run +/// (2p, 2p, 2p, 2p) = [P_TIMES_2_LO, P_TIMES_2_HI, P_TIMES_2_HI, P_TIMES_2_HI, P_TIMES_2_HI] +/// ``` +pub(crate) static P_TIMES_2_HI: u32x8 = + u32x8::new(67108863 << 1, 67108863 << 1, 33554431 << 1, 33554431 << 1, 67108863 << 1, 67108863 << 1, 33554431 << 1, 33554431 << 1); + +/// The low limbs of (16p, 16p, 16p, 16p), so that +/// ```no_run +/// (16p, 16p, 16p, 16p) = [P_TIMES_16_LO, P_TIMES_16_HI, P_TIMES_16_HI, P_TIMES_16_HI, P_TIMES_16_HI] +/// ``` +pub(crate) static P_TIMES_16_LO: u32x8 = + u32x8::new(67108845 << 4, 67108845 << 4, 33554431 << 4, 33554431 << 4, 67108845 << 4, 67108845 << 4, 33554431 << 4, 33554431 << 4); + +/// The high limbs of (16p, 16p, 16p, 16p), so that +/// ```no_run +/// (16p, 16p, 16p, 16p) = [P_TIMES_16_LO, P_TIMES_16_HI, P_TIMES_16_HI, P_TIMES_16_HI, P_TIMES_16_HI] +/// ``` +pub(crate) static P_TIMES_16_HI: u32x8 = + u32x8::new(67108863 << 4, 67108863 << 4, 33554431 << 4, 33554431 << 4, 67108863 << 4, 67108863 << 4, 33554431 << 4, 33554431 << 4); + +pub(crate) static P_TIMES_2_MASKED: FieldElement32x4 = FieldElement32x4([ + u32x8::new( 0, 134217690, 0, 67108862, 134217690, 0, 67108862, 0), + u32x8::new( 0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), + u32x8::new( 0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), + u32x8::new( 0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0), + u32x8::new( 0, 134217726, 0, 67108862, 134217726, 0, 67108862, 0) +]); + +/// Odd multiples of the Ed25519 basepoint: +pub static ODD_MULTIPLES_OF_BASEPOINT: [ExtendedPoint; 8] = [ + ExtendedPoint(FieldElement32x4([ + u32x8::new(52811034, 40265304, 25909283, 26843545, 1, 28827043, 0, 27438313), + u32x8::new(16144682, 13421772, 17082669, 20132659, 0, 39759291, 0, 244362), + u32x8::new(27570973, 26843545, 30858332, 6710886, 0, 8635006, 0, 11264893), + u32x8::new(40966398, 53687091, 8378388, 13421772, 0, 19351346, 0, 13413597), + u32x8::new(20764389, 40265318, 8758491, 26843545, 0, 16611511, 0, 27139452), + ])), + ExtendedPoint(FieldElement32x4([ + u32x8::new(63703867, 19156774, 608100, 2486757, 12685460, 3173753, 21649412, 16313381), + u32x8::new(52397038, 65858675, 26775664, 16661035, 14269998, 9080558, 1059463, 28938752), + u32x8::new( 5461635, 28034025, 23358301, 1245198, 1367765, 20288887, 31111942, 18395221), + u32x8::new( 1886934, 32436996, 681756, 18977693, 8129860, 40112764, 25764567, 11876840), + u32x8::new(63042604, 52399761, 22087481, 29829870, 8565820, 33723612, 28645162, 8502864), + ])), + ExtendedPoint(FieldElement32x4([ + u32x8::new(14879397, 3951036, 9454671, 16606238, 23529732, 44147004, 11890541, 17067526), + u32x8::new(58509479, 57216664, 9671992, 32001147, 60966207, 11801823, 10808378, 15115613), + u32x8::new(54854992, 39210911, 8112050, 1353604, 1337416, 35520540, 32967851, 17786030), + u32x8::new(59007462, 40864509, 26240923, 30403852, 28456403, 21546582, 32732450, 21005910), + u32x8::new(40711675, 22446613, 9664668, 12483629, 26142305, 56254715, 15439904, 214849), + ])), + ExtendedPoint(FieldElement32x4([ + u32x8::new(52231579, 51632644, 173613, 7677257, 26374424, 45994428, 5303371, 1425942), + u32x8::new(38126791, 48854506, 23252518, 30611978, 49977504, 66706952, 1076178, 27100873), + u32x8::new(26349427, 63077566, 20258199, 3884787, 33226507, 2371423, 5787271, 18628170), + u32x8::new(15005754, 22729577, 4978944, 2522289, 1404784, 56367795, 22517039, 29271243), + u32x8::new(22748934, 35977548, 25561257, 31734126, 22775284, 32000077, 927866, 2278697), + ])), + ExtendedPoint(FieldElement32x4([ + u32x8::new(66090281, 61980626, 23780289, 6519561, 62542590, 47174086, 28818882, 15661068), + u32x8::new(17433715, 12931425, 12232056, 7885877, 44179512, 35590146, 32787344, 22631048), + u32x8::new(43729883, 6870635, 15782399, 11810556, 2652935, 31800505, 23683367, 13638649), + u32x8::new(64007953, 40242373, 32810277, 20180235, 20399465, 48133835, 32913956, 19094667), + u32x8::new(56562708, 40269142, 18953105, 9027935, 35700921, 12896915, 14757156, 22773619), + ])), + ExtendedPoint(FieldElement32x4([ + u32x8::new(65129016, 34709402, 25132940, 13788431, 3661652, 16914498, 27409409, 18941039), + u32x8::new(42488074, 49427602, 6177212, 20812339, 41644653, 2977316, 12162542, 5293661), + u32x8::new( 7981168, 12223605, 6239200, 20403609, 20710415, 4828170, 11627702, 4431044), + u32x8::new(65817142, 96824, 25021652, 16364722, 50410869, 24651857, 6979034, 33176209), + u32x8::new(33008344, 8687253, 27859668, 28796356, 30192014, 11975680, 11991047, 27710707), + ])), + ExtendedPoint(FieldElement32x4([ + u32x8::new(14676653, 50945941, 13489249, 31456262, 47726639, 21761847, 3324839, 7843947), + u32x8::new(53352326, 8688989, 12944061, 12994004, 50113821, 37990636, 1537898, 20483689), + u32x8::new(46786852, 15572264, 24004728, 7566233, 32596174, 34437796, 23201722, 3431551), + u32x8::new(49025674, 52497128, 13273618, 10266201, 66795206, 2887684, 30966565, 33449990), + u32x8::new(53210238, 65839385, 15458877, 18409918, 24777464, 25586795, 15335748, 12323382), + ])), + ExtendedPoint(FieldElement32x4([ + u32x8::new(57816016, 23106045, 24948505, 27413507, 32551424, 26145165, 22632568, 27527446), + u32x8::new(53022711, 40974949, 14110533, 30646997, 51399118, 53289754, 32528560, 15822835), + u32x8::new(23810949, 51779690, 17532625, 21326637, 60314333, 43761996, 4852905, 3474945), + u32x8::new(13323962, 10752742, 16431634, 26425049, 24258356, 53260846, 19756601, 19546842), + u32x8::new(17403634, 52199608, 32323720, 5313255, 48522162, 33376516, 31903659, 15291466), + ])), +]; diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs new file mode 100644 index 0000000..d995859 --- /dev/null +++ b/src/backend/avx2/edwards.rs @@ -0,0 +1,1057 @@ +// -*- mode: rust; -*- +// +// This file is part of curve25519-dalek. +// Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence +// See LICENSE for licensing information. +// +// Authors: +// - Isis Agora Lovecruft +// - Henry de Valence + +//! Extended Twisted Edwards for Curve25519, using AVX2. + +// just going to own it +#![allow(bad_style)] + +use std::convert::From; +use std::ops::{Index, Add, Sub, Mul, Neg}; + +use stdsimd::simd::{u32x8, i32x8}; + +use subtle::ConditionallyAssignable; + +use edwards; +use scalar::Scalar; +use curve_models::window::LookupTable; + +use traits::Identity; + +use backend::avx2::field::FieldElement32x4; + +use backend::avx2::field::{A_LANES, B_LANES, C_LANES, D_LANES, ALL_LANES}; +use backend::avx2::field::D_LANES64; + +use backend::avx2; + +/// A point on Curve25519, represented in an AVX2-friendly format. +#[derive(Copy, Clone, Debug)] +pub struct ExtendedPoint(pub(super) FieldElement32x4); + +impl From for ExtendedPoint { + fn from(P: edwards::ExtendedPoint) -> ExtendedPoint { + ExtendedPoint(FieldElement32x4::new(&P.X, &P.Y, &P.Z, &P.T)) + } +} + +impl From for edwards::ExtendedPoint { + fn from(P: ExtendedPoint) -> edwards::ExtendedPoint { + let tmp = P.0.split(); + edwards::ExtendedPoint{X: tmp[0], Y: tmp[1], Z: tmp[2], T: tmp[3]} + } +} + +impl ConditionallyAssignable for ExtendedPoint { + fn conditional_assign(&mut self, other: &ExtendedPoint, choice: u8) { + self.0.conditional_assign(&other.0, choice); + } +} + +impl Default for ExtendedPoint { + fn default() -> ExtendedPoint { + ExtendedPoint::identity() + } +} + +impl Identity for ExtendedPoint { + fn identity() -> ExtendedPoint { + ExtendedPoint(FieldElement32x4([ + u32x8::new(0,1,0,0,1,0,0,0), + u32x8::splat(0), + u32x8::splat(0), + u32x8::splat(0), + u32x8::splat(0), + ])) + } +} + +/// A cached point with some precomputed variables used for readdition. +#[derive(Copy, Clone, Debug)] +pub struct CachedPoint(pub(super) FieldElement32x4); + +impl From for CachedPoint { + fn from(P: ExtendedPoint) -> CachedPoint { + let mut x = P.0; + + // x = (S2 S3 Z2 T2) + x.diff_sum(0b00001111); + + // x = (121666*S2 121666*S3 2*121666*Z2 2*121665*T2) + x.scale_by_curve_constants(); + + // x = (121666*S2 121666*S3 2*121666*Z2 -2*121665*T2) + x.negate(D_LANES); + + CachedPoint(x) + } +} + +impl Default for CachedPoint { + fn default() -> CachedPoint { + CachedPoint::identity() + } +} + +impl Identity for CachedPoint { + fn identity() -> CachedPoint { + CachedPoint(FieldElement32x4([ + u32x8::new(121647, 121666, 0, 0, 243332, 67108845, 0, 33554431), + u32x8::new(67108864, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), + u32x8::new(67108863, 0, 33554431, 0, 0, 67108863, 0, 33554431), + ])) + } +} + +impl ConditionallyAssignable for CachedPoint { + fn conditional_assign(&mut self, other: &CachedPoint, choice: u8) { + self.0.conditional_assign(&other.0, choice); + } +} + +impl<'a> Neg for &'a CachedPoint { + type Output = CachedPoint; + + fn neg(self) -> CachedPoint { + let mut neg = *self; + neg.0.swap_AB(); + neg.0.negate_lazy(D_LANES); + neg + } +} + +impl<'a> Neg for &'a ExtendedPoint { + type Output = ExtendedPoint; + + fn neg(self) -> ExtendedPoint { + let mut neg = *self; + // (X Y Z T) -> (-X Y Z -T) + neg.0.negate(A_LANES | D_LANES); + neg + } +} + +impl ExtendedPoint { + fn double(&self) -> ExtendedPoint { + unsafe { + use stdsimd::vendor::_mm256_permute2x128_si256; + use stdsimd::vendor::_mm256_permutevar8x32_epi32; + use stdsimd::vendor::_mm256_blend_epi32; + use stdsimd::vendor::_mm256_shuffle_epi32; + + let P = &self.0; + + let mut t0 = FieldElement32x4::zero(); + let mut t1 = FieldElement32x4::zero(); + + // Want to compute (X1 Y1 Z1 X1+Y1). + // Not sure how to do this less expensively than computing + // (X1 Y1 Z1 T1) --(256bit shuffle)--> (X1 Y1 X1 Y1) + // (X1 Y1 X1 Y1) --(2x128b shuffle)--> (Y1 X1 Y1 X1) + // and then adding. + + // Set t0 = (X1 Y1 X1 Y1) + t0.0[0] = _mm256_permute2x128_si256(P.0[0].into(), P.0[0].into(), 0b0000_0000).into(); + t0.0[1] = _mm256_permute2x128_si256(P.0[1].into(), P.0[1].into(), 0b0000_0000).into(); + t0.0[2] = _mm256_permute2x128_si256(P.0[2].into(), P.0[2].into(), 0b0000_0000).into(); + t0.0[3] = _mm256_permute2x128_si256(P.0[3].into(), P.0[3].into(), 0b0000_0000).into(); + t0.0[4] = _mm256_permute2x128_si256(P.0[4].into(), P.0[4].into(), 0b0000_0000).into(); + + // Set t1 = (Y1 X1 Y1 X1) + t1.0[0] = _mm256_shuffle_epi32(t0.0[0].into(), 0b10_11_00_01).into(); + t1.0[1] = _mm256_shuffle_epi32(t0.0[1].into(), 0b10_11_00_01).into(); + t1.0[2] = _mm256_shuffle_epi32(t0.0[2].into(), 0b10_11_00_01).into(); + t1.0[3] = _mm256_shuffle_epi32(t0.0[3].into(), 0b10_11_00_01).into(); + t1.0[4] = _mm256_shuffle_epi32(t0.0[4].into(), 0b10_11_00_01).into(); + + // Set t0 = (X1+Y1 X1+Y1 X1+Y1 X1+Y1) + t0.0[0] = t0.0[0] + t1.0[0]; + t0.0[1] = t0.0[1] + t1.0[1]; + t0.0[2] = t0.0[2] + t1.0[2]; + t0.0[3] = t0.0[3] + t1.0[3]; + t0.0[4] = t0.0[4] + t1.0[4]; + + // Set t0 = (X1 Y1 Z1 X1+Y1) + // why does this intrinsic take an i32 for the imm8 ??? + t0.0[0] = _mm256_blend_epi32(P.0[0].into(), t0.0[0].into(), D_LANES as i32).into(); + t0.0[1] = _mm256_blend_epi32(P.0[1].into(), t0.0[1].into(), D_LANES as i32).into(); + t0.0[2] = _mm256_blend_epi32(P.0[2].into(), t0.0[2].into(), D_LANES as i32).into(); + t0.0[3] = _mm256_blend_epi32(P.0[3].into(), t0.0[3].into(), D_LANES as i32).into(); + t0.0[4] = _mm256_blend_epi32(P.0[4].into(), t0.0[4].into(), D_LANES as i32).into(); + + // Set t1 = t0^2, negating the D values + t1 = t0.square(D_LANES64); + + // Now t1 = (S1 S2 S3 -S4) + + let c0 = u32x8::new(0,0,2,2,0,0,2,2); // (ABCD) -> (AAAA) + let c1 = u32x8::new(1,1,3,3,1,1,3,3); // (ABCD) -> (BBBB) + + // See discussion of bounds in the module-level documentation. + // + // We want to compute + // + // + | S1 | S1 | S1 | S1 | + // + | S2 | | | S2 | + // + | | | S3 | | + // + | | | S3 | | + // + | | | |-S4 | + // + | | 2p | 2p | | + // - | | S2 | S2 | | + // ======================= + // S5 S6 S8 S9 + // + for i in 0..5 { + let zero = i32x8::splat(0); + let S1 = _mm256_permutevar8x32_epi32(t1.0[i], c0); + let S2 = _mm256_permutevar8x32_epi32(t1.0[i], c1); + let S3_2: u32x8 = _mm256_blend_epi32(zero, (t1.0[i] + t1.0[i]).into(), 0b01010000).into(); + // tmp0 = (0 0 2*S3 -S4) + let tmp0: u32x8 = _mm256_blend_epi32(S3_2.into(), t1.0[i].into(), 0b10100000).into(); + t0.0[i] = (avx2::constants::P_TIMES_2_MASKED.0[i] + tmp0) + S1; + t0.0[i] = t0.0[i] + _mm256_blend_epi32(zero, S2.into(), 0b10100101).into(); + t0.0[i] = t0.0[i] - _mm256_blend_epi32(S2.into(), zero, 0b10100101).into(); + } + + let c0 = u32x8::new(4,0,6,2,4,0,6,2); // (ABCD) -> (CACA) + let c1 = u32x8::new(5,1,7,3,1,5,3,7); // (ABCD) -> (DBBD) + + for i in 0..5 { + let tmp = t0.0[i]; + t0.0[i] = _mm256_permutevar8x32_epi32(tmp, c0); + t1.0[i] = _mm256_permutevar8x32_epi32(tmp, c1); + } + + ExtendedPoint(&t0 * &t1) + } + } + + pub fn mult_by_pow_2(&self, k: u32) -> ExtendedPoint { + let mut tmp: ExtendedPoint = *self; + for _ in 0..k { + tmp = tmp.double(); + } + tmp + } +} + +impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint { + type Output = ExtendedPoint; + + /// Uses a slight tweak of the parallel unified formulas of HWCD'08 + fn add(self, other: &'b CachedPoint) -> ExtendedPoint { + unsafe { + use stdsimd::vendor::_mm256_permutevar8x32_epi32; + + let mut tmp = self.0; + + // tmp = (Y1-X1 Y1+X1 Z1 T1) = (S0 S1 Z1 T1) + tmp.diff_sum(A_LANES | B_LANES); + + // tmp = (S0*S2' S1*S3' Z1*Z2' T1*T2') = (S8 S9 S10 S11) + tmp = &tmp * &other.0; + + // tmp = (S8 S9 S11 S10) + tmp.swap_CD(); + + // tmp = (S9-S8 S9+S8 S10-S11 S10+S11) = (S12 S13 S14 S15) + tmp.diff_sum(ALL_LANES); + + let c0 = u32x8::new(0,5,2,7,5,0,7,2); // (ABCD) -> (ADDA) + let c1 = u32x8::new(4,1,6,3,4,1,6,3); // (ABCD) -> (CBCB) + + // set t0 = (S12 S15 S15 S12) + // set t1 = (S14 S13 S14 S13) + let mut t0 = FieldElement32x4::zero(); + let mut t1 = FieldElement32x4::zero(); + for i in 0..5 { + t0.0[i] = _mm256_permutevar8x32_epi32(tmp.0[i], c0); + t1.0[i] = _mm256_permutevar8x32_epi32(tmp.0[i], c1); + } + + // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) + ExtendedPoint(&t0 * &t1) + } + } +} + +impl<'a, 'b> Add<&'b ExtendedPoint> for &'a ExtendedPoint { + type Output = ExtendedPoint; + + /// Uses a slight tweak of the parallel unified formulas of HWCD'08 + fn add(self, other: &'b ExtendedPoint) -> ExtendedPoint { + unsafe { + use stdsimd::vendor::_mm256_permute2x128_si256; + use stdsimd::vendor::_mm256_permutevar8x32_epi32; + use stdsimd::vendor::_mm256_blend_epi32; + + let P: &FieldElement32x4 = &self.0; + let Q: &FieldElement32x4 = &other.0; + + let mut t0 = FieldElement32x4::zero(); + let mut t1 = FieldElement32x4::zero(); + + // set t0 = (X1 Y1 X2 Y2) + for i in 0..5 { + t0.0[i] = _mm256_permute2x128_si256(P.0[i].into(), Q.0[i].into(), 32).into(); + } + + // set t0 = (Y1-X1 Y1+X1 Y2-X2 Y2+X2) = (S0 S1 S2 S3) + t0.diff_sum(0xff); + + // set t1 = (S0 S1 Z1 T1) + // set t0 = (S2 S3 Z2 T2) + for i in 0..5 { + // why does this intrinsic take an i32 for the imm8 ??? + t1.0[i] = _mm256_blend_epi32(t0.0[i].into(), P.0[i].into(), (C_LANES | D_LANES) as i32).into(); + t0.0[i] = _mm256_permute2x128_si256(t0.0[i].into(), Q.0[i].into(), 49).into(); + } + + // set t2 = (S0*S2 S1*S3 Z1*Z2 T1*T2) = (S4 S5 S6 S7) + let mut t2 = &t0 * &t1; + + //// set t2 = (S8 S9 S10 S11) + // set t2 = (121666*S4 121666*S5 2*121666*S6 2*121665*S7) + // = ( S8 S9 S10 -S11) + t2.scale_by_curve_constants(); + + // set t2 = (S8 S9 -S11 S10) + t2.swap_CD(); + + // set t2 = (S9-S8 S9+S8 S10+S11 S10-S11) = (S12 S13 S15 S14) + t2.diff_sum(0xff); + + let c0 = u32x8::new(0,4,2,6,4,0,6,2); // (ABCD) -> (ACCA) + let c1 = u32x8::new(5,1,7,3,5,1,7,3); // (ABCD) -> (DBDB) + + // set t0 = (S12 S15 S15 S12) + // set t1 = (S14 S13 S14 S13) + for i in 0..5 { + t0.0[i] = _mm256_permutevar8x32_epi32(t2.0[i], c0); + t1.0[i] = _mm256_permutevar8x32_epi32(t2.0[i], c1); + } + + // return (S12*S14 S15*S13 S15*S14 S12*S13) = (X3 Y3 Z3 T3) + ExtendedPoint(&t0 * &t1) + } + } +} + +impl<'a, 'b> Sub<&'b ExtendedPoint> for &'a ExtendedPoint { + type Output = ExtendedPoint; + + /// Implement subtraction by negating the point and adding. + /// + /// Empirically, this seems about the same cost as a custom subtraction impl (maybe because the + /// benefit is cancelled by increased code size?) + fn sub(self, other: &'b ExtendedPoint) -> ExtendedPoint { + self + &(-other) + } +} + +impl From for LookupTable { + fn from(P: ExtendedPoint) -> Self { + let mut points = [CachedPoint::from(P); 8]; + for i in 0..7 { + points[i+1] = (&P + &points[i]).into(); + } + LookupTable(points) + } +} + +impl<'a, 'b> Mul<&'b Scalar> for &'a ExtendedPoint { + type Output = ExtendedPoint; + /// Scalar multiplication: compute `scalar * self`. + /// + /// Uses a window of size 4. + fn mul(self, scalar: &'b Scalar) -> ExtendedPoint { + // Construct a lookup table of [P,2P,3P,4P,5P,6P,7P,8P] + let lookup_table = LookupTable::::from(*self); + + // Setting s = scalar, compute + // + // s = s_0 + s_1*16^1 + ... + s_63*16^63, + // + // with `-8 ≤ s_i < 8` for `0 ≤ i < 63` and `-8 ≤ s_63 ≤ 8`. + let scalar_digits = scalar.to_radix_16(); + + // Compute s*P as + // + // s*P = P*(s_0 + s_1*16^1 + s_2*16^2 + ... + s_63*16^63) + // s*P = P*s_0 + P*s_1*16^1 + P*s_2*16^2 + ... + P*s_63*16^63 + // 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 = ExtendedPoint::identity(); + for i in (0..64).rev() { + // Q = 16*Q + Q = Q.mult_by_pow_2(4); + // Q += P*s_i + Q = &Q + &lookup_table.select(scalar_digits[i]); + } + Q + } +} + +#[derive(Clone)] +pub struct EdwardsBasepointTable(pub [LookupTable; 32]); + +impl<'a, 'b> Mul<&'b Scalar> for &'a EdwardsBasepointTable { + type Output = ExtendedPoint; + + fn mul(self, scalar: &'b Scalar) -> ExtendedPoint { + let a = scalar.to_radix_16(); + + let tables = &self.0; + let mut P = ExtendedPoint::identity(); + + for i in (0..64).filter(|x| x % 2 == 1) { + P = &P + &tables[i/2].select(a[i]); + } + + P = P.mult_by_pow_2(4); + + for i in (0..64).filter(|x| x % 2 == 0) { + P = &P + &tables[i/2].select(a[i]); + } + + P + } +} + +impl<'a, 'b> Mul<&'a EdwardsBasepointTable> for &'b Scalar { + type Output = ExtendedPoint; + + /// Given `self` a table of precomputed multiples of the point `B`, compute `B * s`. + fn mul(self, basepoint_table: &'a EdwardsBasepointTable) -> ExtendedPoint { + basepoint_table * &self + } +} + +impl EdwardsBasepointTable { + /// Create a table of precomputed multiples of `basepoint`. + pub fn create(basepoint: &ExtendedPoint) -> EdwardsBasepointTable { + // XXX use init_with + let mut table = EdwardsBasepointTable([LookupTable::default(); 32]); + let mut P = *basepoint; + for i in 0..32 { + // P = (16^2)^i * B + table.0[i] = LookupTable::from(P); + P = P.mult_by_pow_2(8); + } + table + } +} + +/// Given a vector of (possibly secret) scalars and a vector of +/// (possibly secret) points, compute `c_1 P_1 + ... + c_n P_n`. +/// +/// This function has the same behaviour as +/// `vartime::multiscalar_mult` but is constant-time. +/// +/// # Input +/// +/// A vector of `Scalar`s and a vector of `ExtendedPoints`. It is an +/// error to call this function with two vectors of different lengths. +/// +/// XXX this takes `edwards::ExtendedPoints` because we have to alloc scratch space here anyways, +/// and we need some space to store the converted points, so we may as well do the conversion here. +/// maybe there's a better way to avoid code duplication... however we can't quite just write a +/// generic `multiscalar_mult` because the non-vectorized code passes between models and this code +/// doesn't. +#[cfg(any(feature = "alloc", feature = "std"))] +pub fn multiscalar_mult<'a, 'b, I, J>(scalars: I, points: J) -> edwards::ExtendedPoint + where I: IntoIterator, + J: IntoIterator +{ + //assert_eq!(scalars.len(), points.len()); + + use clear_on_drop::ClearOnDrop; + let lookup_tables_vec: Vec<_> = points.into_iter() + .map(|P| LookupTable::from(ExtendedPoint::from(*P)) ) + .collect(); + + let lookup_tables = ClearOnDrop::new(lookup_tables_vec); + + // Setting s_i = i-th scalar, compute + // + // s_i = s_{i,0} + s_{i,1}*16^1 + ... + s_{i,63}*16^63, + // + // with `-8 ≤ s_{i,j} < 8` for `0 ≤ j < 63` and `-8 ≤ s_{i,63} ≤ 8`. + let scalar_digits_vec: Vec<_> = scalars.into_iter() + .map(|c| c.to_radix_16()) + .collect(); + + // The above puts the scalar digits into a heap-allocated Vec. + // To ensure that these are erased, pass ownership of the Vec into a + // ClearOnDrop wrapper. + let scalar_digits = ClearOnDrop::new(scalar_digits_vec); + + // Compute s_1*P_1 + ... + s_n*P_n: since + // + // s_i*P_i = P_i*(s_{i,0} + s_{i,1}*16^1 + ... + s_{i,63}*16^63) + // s_i*P_i = P_i*s_{i,0} + P_i*s_{i,1}*16^1 + ... + P_i*s_{i,63}*16^63 + // s_i*P_i = P_i*s_{i,0} + 16*(P_i*s_{i,1} + 16*( ... + 16*P_i*s_{i,63})...) + // + // we have the two-dimensional sum + // + // s_1*P_1 = P_1*s_{1,0} + 16*(P_1*s_{1,1} + 16*( ... + 16*P_1*s_{1,63})...) + // + s_2*P_2 = + P_2*s_{2,0} + 16*(P_2*s_{2,1} + 16*( ... + 16*P_2*s_{2,63})...) + // ... + // + s_n*P_n = + P_n*s_{n,0} + 16*(P_n*s_{n,1} + 16*( ... + 16*P_n*s_{n,63})...) + // + // We sum column-wise top-to-bottom, then right-to-left, + // multiplying by 16 only once per column. + // + // This provides the speedup over doing n independent scalar + // mults: we perform 63 multiplications by 16 instead of 63*n + // multiplications, saving 252*(n-1) doublings. + let mut Q = ExtendedPoint::identity(); + // XXX this algorithm makes no effort to be cache-aware; maybe it could be improved? + for j in (0..64).rev() { + Q = Q.mult_by_pow_2(4); + let it = scalar_digits.iter().zip(lookup_tables.iter()); + for (s_i, lookup_table_i) in it { + // Q = Q + s_{i,j} * P_i + Q = &Q + &lookup_table_i.select(s_i[j]); + } + } + Q.into() +} + +pub mod vartime { + //! Variable-time operations on curve points, useful for non-secret data. + use super::*; + + /// Holds odd multiples 1A, 3A, ..., 15A of a point A. + struct OddMultiples([ExtendedPoint; 8]); + + impl OddMultiples { + fn create(A: ExtendedPoint) -> OddMultiples { + // XXX would be great to skip this initialization + let mut Ai = [A; 8]; + let A2 = A.double(); + for i in 0..7 { + Ai[i+1] = &A2 + &Ai[i]; + } + // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A] + OddMultiples(Ai) + } + } + + impl Index for OddMultiples { + type Output = ExtendedPoint; + + fn index(&self, _index: usize) -> &ExtendedPoint { + &(self.0[_index]) + } + } + + /// Given a point `A` and scalars `a` and `b`, compute the point + /// `aA+bB`, where `B` is the Ed25519 basepoint (i.e., `B = (x,4/5)` + /// with x positive). + /// + /// This is the same as calling the iterator-based function, but slightly faster. + pub fn double_scalar_mult_basepoint(a: &Scalar, + A: &edwards::ExtendedPoint, + b: &Scalar) -> edwards::ExtendedPoint { + let a_naf = a.non_adjacent_form(); + let b_naf = b.non_adjacent_form(); + + // Find starting index + let mut i: usize = 255; + for j in (0..255).rev() { + i = j; + if a_naf[i] != 0 || b_naf[i] != 0 { + break; + } + } + + let odd_multiples_of_A = OddMultiples::create((*A).into()); + let odd_multiples_of_B = &avx2::constants::ODD_MULTIPLES_OF_BASEPOINT; + + let mut Q = ExtendedPoint::identity(); + + loop { + Q = Q.double(); + + if a_naf[i] > 0 { + Q = &Q + &odd_multiples_of_A[( a_naf[i]/2) as usize]; + } else if a_naf[i] < 0 { + Q = &Q - &odd_multiples_of_A[(-a_naf[i]/2) as usize]; + } + + if b_naf[i] > 0 { + Q = &Q + &odd_multiples_of_B[( b_naf[i]/2) as usize]; + } else if b_naf[i] < 0 { + Q = &Q - &odd_multiples_of_B[(-b_naf[i]/2) as usize]; + } + + if i == 0 { + break; + } + i -= 1; + } + + Q.into() + } + + /// Given a vector of public scalars and a vector of public points, compute + /// $$ + /// Q = c\_1 P\_1 + \cdots + c\_n P\_n. + /// $$ + /// + /// # Input + /// + /// A vector of `Scalar`s and a vector of `ExtendedPoints`. It is an + /// error to call this function with two vectors of different lengths. + #[cfg(any(feature = "alloc", feature = "std"))] + pub fn multiscalar_mult<'a, 'b, I, J>(scalars: I, points: J) -> edwards::ExtendedPoint + where I: IntoIterator, + J: IntoIterator + { + //assert_eq!(scalars.len(), points.len()); + + let nafs: Vec<_> = scalars.into_iter() + .map(|c| c.non_adjacent_form()).collect(); + + let odd_multiples: Vec<_> = points.into_iter() + .map(|P| OddMultiples::create((*P).into()) ).collect(); + + let mut Q = ExtendedPoint::identity(); + + for i in (0..255).rev() { + Q = Q.double(); + + for (naf, odd_multiple) in nafs.iter().zip(odd_multiples.iter()) { + if naf[i] > 0 { + Q = &Q + &odd_multiple[( naf[i]/2) as usize]; + } else if naf[i] < 0 { + Q = &Q - &odd_multiple[(-naf[i]/2) as usize]; + } + } + } + Q.into() + } +} + +#[cfg(test)] +mod test { + use super::*; + + use constants; + + fn serial_add(P: edwards::ExtendedPoint, Q: edwards::ExtendedPoint) -> edwards::ExtendedPoint { + use backend::u64::field::FieldElement64; + + let (X1, Y1, Z1, T1) = (P.X, P.Y, P.Z, P.T); + let (X2, Y2, Z2, T2) = (Q.X, Q.Y, Q.Z, Q.T); + + macro_rules! print_var { + ($x:ident) => { + println!("{} = {:?}", stringify!($x), $x.to_bytes()); + } + } + + let S0 = &Y1 - &X1; // R1 + let S1 = &Y1 + &X1; // R3 + let S2 = &Y2 - &X2; // R2 + let S3 = &Y2 + &X2; // R4 + print_var!(S0); + print_var!(S1); + print_var!(S2); + print_var!(S3); + println!(""); + + let S4 = &S0 * &S2; // R5 = R1 * R2 + let S5 = &S1 * &S3; // R6 = R3 * R4 + let S6 = &Z1 * &Z2; // R8 + let S7 = &T1 * &T2; // R7 + print_var!(S4); + print_var!(S5); + print_var!(S6); + print_var!(S7); + println!(""); + + let S8 = &S4 * &FieldElement64([ 121666,0,0,0,0]); // R5 + let S9 = &S5 * &FieldElement64([ 121666,0,0,0,0]); // R6 + let S10 = &S6 * &FieldElement64([2*121666,0,0,0,0]); // R8 + let S11 = &S7 * &(-&FieldElement64([2*121665,0,0,0,0])); // R7 + print_var!(S8 ); + print_var!(S9 ); + print_var!(S10); + print_var!(S11); + println!(""); + + let S12 = &S9 - &S8; // R1 + let S13 = &S9 + &S8; // R4 + let S14 = &S10 - &S11; // R2 + let S15 = &S10 + &S11; // R3 + print_var!(S12); + print_var!(S13); + print_var!(S14); + print_var!(S15); + println!(""); + + let X3 = &S12 * &S14; // R1 * R2 + let Y3 = &S15 * &S13; // R3 * R4 + let Z3 = &S15 * &S14; // R2 * R3 + let T3 = &S12 * &S13; // R1 * R4 + + edwards::ExtendedPoint{X: X3, Y: Y3, Z: Z3, T: T3} + } + + fn addition_test_helper(P: edwards::ExtendedPoint, Q: edwards::ExtendedPoint) { + // Test the serial implementation of the parallel addition formulas + let R_serial: edwards::ExtendedPoint = serial_add(P.into(), Q.into()).into(); + // Test the vector implementation of the parallel addition formulas + let R_vector: edwards::ExtendedPoint = (&ExtendedPoint::from(P) + &ExtendedPoint::from(Q)).into(); + // Test the vector implementation of the parallel subtraction formulas + let S_vector: edwards::ExtendedPoint = (&ExtendedPoint::from(P) - &ExtendedPoint::from(Q)).into(); + + // Test the vector implementation of the parallel readdition formulas + let cached_Q = CachedPoint::from(ExtendedPoint::from(Q)); + let T_vector: edwards::ExtendedPoint = (&ExtendedPoint::from(P) + &cached_Q).into(); + + println!("Testing point addition:"); + println!("P = {:?}", P); + println!("Q = {:?}", Q); + println!("cached Q = {:?}", cached_Q); + println!("R = P + Q = {:?}", &P + &Q); + println!("R_serial = {:?}", R_serial); + println!("R_vector = {:?}", R_vector); + println!("T_vector = {:?}", T_vector); + println!("S = P - Q = {:?}", &P - &Q); + println!("S_vector = {:?}", S_vector); + assert_eq!(R_serial.compress(), (&P + &Q).compress()); + assert_eq!(R_vector.compress(), (&P + &Q).compress()); + assert_eq!(T_vector.compress(), (&P + &Q).compress()); + assert_eq!(S_vector.compress(), (&P - &Q).compress()); + println!("OK!\n"); + } + + #[test] + fn sub_vs_add_minus() { + let P: ExtendedPoint = edwards::ExtendedPoint::identity().into(); + let Q: ExtendedPoint = edwards::ExtendedPoint::identity().into(); + + let mQ = -&Q; + + println!("sub"); + let R1: edwards::ExtendedPoint = (&P - &Q).into(); + println!("add neg"); + let R2: edwards::ExtendedPoint = (&P + &mQ).into(); + + assert_eq!(R2.compress(), edwards::ExtendedPoint::identity().compress()); + assert_eq!(R1.compress(), edwards::ExtendedPoint::identity().compress()); + } + + + #[test] + fn vector_addition_vs_serial_addition_vs_edwards_extendedpoint() { + use constants; + use scalar::Scalar; + + println!("Testing id +- id"); + let P = edwards::ExtendedPoint::identity(); + let Q = edwards::ExtendedPoint::identity(); + addition_test_helper(P, Q); + + println!("Testing id +- B"); + let P = edwards::ExtendedPoint::identity(); + let Q = constants::ED25519_BASEPOINT_POINT; + addition_test_helper(P, Q); + + println!("Testing B +- B"); + let P = constants::ED25519_BASEPOINT_POINT; + let Q = constants::ED25519_BASEPOINT_POINT; + addition_test_helper(P, Q); + + println!("Testing B +- kB"); + let P = constants::ED25519_BASEPOINT_POINT; + let Q = &constants::ED25519_BASEPOINT_TABLE * &Scalar::from_u64(8475983829); + addition_test_helper(P, Q); + } + + fn serial_double(P: edwards::ExtendedPoint) -> edwards::ExtendedPoint { + let (X1, Y1, Z1, T1) = (P.X, P.Y, P.Z, P.T); + + macro_rules! print_var { + ($x:ident) => { + println!("{} = {:?}", stringify!($x), $x.to_bytes()); + } + } + + let S0 = &X1 + &Y1; // R1 + print_var!(S0); + println!(""); + + let S1 = X1.square(); + let S2 = Y1.square(); + let S3 = Z1.square(); + let S4 = S0.square(); + print_var!(S1); + print_var!(S2); + print_var!(S3); + print_var!(S4); + println!(""); + + let S5 = &S1 + &S2; + let S6 = &S1 - &S2; + let S7 = &S3 + &S3; + let S8 = &S7 + &S6; + let S9 = &S5 - &S4; + print_var!(S5); + print_var!(S6); + print_var!(S7); + print_var!(S8); + print_var!(S9); + println!(""); + + let X3 = &S8 * &S9; + let Y3 = &S5 * &S6; + let Z3 = &S8 * &S6; + let T3 = &S5 * &S9; + + edwards::ExtendedPoint{X: X3, Y: Y3, Z: Z3, T: T3} + } + + fn doubling_test_helper(P: edwards::ExtendedPoint) { + let R1: edwards::ExtendedPoint = serial_double(P.into()).into(); + let R2: edwards::ExtendedPoint = ExtendedPoint::from(P).double().into(); + println!("Testing point doubling:"); + println!("P = {:?}", P); + println!("(serial) R1 = {:?}", R1); + println!("(vector) R2 = {:?}", R2); + println!("P + P = {:?}", &P + &P); + assert_eq!(R1.compress(), (&P + &P).compress()); + assert_eq!(R2.compress(), (&P + &P).compress()); + println!("OK!\n"); + } + + #[test] + fn vector_doubling_vs_serial_doubling_vs_edwards_extendedpoint() { + use constants; + use scalar::Scalar; + + println!("Testing [2]id"); + let P = edwards::ExtendedPoint::identity(); + doubling_test_helper(P); + + println!("Testing [2]B"); + let P = constants::ED25519_BASEPOINT_POINT; + doubling_test_helper(P); + + println!("Testing [2]([k]B)"); + let P = &constants::ED25519_BASEPOINT_TABLE * &Scalar::from_u64(8475983829); + doubling_test_helper(P); + } + + #[test] + fn identity_trait_vs_edwards_identity() { + let id1: edwards::ExtendedPoint = ExtendedPoint::identity().into(); + let id2: edwards::ExtendedPoint = edwards::ExtendedPoint::identity(); + assert_eq!(id1.compress(), id2.compress()); + } + + #[test] + fn neg_vs_edwards_neg() { + let B: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.into(); + let Bneg = -&B; + assert_eq!(edwards::ExtendedPoint::from(Bneg).compress(), + (-&constants::ED25519_BASEPOINT_POINT).compress()); + } + + #[test] + fn scalar_mult_vs_edwards_scalar_mult() { + let B: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.into(); + // some random bytes + let s = Scalar::from_bits([233, 1, 233, 147, 113, 78, 244, 120, 40, 45, 103, 51, 224, 199, 189, 218, 96, 140, 211, 112, 39, 194, 73, 216, 173, 33, 102, 93, 76, 200, 84, 12]); + + let R1 = edwards::ExtendedPoint::from(&B * &s); + let R2 = &constants::ED25519_BASEPOINT_TABLE * &s; + + assert_eq!(R1.compress(), R2.compress()); + } + + #[test] + fn scalar_mult_vs_basepoint_table_scalar_mult() { + let B: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.into(); + let B_table = EdwardsBasepointTable::create(&B); + // some random bytes + let s = Scalar::from_bits([233, 1, 233, 147, 113, 78, 244, 120, 40, 45, 103, 51, 224, 199, 189, 218, 96, 140, 211, 112, 39, 194, 73, 216, 173, 33, 102, 93, 76, 200, 84, 12]); + + let P1 = &B * &s; + let P2 = &B_table * &s; + + assert_eq!(edwards::ExtendedPoint::from(P1).compress(), + edwards::ExtendedPoint::from(P2).compress()); + } + + #[test] + fn multiscalar_mult_vs_adding_scalar_mults() { + let B: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.into(); + let s1 = Scalar::from_bits([233, 1, 233, 147, 113, 78, 244, 120, 40, 45, 103, 51, 224, 199, 189, 218, 96, 140, 211, 112, 39, 194, 73, 216, 173, 33, 102, 93, 76, 200, 84, 12]); + let s2 = Scalar::from_bits([165, 30, 79, 89, 58, 24, 195, 245, 248, 146, 203, 236, 119, 43, 64, 119, 196, 111, 188, 251, 248, 53, 234, 59, 215, 28, 218, 13, 59, 120, 14, 4]); + + let P1 = &B * &s2; + let P2 = &B * &s1; + + let R = &(&P1 * &s1) + &(&P2 * &s2); + + let R_multiscalar = multiscalar_mult(&[s1, s2], &[P1.into(), P2.into()]); + + assert_eq!(edwards::ExtendedPoint::from(R).compress(), + R_multiscalar.compress()); + } + + mod vartime { + use super::*; + + #[test] + fn multiscalar_mult_vs_adding_scalar_mults() { + let B: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.into(); + let s1 = Scalar::from_bits([233, 1, 233, 147, 113, 78, 244, 120, 40, 45, 103, 51, 224, 199, 189, 218, 96, 140, 211, 112, 39, 194, 73, 216, 173, 33, 102, 93, 76, 200, 84, 12]); + let s2 = Scalar::from_bits([165, 30, 79, 89, 58, 24, 195, 245, 248, 146, 203, 236, 119, 43, 64, 119, 196, 111, 188, 251, 248, 53, 234, 59, 215, 28, 218, 13, 59, 120, 14, 4]); + + let P1 = &B * &s2; + let P2 = &B * &s1; + + let R = &(&P1 * &s1) + &(&P2 * &s2); + + let R_multiscalar = vartime::multiscalar_mult(&[s1, s2], &[P1.into(), P2.into()]); + + assert_eq!(edwards::ExtendedPoint::from(R).compress(), + R_multiscalar.compress()); + } + } +} + +#[cfg(all(test, feature = "bench"))] +mod bench { + use test::Bencher; + use rand::OsRng; + use super::*; + + use constants; + use scalar::Scalar; + + #[bench] + fn conversion_into__avx2_format(b: &mut Bencher) { + let B = constants::ED25519_BASEPOINT_POINT; + + b.iter(|| ExtendedPoint::from(B)); + } + + #[bench] + fn conversion_outof_avx2_format(b: &mut Bencher) { + let B = constants::ED25519_BASEPOINT_POINT; + let B_avx2 = ExtendedPoint::from(B); + + b.iter(|| edwards::ExtendedPoint::from(B_avx2)); + } + + #[bench] + fn point_readdition(b: &mut Bencher) { + let B = &constants::ED25519_BASEPOINT_TABLE; + let P = ExtendedPoint::from(B * &Scalar::from_u64(83973422)); + let Q = ExtendedPoint::from(B * &Scalar::from_u64(98932328)); + let Q_cached = CachedPoint::from(Q); + + b.iter(|| &P + &Q_cached ); + } + + #[bench] + fn point_addition(b: &mut Bencher) { + let B = &constants::ED25519_BASEPOINT_TABLE; + let P = ExtendedPoint::from(B * &Scalar::from_u64(83973422)); + let Q = ExtendedPoint::from(B * &Scalar::from_u64(98932328)); + + b.iter(|| &P + &Q ); + } + + #[bench] + fn point_doubling(b: &mut Bencher) { + let B = &constants::ED25519_BASEPOINT_TABLE; + let P = ExtendedPoint::from(B * &Scalar::from_u64(83973422)); + + b.iter(|| P.double() ); + } + + #[bench] + fn scalar_mult(b: &mut Bencher) { + let B = &constants::ED25519_BASEPOINT_TABLE; + let P = ExtendedPoint::from(B * &Scalar::from_u64(83973422)); + let s = Scalar::from_bits([233, 1, 233, 147, 113, 78, 244, 120, 40, 45, 103, 51, 224, 199, 189, 218, 96, 140, 211, 112, 39, 194, 73, 216, 173, 33, 102, 93, 76, 200, 84, 12]); + + b.iter(|| &P * &s ); + } + + #[bench] + fn basepoint_table_creation(b: &mut Bencher) { + let B = ExtendedPoint::from(constants::ED25519_BASEPOINT_POINT); + + b.iter(|| EdwardsBasepointTable::create(&B) ); + } + + #[bench] + fn basepoint_mult(b: &mut Bencher) { + let B = ExtendedPoint::from(constants::ED25519_BASEPOINT_POINT); + let table = EdwardsBasepointTable::create(&B); + let s = Scalar::from_bits([233, 1, 233, 147, 113, 78, 244, 120, 40, 45, 103, 51, 224, 199, 189, 218, 96, 140, 211, 112, 39, 194, 73, 216, 173, 33, 102, 93, 76, 200, 84, 12]); + + b.iter(|| &table * &s ); + } + + #[bench] + fn ten_fold_scalar_mult(b: &mut Bencher) { + let mut csprng: OsRng = OsRng::new().unwrap(); + // Create 10 random scalars + let scalars: Vec<_> = (0..10).map(|_| Scalar::random(&mut csprng)).collect(); + // Create 10 points (by doing scalar mults) + let B = &constants::ED25519_BASEPOINT_POINT; + let points: Vec<_> = scalars.iter().map(|s| B * &s).collect(); + + b.iter(|| multiscalar_mult(&scalars, &points)); + } + + mod vartime { + use super::super::*; + use super::{constants, Bencher, OsRng}; + + #[bench] + fn double_scalar_mult(b: &mut Bencher) { + let mut csprng: OsRng = OsRng::new().unwrap(); + // Create 2 random scalars + let s1 = Scalar::random(&mut csprng); + let s2 = Scalar::random(&mut csprng); + let B = constants::ED25519_BASEPOINT_POINT; + let P = &B * &s1; + + b.iter(|| vartime::double_scalar_mult_basepoint(&s2, &P, &s1) ); + } + + #[bench] + fn ten_fold_scalar_mult(b: &mut Bencher) { + let mut csprng: OsRng = OsRng::new().unwrap(); + // Create 10 random scalars + let scalars: Vec<_> = (0..10).map(|_| Scalar::random(&mut csprng)).collect(); + // Create 10 points (by doing scalar mults) + let B = &constants::ED25519_BASEPOINT_POINT; + let points: Vec<_> = scalars.iter().map(|s| B * &s).collect(); + + b.iter(|| vartime::multiscalar_mult(&scalars, &points)); + } + } +} + diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs new file mode 100644 index 0000000..7f422e0 --- /dev/null +++ b/src/backend/avx2/field.rs @@ -0,0 +1,689 @@ +// -*- mode: rust; coding: utf-8; -*- +// +// This file is part of curve25519-dalek. +// Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence +// See LICENSE for licensing information. +// +// Authors: +// - Isis Agora Lovecruft +// - Henry de Valence + +//! 4-way vectorized 32bit field arithmetic using AVX2. +//! + +#![allow(bad_style)] + +pub const A_LANES: u8 = 0b0000_0101; +pub const B_LANES: u8 = 0b0000_1010; +pub const C_LANES: u8 = 0b0101_0000; +pub const D_LANES: u8 = 0b1010_0000; + +pub const A_LANES64: u8 = 0b00_00_00_11; +pub const B_LANES64: u8 = 0b00_00_11_00; +pub const C_LANES64: u8 = 0b00_11_00_00; +pub const D_LANES64: u8 = 0b11_00_00_00; + +pub const ALL_LANES: u8 = A_LANES | B_LANES | C_LANES | D_LANES; + +use std::ops::Mul; + +use stdsimd::simd::{u32x8, i32x8, u64x4}; + +use backend::u64::field::FieldElement64; + +use backend::avx2::constants::{P_TIMES_2_LO, P_TIMES_2_HI, P_TIMES_16_LO, P_TIMES_16_HI}; + +/// A vector of four `FieldElements`, implemented using AVX2. +#[derive(Clone, Copy, Debug)] +pub(crate) struct FieldElement32x4(pub(crate) [u32x8; 5]); + +use subtle::ConditionallyAssignable; + +impl ConditionallyAssignable for FieldElement32x4 { + fn conditional_assign(&mut self, other: &FieldElement32x4, choice: u8) { + let mask = (-(choice as i32)) as u32; + let mask_vec = u32x8::splat(mask); + for i in 0..5 { + self.0[i] = self.0[i] ^ (mask_vec & (self.0[i] ^ other.0[i])); + } + } +} + + +impl FieldElement32x4 { + pub(crate) fn split(&self) -> [FieldElement64; 4] { + let mut out = [FieldElement64::zero(); 4]; + for i in 0..5 { + + let a_2i = self.0[i].extract(0) as u64; // + let b_2i = self.0[i].extract(1) as u64; // + let a_2i_1 = self.0[i].extract(2) as u64; // `. + let b_2i_1 = self.0[i].extract(3) as u64; // | pre-swapped to avoid + let c_2i = self.0[i].extract(4) as u64; // | a cross lane shuffle + let d_2i = self.0[i].extract(5) as u64; // .' + let c_2i_1 = self.0[i].extract(6) as u64; // + let d_2i_1 = self.0[i].extract(7) as u64; // + + out[0].0[i] = a_2i + (a_2i_1 << 26); + out[1].0[i] = b_2i + (b_2i_1 << 26); + out[2].0[i] = c_2i + (c_2i_1 << 26); + out[3].0[i] = d_2i + (d_2i_1 << 26); + } + + out + } + + pub fn zero() -> FieldElement32x4 { + FieldElement32x4([u32x8::splat(0);5]) + } + + pub fn splat(x: &FieldElement64) -> FieldElement32x4 { + FieldElement32x4::new(x,x,x,x) + } + + pub fn new( + x0: &FieldElement64, + x1: &FieldElement64, + x2: &FieldElement64, + x3: &FieldElement64, + ) -> FieldElement32x4 { + let mut buf = [u32x8::splat(0); 5]; + let low_26_bits = (1 << 26) - 1; + for i in 0..5 { + let a_2i = (x0.0[i] & low_26_bits) as u32; + let a_2i_1 = (x0.0[i] >> 26) as u32; + let b_2i = (x1.0[i] & low_26_bits) as u32; + let b_2i_1 = (x1.0[i] >> 26) as u32; + let c_2i = (x2.0[i] & low_26_bits) as u32; + let c_2i_1 = (x2.0[i] >> 26) as u32; + let d_2i = (x3.0[i] & low_26_bits) as u32; + let d_2i_1 = (x3.0[i] >> 26) as u32; + + buf[i] = u32x8::new(a_2i, b_2i, a_2i_1, b_2i_1, c_2i, d_2i, c_2i_1, d_2i_1); + } + + let mut out = FieldElement32x4(buf); + out.reduce32(); + return out; + } + + pub fn negate_lazy(&mut self, mask: u8) { + let mask = mask as i32; + unsafe { + use stdsimd::vendor::_mm256_blend_epi32; + self.0[0] = _mm256_blend_epi32(self.0[0].into(), (P_TIMES_2_LO - self.0[0]).into(), mask).into(); + self.0[1] = _mm256_blend_epi32(self.0[1].into(), (P_TIMES_2_HI - self.0[1]).into(), mask).into(); + self.0[2] = _mm256_blend_epi32(self.0[2].into(), (P_TIMES_2_HI - self.0[2]).into(), mask).into(); + self.0[3] = _mm256_blend_epi32(self.0[3].into(), (P_TIMES_2_HI - self.0[3]).into(), mask).into(); + self.0[4] = _mm256_blend_epi32(self.0[4].into(), (P_TIMES_2_HI - self.0[4]).into(), mask).into(); + } + } + + /// Negate variables in lanes where mask is set + pub fn negate(&mut self, mask: u8) { + let mask = mask as i32; + unsafe { + use stdsimd::vendor::_mm256_blend_epi32; + self.0[0] = _mm256_blend_epi32(self.0[0].into(), (P_TIMES_16_LO - self.0[0]).into(), mask).into(); + self.0[1] = _mm256_blend_epi32(self.0[1].into(), (P_TIMES_16_HI - self.0[1]).into(), mask).into(); + self.0[2] = _mm256_blend_epi32(self.0[2].into(), (P_TIMES_16_HI - self.0[2]).into(), mask).into(); + self.0[3] = _mm256_blend_epi32(self.0[3].into(), (P_TIMES_16_HI - self.0[3]).into(), mask).into(); + self.0[4] = _mm256_blend_epi32(self.0[4].into(), (P_TIMES_16_HI - self.0[4]).into(), mask).into(); + } + self.reduce32(); + } + + /// Given `self = (A,B,C,D)`, set `self = (B,A,C,D)` + pub fn swap_AB(&mut self) { + unsafe { + use stdsimd::vendor::_mm256_shuffle_epi32; + use stdsimd::vendor::_mm256_blend_epi32; + for i in 0..5 { + let swapped = _mm256_shuffle_epi32(self.0[i].into(), 0b10_11_00_01); + self.0[i] = _mm256_blend_epi32(self.0[i].into(), swapped, 0b00001111).into(); + } + } + } + + /// Given `self = (A,B,C,D)`, set `self = (A,B,D,C)` + pub fn swap_CD(&mut self) { + unsafe { + use stdsimd::vendor::_mm256_shuffle_epi32; + use stdsimd::vendor::_mm256_blend_epi32; + for i in 0..5 { + let swapped = _mm256_shuffle_epi32(self.0[i].into(), 0b10_11_00_01); + self.0[i] = _mm256_blend_epi32(self.0[i].into(), swapped, 0b11110000).into(); + } + } + } + + /// Given `self = (A,B,C,D)`, set `self = (B - A, B + A, D - C, D + C)` according to `mask`. + pub fn diff_sum(&mut self, mask: u8) { + let mask = mask as i32; + unsafe { + use stdsimd::vendor::{_mm256_shuffle_epi32, _mm256_blend_epi32}; + + let x01 = self.0[0]; + let x01_shuf = _mm256_shuffle_epi32(x01.as_i32x8(), 0b10_11_00_01).as_u32x8(); + let v1 = (x01_shuf + P_TIMES_2_LO) - x01; + let v2 = x01_shuf + x01; + let diffsum01 = _mm256_blend_epi32(v1.into(), v2.into(), 0b10101010).as_u32x8(); + self.0[0] = _mm256_blend_epi32(x01.into(), diffsum01.into(), mask).into(); + + let x23 = self.0[1]; + let x23_shuf = _mm256_shuffle_epi32(x23.as_i32x8(), 0b10_11_00_01).as_u32x8(); + let v1 = (x23_shuf + P_TIMES_2_HI) - x23; + let v2 = x23_shuf + x23; + let diffsum23 = _mm256_blend_epi32(v1.into(), v2.into(), 0b10101010).as_u32x8(); + self.0[1] = _mm256_blend_epi32(x23.into(), diffsum23.into(), mask).into(); + + let x45 = self.0[2]; + let x45_shuf = _mm256_shuffle_epi32(x45.as_i32x8(), 0b10_11_00_01).as_u32x8(); + let v1 = (x45_shuf + P_TIMES_2_HI) - x45; + let v2 = x45_shuf + x45; + let diffsum45 = _mm256_blend_epi32(v1.into(), v2.into(), 0b10101010).as_u32x8(); + self.0[2] = _mm256_blend_epi32(x45.into(), diffsum45.into(), mask).into(); + + let x67 = self.0[3]; + let x67_shuf = _mm256_shuffle_epi32(x67.as_i32x8(), 0b10_11_00_01).as_u32x8(); + let v1 = (x67_shuf + P_TIMES_2_HI) - x67; + let v2 = x67_shuf + x67; + let diffsum67 = _mm256_blend_epi32(v1.into(), v2.into(), 0b10101010).as_u32x8(); + self.0[3] = _mm256_blend_epi32(x67.into(), diffsum67.into(), mask).into(); + + let x89 = self.0[4]; + let x89_shuf = _mm256_shuffle_epi32(x89.as_i32x8(), 0b10_11_00_01).as_u32x8(); + let v1 = (x89_shuf + P_TIMES_2_HI) - x89; + let v2 = x89_shuf + x89; + let diffsum89 = _mm256_blend_epi32(v1.into(), v2.into(), 0b10101010).as_u32x8(); + self.0[4] = _mm256_blend_epi32(x89.into(), diffsum89.into(), mask).into(); + } + } + + /// Let `self` \\(= (A, B, C, D) \\). + /// + /// Compute + /// $$( 121666A, 121666B, 2\cdot 121666C, 2\cdot 121665 D).$$ + pub fn scale_by_curve_constants(&mut self) { + let mut b = [u64x4::splat(0); 10]; + + let consts = u32x8::new(121666, 0, 121666, 0, 2*121666, 0, 2*121665, 0); + + unsafe { + use stdsimd::vendor::_mm256_mul_epu32; + + let (b0, b1) = unpack_pair(self.0[0]); + b[0] = _mm256_mul_epu32(b0, consts); + b[1] = _mm256_mul_epu32(b1, consts); + + let (b2, b3) = unpack_pair(self.0[1]); + b[2] = _mm256_mul_epu32(b2, consts); + b[3] = _mm256_mul_epu32(b3, consts); + + let (b4, b5) = unpack_pair(self.0[2]); + b[4] = _mm256_mul_epu32(b4, consts); + b[5] = _mm256_mul_epu32(b5, consts); + + let (b6, b7) = unpack_pair(self.0[3]); + b[6] = _mm256_mul_epu32(b6, consts); + b[7] = _mm256_mul_epu32(b7, consts); + + let (b8, b9) = unpack_pair(self.0[4]); + b[8] = _mm256_mul_epu32(b8, consts); + b[9] = _mm256_mul_epu32(b9, consts); + } + + *self = FieldElement32x4::reduce64(b); + } + + pub fn reduce32(&mut self) { + + let shifts = i32x8::new(26,26,25,25,26,26,25,25); + let masks = u32x8::new((1<<26)-1, (1<<26)-1, (1<<25)-1, (1<<25)-1, + (1<<26)-1, (1<<26)-1, (1<<25)-1, (1<<25)-1); + + let carry = |v: u32x8| -> u32x8 { + unsafe { + use stdsimd::vendor::_mm256_srlv_epi32; + _mm256_srlv_epi32(v.into(), shifts).into() + } + }; + + let swap_lanes = |v: u32x8| -> u32x8 { + unsafe { + use stdsimd::vendor::_mm256_shuffle_epi32; + _mm256_shuffle_epi32(v.into(), 0b01_00_11_10).into() + } + }; + + let combine = |v_lo: u32x8, v_hi: u32x8| -> u32x8 { + unsafe { + use stdsimd::vendor::_mm256_blend_epi32; + _mm256_blend_epi32(v_lo.into(), v_hi.into(), 0b11_00_11_00).into() + } + }; + + let v = &mut self.0; + + let c10 = swap_lanes(carry(v[0])); + v[0] = (v[0] & masks) + combine(u32x8::splat(0), c10); + let c32 = swap_lanes(carry(v[1])); + v[1] = (v[1] & masks) + combine(c10, c32); + let c54 = swap_lanes(carry(v[2])); + v[2] = (v[2] & masks) + combine(c32, c54); + let c76 = swap_lanes(carry(v[3])); + v[3] = (v[3] & masks) + combine(c54, c76); + let c98 = swap_lanes(carry(v[4])); + v[4] = (v[4] & masks) + combine(c76, c98); + + // Still need to account for c9 + // c98 = (c9, c9, c8, c8, c9, c9, c8, c8) + // + let c9_19: u32x8; + unsafe { + use stdsimd::vendor::_mm256_mul_epu32; + use stdsimd::vendor::_mm256_shuffle_epi32; + let c9_spread: u32x8 = _mm256_shuffle_epi32(c98.into(), 0b11_01_10_00).into(); + let c9_19_spread: u32x8 = _mm256_mul_epu32(c9_spread, u64x4::splat(19).into()).into(); + c9_19 = _mm256_shuffle_epi32(c9_19_spread.into(), 0b11_01_10_00).into(); + } + + v[0] = v[0] + c9_19; + } + + pub fn reduce64(mut z: [u64x4; 10]) -> FieldElement32x4 { + // These aren't const because splat isn't a const fn + let LOW_25_BITS: u64x4 = u64x4::splat((1<<25)-1); + let LOW_26_BITS: u64x4 = u64x4::splat((1<<26)-1); + + // Carry the value from limb i = 0..8 to limb i+1 + let carry = |z: &mut [u64x4; 10], i: usize| { + debug_assert!(i < 9); + if i % 2 == 0 { + // Even limbs have 26 bits + z[i+1] = z[i+1] + (z[i] >> 26); + z[i] = z[i] & LOW_26_BITS; + } else { + // Odd limbs have 25 bits + z[i+1] = z[i+1] + (z[i] >> 25); + z[i] = z[i] & LOW_25_BITS; + } + }; + + // 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) + + // Last carry has a multiplication by 19. In the serial case we + // do a 64-bit multiplication by 19, but here we want to do a + // 32-bit multiplication. However, if we only know z[9] < 2^64, + // the carry is bounded as c < 2^(64-25) = 2^39, which is too + // big. To ensure c < 2^32, we would need z[9] < 2^57. + // Instead, we split the carry in two, with c = c_0 + c_1*2^26. + + let c = z[9] >> 25; + z[9] = z[9] & LOW_25_BITS; + let mut c0 = c & LOW_26_BITS; // c0 < 2^26; + let mut c1 = c >> 26; // c1 < 2^(39-26) = 2^13; + + unsafe { + use stdsimd::vendor::_mm256_mul_epu32; + let x19 = u32x8::from(u64x4::splat(19)); + c0 = _mm256_mul_epu32(u32x8::from(c0), x19); // c0 < 2^30.25 + c1 = _mm256_mul_epu32(u32x8::from(c1), x19); // c1 < 2^17.25 + } + + z[0] = z[0] + c0; // z0 < 2^26 + 2^30.25 < 2^30.33 + z[1] = z[1] + c1; // z1 < 2^25 + 2^17.25 < 2^25.0067 + carry(&mut z, 0); // z0 < 2^26, z1 < 2^25.0067 + 2^4.33 = 2^25.007 + + // Now repack the [u64x4; 10] into a FieldElement32x4 + + FieldElement32x4([ + repack_pair(z[0].into(), z[1].into()), + repack_pair(z[2].into(), z[3].into()), + repack_pair(z[4].into(), z[5].into()), + repack_pair(z[6].into(), z[7].into()), + repack_pair(z[8].into(), z[9].into()), + ]) + } +} + +#[inline(always)] +pub fn unpack_pair(src: u32x8) -> (u32x8, u32x8) { + let a: u32x8; + let b: u32x8; + let zero = i32x8::new(0,0,0,0,0,0,0,0); + unsafe { + use stdsimd::vendor::_mm256_unpackhi_epi32; + use stdsimd::vendor::_mm256_unpacklo_epi32; + a = _mm256_unpacklo_epi32(src.as_i32x8(), zero).as_u32x8(); + b = _mm256_unpackhi_epi32(src.as_i32x8(), zero).as_u32x8(); + } + (a,b) +} + +#[inline(always)] +pub fn repack_pair(x: u32x8, y: u32x8) -> u32x8 { + unsafe { + use stdsimd::vendor::_mm256_shuffle_epi32; + use stdsimd::vendor::_mm256_blend_epi32; + + // Input: x = (a0, 0, b0, 0, c0, 0, d0) + // Input: y = (a1, 0, b1, 0, c1, 0, d1) + + let x_shuffled = _mm256_shuffle_epi32(x.into(), 0b11_01_10_00); + let y_shuffled = _mm256_shuffle_epi32(y.into(), 0b10_00_11_01); + + // x' = (a0, b0, 0, 0, c0, d0, 0, 0) + // y' = ( 0, 0, a1, b1, 0, 0, c1, d1) + + return _mm256_blend_epi32(x_shuffled, y_shuffled, 0b11001100).as_u32x8(); + } +} + +impl FieldElement32x4 { + /// Square this field element, then conditionally negate according to `neg_mask`; for instance, + /// `neg_mask = 0b11_00_00_00` negates the \\( D \\) value. + /// + /// # Precondition + /// + /// Limbs must be bounded by bit-excess \\( b < 2.0 \\). + pub fn square(&self, neg_mask: u8) -> FieldElement32x4 { + #[inline(always)] + fn m(x: u32x8, y: u32x8) -> u64x4 { + use stdsimd::vendor::_mm256_mul_epu32; + unsafe { _mm256_mul_epu32(x,y) } + } + + #[inline(always)] + fn m_lo(x: u32x8, y: u32x8) -> u32x8 { + use stdsimd::vendor::_mm256_mul_epu32; + unsafe { u32x8::from(_mm256_mul_epu32(x,y)) } + } + + let v19 = u32x8::new(19,0,19,0,19,0,19,0); + + let (x0, x1) = unpack_pair(self.0[0]); + let (x2, x3) = unpack_pair(self.0[1]); + let (x4, x5) = unpack_pair(self.0[2]); + let (x6, x7) = unpack_pair(self.0[3]); + let (x8, x9) = unpack_pair(self.0[4]); + + let x0_2 = x0 << 1; + let x1_2 = x1 << 1; + let x2_2 = x2 << 1; + let x3_2 = x3 << 1; + let x4_2 = x4 << 1; + let x5_2 = x5 << 1; + let x6_2 = x6 << 1; + let x7_2 = x7 << 1; + + let x5_19 = m_lo(v19, x5); + let x6_19 = m_lo(v19, x6); + let x7_19 = m_lo(v19, x7); + let x8_19 = m_lo(v19, x8); + let x9_19 = m_lo(v19, x9); + + let mut z0 = m(x0, x0) + m(x2_2,x8_19) + m(x4_2,x6_19) + ((m(x1_2,x9_19) + m(x3_2,x7_19) + m(x5,x5_19)) << 1); + let mut z1 = m(x0_2,x1) + m(x3_2,x8_19) + m(x5_2,x6_19) + ((m(x2,x9_19) + m(x4,x7_19)) << 1); + let mut z2 = m(x0_2,x2) + m(x1_2,x1) + m(x4_2,x8_19) + m(x6,x6_19) + ((m(x3_2,x9_19) + m(x5_2,x7_19)) << 1); + let mut z3 = m(x0_2,x3) + m(x1_2,x2) + m(x5_2,x8_19) + ((m(x4,x9_19) + m(x6,x7_19)) << 1); + let mut z4 = m(x0_2,x4) + m(x1_2,x3_2) + m(x2, x2) + m(x6_2,x8_19) + ((m(x5_2,x9_19) + m(x7,x7_19)) << 1); + let mut z5 = m(x0_2,x5) + m(x1_2,x4) + m(x2_2,x3) + m(x7_2,x8_19) + ((m(x6,x9_19)) << 1); + let mut z6 = m(x0_2,x6) + m(x1_2,x5_2) + m(x2_2,x4) + m(x3_2,x3) + m(x8,x8_19) + ((m(x7_2,x9_19)) << 1); + let mut z7 = m(x0_2,x7) + m(x1_2,x6) + m(x2_2,x5) + m(x3_2,x4) + ((m(x8,x9_19)) << 1); + let mut z8 = m(x0_2,x8) + m(x1_2,x7_2) + m(x2_2,x6) + m(x3_2,x5_2) + m(x4,x4) + ((m(x9,x9_19)) << 1); + let mut z9 = m(x0_2,x9) + m(x1_2,x8) + m(x2_2,x7) + m(x3_2,x6) + m(x4_2,x5); + + #[inline(always)] + fn mask_neg(x: u64x4, p: u64x4, mask: u8) -> u64x4 { + unsafe { + use stdsimd::vendor::_mm256_blend_epi32; + _mm256_blend_epi32(x.into(), (p - x).into(), mask as i32).into() + } + } + + // The biggest z_i is bounded as z_i < 249*2^(51 + 2*b); + // if b < 1.5 we get z_i < 4485585228861014016. + // + // The limbs of the multiples of p are bounded above by + // + // 0x3fffffff << 37 = 9223371899415822336 < 2^63 + // + // and below by + // + // 0x1fffffff << 37 = 4611685880988434432 + // > 4485585228861014016 + // + // So these multiples of p are big enough to avoid underflow + // in subtraction, and small enough to fit within u64 + // with room for a carry. + + let low__p37 = u64x4::splat(0x3ffffed << 37); + let even_p37 = u64x4::splat(0x3ffffff << 37); + let odd__p37 = u64x4::splat(0x1ffffff << 37); + + z0 = mask_neg(z0, low__p37, neg_mask); + z1 = mask_neg(z1, odd__p37, neg_mask); + z2 = mask_neg(z2, even_p37, neg_mask); + z3 = mask_neg(z3, odd__p37, neg_mask); + z4 = mask_neg(z4, even_p37, neg_mask); + z5 = mask_neg(z5, odd__p37, neg_mask); + z6 = mask_neg(z6, even_p37, neg_mask); + z7 = mask_neg(z7, odd__p37, neg_mask); + z8 = mask_neg(z8, even_p37, neg_mask); + z9 = mask_neg(z9, odd__p37, neg_mask); + + FieldElement32x4::reduce64([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9]) + } +} + +impl<'a, 'b> Mul<&'b FieldElement32x4> for &'a FieldElement32x4 { + type Output = FieldElement32x4; + fn mul(self, _rhs: &'b FieldElement32x4) -> FieldElement32x4 { + + #[inline(always)] + fn m(x: u32x8, y: u32x8) -> u64x4 { + use stdsimd::vendor::_mm256_mul_epu32; + unsafe { _mm256_mul_epu32(x,y) } + } + + #[inline(always)] + fn m_lo(x: u32x8, y: u32x8) -> u32x8 { + use stdsimd::vendor::_mm256_mul_epu32; + unsafe { u32x8::from(_mm256_mul_epu32(x,y)) } + } + + let (x0, x1) = unpack_pair(self.0[0]); + let (x2, x3) = unpack_pair(self.0[1]); + let (x4, x5) = unpack_pair(self.0[2]); + let (x6, x7) = unpack_pair(self.0[3]); + let (x8, x9) = unpack_pair(self.0[4]); + + let (y0, y1) = unpack_pair(_rhs.0[0]); + let (y2, y3) = unpack_pair(_rhs.0[1]); + let (y4, y5) = unpack_pair(_rhs.0[2]); + let (y6, y7) = unpack_pair(_rhs.0[3]); + let (y8, y9) = unpack_pair(_rhs.0[4]); + + let v19 = u32x8::new(19,0,19,0,19,0,19,0); + + let y1_19 = m_lo(v19, y1); // This fits in a u32 + let y2_19 = m_lo(v19, y2); // iff 26 + b + lg(19) < 32 + let y3_19 = m_lo(v19, y3); // if b < 32 - 26 - 4.248 = 1.752 + let y4_19 = m_lo(v19, y4); + let y5_19 = m_lo(v19, y5); // below, b<2.5: this is a bottleneck, + let y6_19 = m_lo(v19, y6); // could be avoided by promoting to + let y7_19 = m_lo(v19, y7); // u64 here instead of in m() + let y8_19 = m_lo(v19, y8); + let y9_19 = m_lo(v19, y9); + + let x1_2 = x1 + x1; // This fits in a u32 iff 25 + b + 1 < 32 + let x3_2 = x3 + x3; // iff b < 6 + let x5_2 = x5 + x5; + let x7_2 = x7 + x7; + let x9_2 = x9 + x9; + + let z0 = m(x0,y0) + m(x1_2,y9_19) + m(x2,y8_19) + m(x3_2,y7_19) + m(x4,y6_19) + m(x5_2,y5_19) + m(x6,y4_19) + m(x7_2,y3_19) + m(x8,y2_19) + m(x9_2,y1_19); + let z1 = m(x0,y1) + m(x1,y0) + m(x2,y9_19) + m(x3,y8_19) + m(x4,y7_19) + m(x5,y6_19) + m(x6,y5_19) + m(x7,y4_19) + m(x8,y3_19) + m(x9,y2_19); + let z2 = m(x0,y2) + m(x1_2,y1) + m(x2,y0) + m(x3_2,y9_19) + m(x4,y8_19) + m(x5_2,y7_19) + m(x6,y6_19) + m(x7_2,y5_19) + m(x8,y4_19) + m(x9_2,y3_19); + let z3 = m(x0,y3) + m(x1,y2) + m(x2,y1) + m(x3,y0) + m(x4,y9_19) + m(x5,y8_19) + m(x6,y7_19) + m(x7,y6_19) + m(x8,y5_19) + m(x9,y4_19); + let z4 = m(x0,y4) + m(x1_2,y3) + m(x2,y2) + m(x3_2,y1) + m(x4,y0) + m(x5_2,y9_19) + m(x6,y8_19) + m(x7_2,y7_19) + m(x8,y6_19) + m(x9_2,y5_19); + let z5 = m(x0,y5) + m(x1,y4) + m(x2,y3) + m(x3,y2) + m(x4,y1) + m(x5,y0) + m(x6,y9_19) + m(x7,y8_19) + m(x8,y7_19) + m(x9,y6_19); + let z6 = m(x0,y6) + m(x1_2,y5) + m(x2,y4) + m(x3_2,y3) + m(x4,y2) + m(x5_2,y1) + m(x6,y0) + m(x7_2,y9_19) + m(x8,y8_19) + m(x9_2,y7_19); + let z7 = m(x0,y7) + m(x1,y6) + m(x2,y5) + m(x3,y4) + m(x4,y3) + m(x5,y2) + m(x6,y1) + m(x7,y0) + m(x8,y9_19) + m(x9,y8_19); + let z8 = m(x0,y8) + m(x1_2,y7) + m(x2,y6) + m(x3_2,y5) + m(x4,y4) + m(x5_2,y3) + m(x6,y2) + m(x7_2,y1) + m(x8,y0) + m(x9_2,y9_19); + let z9 = m(x0,y9) + m(x1,y8) + m(x2,y7) + m(x3,y6) + m(x4,y5) + m(x5,y4) + m(x6,y3) + m(x7,y2) + m(x8,y1) + m(x9,y0); + + FieldElement32x4::reduce64([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9]) + } +} + + +#[cfg(test)] +mod test { + use super::*; + + #[test] + fn scale_by_curve_constants() { + let mut x = FieldElement32x4::splat(&FieldElement64::one()); + x.scale_by_curve_constants(); + + let xs = x.split(); + assert_eq!(xs[0], FieldElement64([ 121666,0,0,0,0])); + assert_eq!(xs[1], FieldElement64([ 121666,0,0,0,0])); + assert_eq!(xs[2], FieldElement64([2*121666,0,0,0,0])); + assert_eq!(xs[3], FieldElement64([2*121665,0,0,0,0])); + } + + #[test] + fn diff_sum_vs_serial() { + let x0 = FieldElement64([10000, 10001, 10002, 10003, 10004]); + let x1 = FieldElement64([10100, 10101, 10102, 10103, 10104]); + let x2 = FieldElement64([10200, 10201, 10202, 10203, 10204]); + let x3 = FieldElement64([10300, 10301, 10302, 10303, 10304]); + + let mut vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); + vec.diff_sum(0xff); + + let result = vec.split(); + + assert_eq!(result[0], &x1 - &x0); + assert_eq!(result[1], &x1 + &x0); + assert_eq!(result[2], &x3 - &x2); + assert_eq!(result[3], &x3 + &x2); + + let mut vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); + vec.diff_sum(0b01011111); // leave D unchanged + + let result = vec.split(); + + assert_eq!(result[0], &x1 - &x0); + assert_eq!(result[1], &x1 + &x0); + assert_eq!(result[2], &x3 - &x2); + assert_eq!(result[3], x3); + } + + #[test] + fn square_vs_serial() { + let x0 = FieldElement64([10000, 10001, 10002, 10003, 10004]); + let x1 = FieldElement64([10100, 10101, 10102, 10103, 10104]); + let x2 = FieldElement64([10200, 10201, 10202, 10203, 10204]); + let x3 = FieldElement64([10300, 10301, 10302, 10303, 10304]); + + let vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); + + let neg_mask = 0b11_00_00_00; + + let result = vec.square(neg_mask).split(); + + assert_eq!(result[0], &x0 * &x0); + assert_eq!(result[1], &x1 * &x1); + assert_eq!(result[2], &x2 * &x2); + assert_eq!(result[3], -&(&x3 * &x3)); + } + + + #[test] + fn multiply_vs_serial() { + let x0 = FieldElement64([10000, 10001, 10002, 10003, 10004]); + let x1 = FieldElement64([10100, 10101, 10102, 10103, 10104]); + let x2 = FieldElement64([10200, 10201, 10202, 10203, 10204]); + let x3 = FieldElement64([10300, 10301, 10302, 10303, 10304]); + + let vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); + let vecprime = vec.clone(); + + let result = (&vec * &vecprime).split(); + + assert_eq!(result[0], &x0 * &x0); + assert_eq!(result[1], &x1 * &x1); + assert_eq!(result[2], &x2 * &x2); + assert_eq!(result[3], &x3 * &x3); + } + + #[test] + fn test_unpack_repack_pair() { + let x0 = FieldElement64([10000 + (10001 << 26), 0, 0, 0, 0]); + let x1 = FieldElement64([10100 + (10101 << 26), 0, 0, 0, 0]); + let x2 = FieldElement64([10200 + (10201 << 26), 0, 0, 0, 0]); + let x3 = FieldElement64([10300 + (10301 << 26), 0, 0, 0, 0]); + + let vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); + + let src = vec.0[0]; + + let (a,b) = unpack_pair(src); + + let expected_a = u32x8::new(10000, 0, 10100, 0, 10200, 0, 10300, 0); + let expected_b = u32x8::new(10001, 0, 10101, 0, 10201, 0, 10301, 0); + + assert_eq!(a, expected_a); + assert_eq!(b, expected_b); + + let expected_src = repack_pair(a,b); + + assert_eq!(src, expected_src); + } + + #[test] + fn new_split_roundtrips() { + let x0 = FieldElement64::from_bytes(&[0x10; 32]); + let x1 = FieldElement64::from_bytes(&[0x11; 32]); + let x2 = FieldElement64::from_bytes(&[0x12; 32]); + let x3 = FieldElement64::from_bytes(&[0x13; 32]); + + let vec = FieldElement32x4::new(&x0, &x1, &x2, &x3); + + let splits = vec.split(); + + assert_eq!(x0, splits[0]); + assert_eq!(x1, splits[1]); + assert_eq!(x2, splits[2]); + assert_eq!(x3, splits[3]); + } + +} + +#[cfg(all(test, feature = "bench"))] +mod bench { + use test::Bencher; + use super::*; + + #[bench] + fn multiply(b: &mut Bencher) { + let vec = FieldElement32x4::splat(&FieldElement64::zero()); + let vecprime = vec.clone(); + + b.iter(|| &vec * &vecprime ); + } +} + diff --git a/src/backend/avx2/mod.rs b/src/backend/avx2/mod.rs new file mode 100644 index 0000000..04462b1 --- /dev/null +++ b/src/backend/avx2/mod.rs @@ -0,0 +1,487 @@ +// -*- mode: rust; -*- +// +// This file is part of curve25519-dalek. +// Copyright (c) 2016-2017 Isis Lovecruft, Henry de Valence +// See LICENSE for licensing information. +// +// Authors: +// - Isis Agora Lovecruft +// - Henry de Valence + +//! An implementation of group operations on the twisted Edwards form of +//! Curve25519, using AVX2 to implement the 4-way parallel formulas of +//! Hisil, Wong, Carter, and Dawson (HWCD). +//! +//! Their 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08], which +//! introduced the extended coordinates used in other parts of `-dalek`, +//! also describes 4-way parallel formulas for point addition and +//! doubling: +//! +//! * a unified addition algorithm taking an effective \\(2\mathbf M + +//! 1\mathbf D\\); +//! +//! * a doubling algorithm taking an effective \\(1\mathbf M + 1\mathbf +//! S\\); +//! +//! * a dedicated (i.e., for distinct points) addition algorithm taking +//! an effective \\(2 \mathbf M \\). +//! +//! Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of +//! multiplication and squaring of generic field elements and \\(\mathbf +//! D\\) represents the cost of multiplication by a curve constant. +//! +//! Currently, this implementation uses only the first two algorithms. +//! +//! # Parallel formulas +//! +//! The doubling formula is presented in the HWCD paper as follows: +//! +//! | Cost | Processor 1 | Processor 2 | Processor 3 | Processor 4 | +//! |------------------|--------------------------------|--------------------------------|--------------------------------|--------------------------------| +//! | | idle | idle | idle | \\( R\_1 \gets X\_1 + Y\_1 \\) | +//! | \\(1\mathbf S\\) | \\( R\_2 \gets X\_1\^2 \\) | \\( R\_3 \gets Y\_1\^2 \\) | \\( R\_4 \gets Z\_1\^2 \\) | \\( R\_5 \gets R\_1\^2 \\) | +//! | | \\( R\_6 \gets R\_2 + R\_3 \\) | \\( R\_7 \gets R\_2 - R\_3 \\) | \\( R\_4 \gets 2 R\_4 \\) | idle | +//! | | idle | \\( R\_1 \gets R\_4 + R\_7 \\) | idle | \\( R\_2 \gets R\_6 - R\_5 \\) | +//! | \\(1\mathbf M\\) | \\( X\_3 \gets R\_1 R\_2 \\) | \\( Y\_3 \gets R\_6 R\_7 \\) | \\( T\_3 \gets R\_2 R\_6 \\) | \\( Z\_3 \gets R\_1 R\_7 \\) | +//! +//! and the unified addition algorithm is presented as follows: +//! +//! | Cost | Processor 1 | Processor 2 | Processor 3 | Processor 4 | +//! |------------------|--------------------------------|--------------------------------|--------------------------------|--------------------------------| +//! | | \\( R\_1 \gets Y\_1 - X\_1 \\) | \\( R\_2 \gets Y\_2 - X\_2 \\) | \\( R\_3 \gets Y\_1 + X\_1 \\) | \\( R\_4 \gets Y\_2 + X\_2 \\) | +//! | \\(1\mathbf M\\) | \\( R\_5 \gets R\_1 R\_2 \\) | \\( R\_6 \gets R\_3 R\_4 \\) | \\( R\_7 \gets T\_1 T\_2 \\) | \\( R\_8 \gets Z\_1 Z\_2 \\) | +//! | \\(1\mathbf D\\) | idle | idle | \\( R\_7 \gets k R\_7 \\) | \\( R\_8 \gets 2 R\_8 \\) | +//! | | \\( R\_1 \gets R\_6 - R\_5 \\) | \\( R\_2 \gets R\_8 - R\_7 \\) | \\( R\_3 \gets R\_8 + R\_7 \\) | \\( R\_4 \gets R\_6 + R\_5 \\) | +//! | \\(1\mathbf M\\) | \\( X\_3 \gets R\_1 R\_2 \\) | \\( Y\_3 \gets R\_3 R\_4 \\) | \\( T\_3 \gets R\_1 R\_4 \\) | \\( Z\_3 \gets R\_2 R\_3 \\) | +//! +//! Here \\( k = 2d \\) is a curve constant. +//! +//! # Implementation strategy +//! +//! For a software implementation, each "processor"'s operations are too +//! low-latency to parallelize across threads. However, the main cost +//! is in the multiplication and squaring steps, which share a single +//! instruction. +//! +//! Our strategy is to implement 4-wide multiplication and squaring +//! using one 64-bit AVX2 lane for each field element. Field elements +//! are represented in the usual way as 10 `u32` limbs in radix +//! \\(25.5\\) (i.e., alternating between \\(2\^{26}\\) for even limbs +//! and \\(2\^{25}\\) for odd limbs). This has the effect that passing +//! between the parallel 32-bit AVX2 representation and the serial +//! 64-bit representation amounts to regrouping digits. +//! +//! The addition and subtraction steps are done largely serially, using +//! masking to handle the instruction divergence. The remaining +//! obstacle to parallelism is the multiplication by the curve constant +//! \\(k = 2d\\). In the Curve25519 case, this is +//! +//! $$ k \equiv 2 \frac{-121665}{121666} \\ \equiv 16295367250680780974490674513165176452449235426866156013048779062215315747161 \pmod p. $$ +//! +//! HWCD suggest parallelising this step by breaking \\(k\\) into four +//! parts as \\(k = k_0 + 2\^n k_1 + 2\^{2n} k_2 + 2\^{3n} k_3 \\) and +//! computing \\(k_i R_7 \\) in parallel. However, this would be +//! somewhat awkward in our case, since we would normally represent +//! \\(k\\) as \\( 10 \\) 32-bit limbs, and \\(10 \\) is not divisible +//! by \\(4\\), so we would need a specialized routine to perform a +//! vectorized multiplication by 64-bit constants. +//! +//! Instead, since we are working projectively, we can multiply +//! \\(R_7\\) by \\( -2\cdot 121665 \\) and multiply the other three +//! variables by \\(121666\\). This trick was suggested by Mike +//! Hamburg. Ignoring the sign for the moment, since +//! \\(2 \cdot 121666 < 2\^{18}\\), all these constants fit in 32 bits, +//! so (up to sign) this can be done in parallel as four multiplications +//! by small constants \\( (121666, 121666, 2\cdot 121665, 2\cdot 121666) \\). +//! +//! How do we handle the sign? +//! Since we're primarily interested in Ristretto performance, not +//! Curve25519 performance, we could alternately work on the +//! \\(4\\)-isogenous "IsoEd25519" curve, which has \\(d = 121665\\). +//! However, this would only save the negation step, since multiplying +//! one field element by a 32-bit constant is not much easier than +//! multiplying four field elements by 32-bit constants, and it would +//! prevent accelerating Curve25519, so we don't make this choice. +//! Instead, we just negate one lane, and move the \\(1 \mathbf D\\) +//! into precomputation (see below). +//! +//! The 4-wide formulas of the HWCD paper do not seem to have been +//! implemented using SIMD before. The HWCD paper also describes and +//! analyzes a 2-wide variant of the Montgomery ladder (for comparison +//! with parallel Edwards formulas); this strategy was used in 2015 by +//! Tung Chou's `sandy2x` implementation, which used a 2-wide field +//! implementation in 128-bit vector registers. +//! +//! Curiously, however, although the [`sandy2x` paper][sandy2x] also +//! implements Edwards arithmetic, and cites the HWCD paper, it doesn't +//! mention or discuss the parallel formulas from HWCD, or that the +//! 2-wide Montgomery formulas it uses were previously published there. +//! There is also a 2015 paper by Hernández and López on using AVX2 for +//! the X25519 Montgomery ladder, but neither the paper nor the code are +//! publicly available, and it apparently gives only a [slight +//! speedup][avx2trac], suggesting that it also overlooked the +//! HWCD formulas. +//! +//! HWCD also suggest using a mixed representation, passing between \\( +//! \mathbb P\^3 \\) "extended" coordinates and \\( \mathbb P\^2 \\) +//! "projective" coordinates, where doubling is slightly cheaper (saving +//! about \\(\mathbf 1M\\). This approach is used for the +//! non-vectorized `u32` and `u64` backends, and more +//! details on the different coordinate systems can be found in the +//! `curve_models` module documentation. +//! +//! This optimization is not compatible with the parallel formulas, which are +//! therefore slightly less efficient when counting the total number of +//! field multiplications and squarings. In particular, vectorized doublings +//! are less efficient than serial doublings. +//! In addition, the parallel formulas can only use a \\( 32 \times 32 +//! \rightarrow 64 \\)-bit integer multiplier, so the speedup from +//! vectorization must overcome the disadvantage of losing the \\( 64 +//! \times 64 \rightarrow 128\\)-bit (serial) integer multiplier. +//! +//! # Tweaked formulas +//! +//! After tweaking the formulas as described above, we obtain the +//! following. To avoid confusion with the original HWCD formulas, +//! temporary variables are named \\(S\\) instead of \\(R\\) and are in +//! static single-assignment (SSA) form. +//! +//! ## Addition +//! +//! To add points \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and \\(P_2 = (X_2 +//! : Y_2 : Z_2 : T_2 ) \\), we compute +//! +//! $$ +//! \begin{aligned} +//! S\_0 &\gets Y\_1 - X\_1 \\\\ +//! S\_1 &\gets Y\_1 + X\_1 \\\\ +//! S\_2 &\gets Y\_2 - X\_2 \\\\ +//! S\_3 &\gets Y\_2 + X\_2 +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! S\_4 &\gets S\_0 S\_2 \\\\ +//! S\_5 &\gets S\_1 S\_3 \\\\ +//! S\_6 &\gets Z\_1 Z\_2 \\\\ +//! S\_7 &\gets T\_1 T\_2 +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! S\_8 &\gets S\_4 \cdot 121666 \\\\ +//! S\_9 &\gets S\_5 \cdot 121666 \\\\ +//! S\_{10} &\gets S\_6 \cdot 2 \cdot 121666 \\\\ +//! S\_{11} &\gets S\_7 \cdot -2 \cdot 121665 +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! S\_{12} &\gets S\_9 - S\_8 \\\\ +//! S\_{13} &\gets S\_9 + S\_8 \\\\ +//! S\_{14} &\gets S\_{10} - S\_{11} \\\\ +//! S\_{15} &\gets S\_{10} + S\_{11} +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! X\_3 &\gets S\_{12} S\_{14} \\\\ +//! Y\_3 &\gets S\_{15} S\_{13} \\\\ +//! Z\_3 &\gets S\_{15} S\_{14} \\\\ +//! T\_3 &\gets S\_{12} S\_{13} +//! \end{aligned} +//! $$ +//! +//! to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\). +//! +//! ## Readdition +//! +//! If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we can precompute +//! +//! $$ +//! \begin{aligned} +//! S\_2 &\gets Y\_2 - X\_2 \\\\ +//! S\_3 &\gets Y\_2 + X\_2 +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! S\_2' &\gets S\_2 \cdot 121666 \\\\ +//! S\_3' &\gets S\_3 \cdot 121666 \\\\ +//! Z\_2' &\gets Z\_2 \cdot 2 \cdot 121666 \\\\ +//! T\_2' &\gets T\_2 \cdot -2 \cdot 121665 \\\\ +//! \end{aligned} +//! $$ +//! +//! to obtain the `CachedPoint` \\( (S\_2', S\_3', Z\_2', T\_2') \\). +//! This precomputation is essentially the same as that suggested in +//! §3.1 of HWCD, with the difference that the multiplication by the curve +//! constant \\( -121665 / 121666 \\) is spread over all four +//! coordinates, to allow a vectorized computation of four +//! multiplications of small constants instead of a serial computation +//! of multiplication by a large constant. +//! +//! To perform readdition of \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and +//! \\(P_2 = (S\_2', S\_3', Z\_2', T\_2') \\), we compute +//! +//! $$ +//! \begin{aligned} +//! S\_0 &\gets Y\_1 - X\_1 \\\\ +//! S\_1 &\gets Y\_1 + X\_1 +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! S\_8 &\gets S\_0 S\_2' \\\\ +//! S\_9 &\gets S\_1 S\_3' \\\\ +//! S\_{10} &\gets Z\_1 Z\_2' \\\\ +//! S\_{11} &\gets T\_1 T\_2' +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! S\_{12} &\gets S\_9 - S\_8 \\\\ +//! S\_{13} &\gets S\_9 + S\_8 \\\\ +//! S\_{14} &\gets S\_{10} - S\_{11} \\\\ +//! S\_{15} &\gets S\_{10} + S\_{11} +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! X\_3 &\gets S\_{12} S\_{14} \\\\ +//! Y\_3 &\gets S\_{15} S\_{13} \\\\ +//! Z\_3 &\gets S\_{15} S\_{14} \\\\ +//! T\_3 &\gets S\_{12} S\_{13} +//! \end{aligned} +//! $$ +//! +//! to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\). +//! +//! Compared to the addition formulas above, this saves \\( 1\mathbf D \\). +//! +//! ## Doubling +//! +//! To double a point \\( P = (X\_1 : Y\_1 : Z\_1 : T\_1) \\), we compute +//! +//! $$ S\_0 \gets X\_1 + Y\_1 $$ +//! +//! $$ +//! \begin{aligned} +//! S\_1 &\gets X\_1\^2 \\\\ +//! S\_2 &\gets Y\_1\^2 \\\\ +//! S\_3 &\gets Z\_1\^2 \\\\ +//! S\_4 &\gets S\_0\^2 +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! S\_5 &\gets S\_1 + S\_2 \\\\ +//! S\_6 &\gets S\_1 - S\_2 \\\\ +//! S\_7 &\gets 2S\_3 \\\\ +//! S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\ +//! S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4 +//! \end{aligned} +//! $$ +//! +//! $$ +//! \begin{aligned} +//! X\_3 &\gets S\_8 S\_9 \\\\ +//! Y\_3 &\gets S\_5 S\_6 \\\\ +//! Z\_3 &\gets S\_8 S\_6 \\\\ +//! T\_3 &\gets S\_5 S\_9 +//! \end{aligned} +//! $$ +//! +//! to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = [2]P\_1 \\). +//! +//! Performing too many intermediate additions and subtractions grows +//! the bounds beyond what is allowed as input to multiplication, +//! forcing an extra carry pass. However, it is just possible to avoid +//! this by rearranging signs. +//! +//! Assume that the bounds on the limbs of each field element are +//! parameterized by \\( b \in \mathbb R \\) representing the excess +//! bits, so that each limb is bounded by either \\( 2\^{25} \\) or \\( +//! 2\^{26} \\). +//! +//! The multiplication routine requires that its inputs are bounded by +//! \\( b < 1.75 \\), in order to fit a multiplication by \\( 19 \\) +//! into 32 bits. Since \\( \lg 19 < 4.25 \\), \\( 19x < 2\^{32} \\) +//! when \\( x < 2\^{27.75} = 2\^{26 + 1.75} \\). However, this is only +//! required for one of the inputs; the other can grow up to \\( b < 2.5 +//! \\). +//! +//! Computing \\( (S\_5, S\_6, S\_8, S\_9 ) \\) as +//! +//! $$ +//! \begin{matrix} +//! & S\_1 & S\_1 & S\_1 & S\_1 \\\\ +//! +& S\_2 & & & S\_2 \\\\ +//! +& & & S\_3 & \\\\ +//! +& & & S\_3 & \\\\ +//! +& & 2p & 2p & 2p \\\\ +//! -& & S\_2 & S\_2 & \\\\ +//! -& & & & S\_4 \\\\ +//! =& S\_5 & S\_6 & S\_8 & S\_9 +//! \end{matrix} +//! $$ +//! +//! results in bit-excesses \\( (1.00, 1.59, 2.33, 2.00)\\) for +//! \\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute +//! are then +//! +//! $$ +//! \begin{aligned} +//! X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.00) \\\\ +//! Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ +//! Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ +//! T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 2.00) +//! \end{aligned} +//! $$ +//! +//! which are too large. However, if we flip the sign of \\( S\_4 = +//! S\_0\^2 \\) during squaring, so that we output \\(S\_4' = -S\_4 +//! \pmod p\\), then we can compute +//! +//! $$ +//! \begin{matrix} +//! & S\_1 & S\_1 & S\_1 & S\_1 \\\\ +//! +& S\_2 & & & S\_2 \\\\ +//! +& & & S\_3 & \\\\ +//! +& & & S\_3 & \\\\ +//! +& & & & S\_4' \\\\ +//! +& & 2p & 2p & \\\\ +//! -& & S\_2 & S\_2 & \\\\ +//! =& S\_5 & S\_6 & S\_8 & S\_9 +//! \end{matrix} +//! $$ +//! +//! resulting in bit-excesses \\( (1.00, 1.59, 2.33, 1.59)\\) for +//! \\( (S\_5, S\_6, S\_8, S\_9 ) \\). The products we want to compute +//! are then +//! +//! $$ +//! \begin{aligned} +//! X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.59) \\\\ +//! Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\ +//! Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\ +//! T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 1.59) +//! \end{aligned} +//! $$ +//! +//! whose right-hand sides are all bounded with \\( b < 1.75 \\) and +//! whose left-hand sides are all bounded with \\( b < 2.5 \\). +//! +//! # Field element representation +//! +//! The field element representation is oriented around the AVX2 +//! `vpmuluqdq` instruction, which multiplies the low 32 bits of each +//! 64-bit lane of each operand to produce a 64-bit result. +//! +//! ```text,no_run +//! (a1 ?? b1 ?? c1 ?? d1 ??) +//! (a2 ?? b2 ?? c2 ?? d2 ??) +//! +//! (a1*a2 b1*b2 c1*c2 d1*d2) +//! ``` +//! +//! To unpack 32-bit values into 64-bit lanes for use in multiplication +//! it would be convenient to use the `vpunpck[lh]dq` instructions, +//! which unpack and interleave the low and high 32-bit lanes of two +//! source vectors. +//! However, the AVX2 versions of these instructions are designed to +//! operate only within 128-bit lanes of the 256-bit vectors, so that +//! interleaving the low lanes of `(a0 b0 c0 d0 a1 b1 c1 d1)` with zero +//! gives `(a0 00 b0 00 a1 00 b1 00)`. Instead, we pre-shuffle the data +//! layout as `(a0 b0 a1 b1 c0 d0 c1 d1)` so that we can unpack the +//! "low" and "high" parts as +//! +//! ```text,no_run +//! (a0 00 b0 00 c0 00 d0 00) +//! (a1 00 b1 00 c1 00 d1 00) +//! ``` +//! +//! The data layout for a vector of four field elements \\( (a,b,c,d) +//! \\) with limbs \\( a_0, a_1, \ldots, a_9 \\) is as `[u32x8; 5]` in +//! the form +//! +//! ```text,no_run +//! (a0 b0 a1 b1 c0 d0 c1 d1) +//! (a2 b2 a3 b3 c2 d2 c3 d3) +//! (a4 b4 a5 b5 c4 d4 c5 d5) +//! (a6 b6 a7 b7 c6 d6 c7 d7) +//! (a8 b8 a9 b9 c8 d8 c9 d9) +//! ``` +//! +//! Since this breaks cleanly into two 128-bit lanes, it may be possible +//! to adapt it to 128-bit vector instructions such as NEON without too +//! much difficulty. +//! +//! Going the other direction, to extend this to AVX512, we could either +//! run two point operations in parallel in lower and upper halves of +//! the registers, or use 2-way parallelism within a field operation. +//! +//! We don't attempt to use AVX2 for serial field element computations +//! such as inversion, since wherever we have AVX2 we also have `mulx`. +//! However, it might be useful for batched inverse square-root +//! computations, which can't be batched in the same way inversions can. +//! +//! # Implementation details +//! +//! The implementation uses the unstable `stdsimd` crate to provide AVX2 +//! intrinsics, and the code is not yet cleanly factored between the +//! field element parts and the point parts. +//! +//! When compiling with AVX512VL, LLVM is able to use the extra +//! `ymm16..ymm31` registers to reduce register pressure, and avoid +//! spills during field multiplication. This gives a small but +//! noticeable speedup. +//! +//! The addition and subtraction steps involve masking, to apply +//! operations to a single lane of the vector. AVX512VL extends the +//! predication features of AVX512 to AVX2 code and would probably be +//! beneficial. Unfortunately, LLVM is currently unable to lower `op + +//! blend` into an AVX512VL masked operation. However, the explicitly +//! masked versions of the intrinsics seem to produce the same LLVM IR +//! as an `op + blend`, so hopefully this will improve as the AVX512 +//! support in LLVM improves. +//! +//! When used for constant-time variable-base scalar multiplication, +//! this strategy (using AVX2) gives a significant speedup over the +//! serial implementation (using the \\(64 \times 64\\) multiplier) of +//! approximately 1.6x for Skylake-X with `target_cpu=skylake` (using AVX2), of +//! approximately 1.8x for Skylake-X with `target_cpu=skylake-avx512` (using the extra +//! `ymm16..ymm31` registers from AVX512VL), and of approximately 1.0x +//! for Ryzen (which implements AVX2 at half rate). +//! +//! When used for variable-time double-base scalar multiplication \\( aA +//! + bB \\) for fixed \\(B\\) (as in, e.g., signature verification), +//! this strategy provides a 1.4x speedup on Skylake-X over the same +//! operation as implemented in `ed25519-donna`, the fastest +//! production-quality Ed25519 implementation. +//! +//! (Note: since testing this, the experimental `llvm50` Rust branch +//! used to compile the experimental `stdsimd` intrinsics have fallen +//! out of sync and it is no longer possible to compile for +//! `skylake-avx512`. This is why all of this branch is part of the +//! `yolocrypto` feature, pending upstream work.) +//! +//! [sandy2x]: https://eprint.iacr.org/2015/943.pdf +//! [avx2trac]: https://trac.torproject.org/projects/tor/ticket/8897#comment:28 +//! [hwcd08]: https://www.iacr.org/archive/asiacrypt2008/53500329/53500329.pdf + + +pub(crate) mod field; + +pub(crate) mod edwards; + +pub(crate) mod constants; diff --git a/src/backend/mod.rs b/src/backend/mod.rs index be3a347..0e2ce3e 100644 --- a/src/backend/mod.rs +++ b/src/backend/mod.rs @@ -28,3 +28,7 @@ pub mod u32; #[cfg(feature="radix_51")] pub mod u64; +/// Code using AVX2. +#[cfg(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2")))] +pub mod avx2; + diff --git a/src/edwards.rs b/src/edwards.rs index 1f00d9b..5fca4e9 100644 --- a/src/edwards.rs +++ b/src/edwards.rs @@ -432,32 +432,41 @@ impl<'a, 'b> Mul<&'b Scalar> for &'a ExtendedPoint { /// For scalar multiplication of a basepoint, /// `EdwardsBasepointTable` is approximately 4x faster. fn mul(self, scalar: &'b Scalar) -> ExtendedPoint { - // Construct a lookup table of [P,2P,3P,4P,5P,6P,7P,8P] - let lookup_table = LookupTable::::from(self); - - // Setting s = scalar, compute - // - // s = s_0 + s_1*16^1 + ... + s_63*16^63, - // - // with `-8 ≤ s_i < 8` for `0 ≤ i < 63` and `-8 ≤ s_63 ≤ 8`. - let scalar_digits = scalar.to_radix_16(); - - // Compute s*P as - // - // s*P = P*(s_0 + s_1*16^1 + s_2*16^2 + ... + s_63*16^63) - // s*P = P*s_0 + P*s_1*16^1 + P*s_2*16^2 + ... + P*s_63*16^63 - // 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 = ExtendedPoint::identity(); - for i in (0..64).rev() { - // Q <-- 16*Q - Q = Q.mult_by_pow_2(4); - // Q <-- Q + P * s_i - Q = (&Q + &lookup_table.select(scalar_digits[i])).to_extended() + // If we built with AVX2, use the AVX2 backend. + #[cfg(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2")))] { + use backend::avx2::edwards as edwards_avx2; + let P_avx2 = edwards_avx2::ExtendedPoint::from(*self); + return ExtendedPoint::from(&P_avx2 * scalar); } + // Otherwise, proceed as normal: + #[cfg(not(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2"))))] { + // Construct a lookup table of [P,2P,3P,4P,5P,6P,7P,8P] + let lookup_table = LookupTable::::from(self); - Q + // Setting s = scalar, compute + // + // s = s_0 + s_1*16^1 + ... + s_63*16^63, + // + // with `-8 ≤ s_i < 8` for `0 ≤ i < 63` and `-8 ≤ s_63 ≤ 8`. + let scalar_digits = scalar.to_radix_16(); + + // Compute s*P as + // + // s*P = P*(s_0 + s_1*16^1 + s_2*16^2 + ... + s_63*16^63) + // s*P = P*s_0 + P*s_1*16^1 + P*s_2*16^2 + ... + P*s_63*16^63 + // 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 = ExtendedPoint::identity(); + for i in (0..64).rev() { + // Q <-- 16*Q + Q = Q.mult_by_pow_2(4); + // Q <-- Q + P * s_i + Q = (&Q + &lookup_table.select(scalar_digits[i])).to_extended() + } + + Q + } } } @@ -487,7 +496,6 @@ impl<'a, 'b> Mul<&'b ExtendedPoint> for &'a Scalar { /// A iterable of `Scalar`s and a iterable of `ExtendedPoints`. It is an /// error to call this function with two iterators of different lengths. /// -/// XXX need to clear memory // XXX later when we do more fancy multiscalar mults, we can delegate // based on the iter's size hint -- hdevalence #[cfg(any(feature = "alloc", feature = "std"))] @@ -495,62 +503,71 @@ pub fn multiscalar_mult<'a, 'b, I, J>(scalars: I, points: J) -> ExtendedPoint where I: IntoIterator, J: IntoIterator { - //assert_eq!(scalars.len(), points.len()); - - use clear_on_drop::ClearOnDrop; + // If we built with AVX2, use the AVX2 backend. + #[cfg(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2")))] { + use backend::avx2::edwards as edwards_avx2; - let lookup_tables_vec: Vec<_> = points.into_iter() - .map(|P| LookupTable::::from(P) ) - .collect(); - - let lookup_tables = ClearOnDrop::new(lookup_tables_vec); - - // Setting s_i = i-th scalar, compute - // - // s_i = s_{i,0} + s_{i,1}*16^1 + ... + s_{i,63}*16^63, - // - // with `-8 ≤ s_{i,j} < 8` for `0 ≤ j < 63` and `-8 ≤ s_{i,63} ≤ 8`. - let scalar_digits_vec: Vec<_> = scalars.into_iter() - .map(|c| c.to_radix_16()) - .collect(); - - // This above puts the scalar digits into a heap-allocated Vec. - // To ensure that these are erased, pass ownership of the Vec into a - // ClearOnDrop wrapper. - let scalar_digits = ClearOnDrop::new(scalar_digits_vec); - - // Compute s_1*P_1 + ... + s_n*P_n: since - // - // s_i*P_i = P_i*(s_{i,0} + s_{i,1}*16^1 + ... + s_{i,63}*16^63) - // s_i*P_i = P_i*s_{i,0} + P_i*s_{i,1}*16^1 + ... + P_i*s_{i,63}*16^63 - // s_i*P_i = P_i*s_{i,0} + 16*(P_i*s_{i,1} + 16*( ... + 16*P_i*s_{i,63})...) - // - // we have the two-dimensional sum - // - // s_1*P_1 = P_1*s_{1,0} + 16*(P_1*s_{1,1} + 16*( ... + 16*P_1*s_{1,63})...) - // + s_2*P_2 = + P_2*s_{2,0} + 16*(P_2*s_{2,1} + 16*( ... + 16*P_2*s_{2,63})...) - // ... - // + s_n*P_n = + P_n*s_{n,0} + 16*(P_n*s_{n,1} + 16*( ... + 16*P_n*s_{n,63})...) - // - // We sum column-wise top-to-bottom, then right-to-left, - // multiplying by 16 only once per column. - // - // This provides the speedup over doing n independent scalar - // mults: we perform 63 multiplications by 16 instead of 63*n - // multiplications, saving 252*(n-1) doublings. - let mut Q = ExtendedPoint::identity(); - // XXX this impl makes no effort to be cache-aware; maybe it could be improved? - for j in (0..64).rev() { - Q = Q.mult_by_pow_2(4); - let it = scalar_digits.iter().zip(lookup_tables.iter()); - for (s_i, lookup_table_i) in it { - // R_i = s_{i,j} * P_i - let R_i = lookup_table_i.select(s_i[j]); - // Q = Q + R_i - Q = (&Q + &R_i).to_extended(); - } + edwards_avx2::multiscalar_mult(scalars, points) + } + // Otherwise, proceed as normal: + #[cfg(not(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2"))))] { + //assert_eq!(scalars.len(), points.len()); + + use clear_on_drop::ClearOnDrop; + + let lookup_tables_vec: Vec<_> = points.into_iter() + .map(|P| LookupTable::::from(P) ) + .collect(); + + let lookup_tables = ClearOnDrop::new(lookup_tables_vec); + + // Setting s_i = i-th scalar, compute + // + // s_i = s_{i,0} + s_{i,1}*16^1 + ... + s_{i,63}*16^63, + // + // with `-8 ≤ s_{i,j} < 8` for `0 ≤ j < 63` and `-8 ≤ s_{i,63} ≤ 8`. + let scalar_digits_vec: Vec<_> = scalars.into_iter() + .map(|c| c.to_radix_16()) + .collect(); + + // This above puts the scalar digits into a heap-allocated Vec. + // To ensure that these are erased, pass ownership of the Vec into a + // ClearOnDrop wrapper. + let scalar_digits = ClearOnDrop::new(scalar_digits_vec); + + // Compute s_1*P_1 + ... + s_n*P_n: since + // + // s_i*P_i = P_i*(s_{i,0} + s_{i,1}*16^1 + ... + s_{i,63}*16^63) + // s_i*P_i = P_i*s_{i,0} + P_i*s_{i,1}*16^1 + ... + P_i*s_{i,63}*16^63 + // s_i*P_i = P_i*s_{i,0} + 16*(P_i*s_{i,1} + 16*( ... + 16*P_i*s_{i,63})...) + // + // we have the two-dimensional sum + // + // s_1*P_1 = P_1*s_{1,0} + 16*(P_1*s_{1,1} + 16*( ... + 16*P_1*s_{1,63})...) + // + s_2*P_2 = + P_2*s_{2,0} + 16*(P_2*s_{2,1} + 16*( ... + 16*P_2*s_{2,63})...) + // ... + // + s_n*P_n = + P_n*s_{n,0} + 16*(P_n*s_{n,1} + 16*( ... + 16*P_n*s_{n,63})...) + // + // We sum column-wise top-to-bottom, then right-to-left, + // multiplying by 16 only once per column. + // + // This provides the speedup over doing n independent scalar + // mults: we perform 63 multiplications by 16 instead of 63*n + // multiplications, saving 252*(n-1) doublings. + let mut Q = ExtendedPoint::identity(); + // XXX this impl makes no effort to be cache-aware; maybe it could be improved? + for j in (0..64).rev() { + Q = Q.mult_by_pow_2(4); + let it = scalar_digits.iter().zip(lookup_tables.iter()); + for (s_i, lookup_table_i) in it { + // R_i = s_{i,j} * P_i + let R_i = lookup_table_i.select(s_i[j]); + // Q = Q + R_i + Q = (&Q + &R_i).to_extended(); + } + } + Q } - Q } /// A precomputed table of multiples of a basepoint, for accelerating @@ -787,79 +804,99 @@ pub mod vartime { where I: IntoIterator, J: IntoIterator { - //assert_eq!(scalars.len(), points.len()); + // If we built with AVX2, use the AVX2 backend. + #[cfg(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2")))] { + use backend::avx2::edwards as edwards_avx2; - let nafs: Vec<_> = scalars.into_iter() - .map(|c| c.non_adjacent_form()).collect(); - let odd_multiples: Vec<_> = points.into_iter() - .map(|P| OddMultiples::create(P)).collect(); + edwards_avx2::vartime::multiscalar_mult(scalars, points) + } + // Otherwise, proceed as normal: + #[cfg(not(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2"))))] { + //assert_eq!(scalars.len(), points.len()); - let mut r = ProjectivePoint::identity(); + let nafs: Vec<_> = scalars.into_iter() + .map(|c| c.non_adjacent_form()).collect(); + let odd_multiples: Vec<_> = points.into_iter() + .map(|P| OddMultiples::create(P)).collect(); - for i in (0..255).rev() { - let mut t = r.double(); + let mut r = ProjectivePoint::identity(); - for (naf, odd_multiple) in nafs.iter().zip(odd_multiples.iter()) { - if naf[i] > 0 { - t = &t.to_extended() + &odd_multiple[( naf[i]/2) as usize]; - } else if naf[i] < 0 { - t = &t.to_extended() - &odd_multiple[(-naf[i]/2) as usize]; + for i in (0..255).rev() { + let mut t = r.double(); + + for (naf, odd_multiple) in nafs.iter().zip(odd_multiples.iter()) { + if naf[i] > 0 { + t = &t.to_extended() + &odd_multiple[( naf[i]/2) as usize]; + } else if naf[i] < 0 { + t = &t.to_extended() - &odd_multiple[(-naf[i]/2) as usize]; + } } + + r = t.to_projective(); } - r = t.to_projective(); + r.to_extended() } - - r.to_extended() } /// Given a point \\(A\\) and scalars \\(a\\) and \\(b\\), compute the point /// \\(aA+bB\\), where \\(B\\) is the Ed25519 basepoint (i.e., \\(B = (x,4/5)\\) /// with x positive). #[cfg(feature="precomputed_tables")] - pub fn double_scalar_mult_basepoint(a: &Scalar, - A: &ExtendedPoint, - b: &Scalar) -> ExtendedPoint { - let a_naf = a.non_adjacent_form(); - let b_naf = b.non_adjacent_form(); + pub fn double_scalar_mult_basepoint( + a: &Scalar, + A: &ExtendedPoint, + b: &Scalar, + ) -> ExtendedPoint { + // If we built with AVX2, use the AVX2 backend. + #[cfg(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2")))] { + use backend::avx2::edwards as edwards_avx2; - // Find starting index - let mut i: usize = 255; - for j in (0..255).rev() { - i = j; - if a_naf[i] != 0 || b_naf[i] != 0 { - break; - } + edwards_avx2::vartime::double_scalar_mult_basepoint(a, A, b) } + // Otherwise, proceed as normal: + #[cfg(not(all(feature="nightly", all(feature="avx2_backend", target_feature="avx2"))))] { + let a_naf = a.non_adjacent_form(); + let b_naf = b.non_adjacent_form(); - let odd_multiples_of_A = OddMultiples::create(A); - let odd_multiples_of_B = &constants::AFFINE_ODD_MULTIPLES_OF_BASEPOINT; - - let mut r = ProjectivePoint::identity(); - loop { - let mut t = r.double(); - - if a_naf[i] > 0 { - t = &t.to_extended() + &odd_multiples_of_A[( a_naf[i]/2) as usize]; - } else if a_naf[i] < 0 { - t = &t.to_extended() - &odd_multiples_of_A[(-a_naf[i]/2) as usize]; + // Find starting index + let mut i: usize = 255; + for j in (0..255).rev() { + i = j; + if a_naf[i] != 0 || b_naf[i] != 0 { + break; + } } - if b_naf[i] > 0 { - t = &t.to_extended() + &odd_multiples_of_B[( b_naf[i]/2) as usize]; - } else if b_naf[i] < 0 { - t = &t.to_extended() - &odd_multiples_of_B[(-b_naf[i]/2) as usize]; + let odd_multiples_of_A = OddMultiples::create(A); + let odd_multiples_of_B = &constants::AFFINE_ODD_MULTIPLES_OF_BASEPOINT; + + let mut r = ProjectivePoint::identity(); + loop { + let mut t = r.double(); + + if a_naf[i] > 0 { + t = &t.to_extended() + &odd_multiples_of_A[( a_naf[i]/2) as usize]; + } else if a_naf[i] < 0 { + t = &t.to_extended() - &odd_multiples_of_A[(-a_naf[i]/2) as usize]; + } + + if b_naf[i] > 0 { + t = &t.to_extended() + &odd_multiples_of_B[( b_naf[i]/2) as usize]; + } else if b_naf[i] < 0 { + t = &t.to_extended() - &odd_multiples_of_B[(-b_naf[i]/2) as usize]; + } + + r = t.to_projective(); + + if i == 0 { + break; + } + i -= 1; } - r = t.to_projective(); - - if i == 0 { - break; - } - i -= 1; + r.to_extended() } - - r.to_extended() } } @@ -1049,7 +1086,7 @@ mod test { /// Test precomputed basepoint mult #[test] - #[cfg(feature="basepoint_table_creation")] + #[cfg(feature="precomputed_tables")] fn test_precomputed_basepoint_mult() { let table = EdwardsBasepointTable::create(&constants::ED25519_BASEPOINT_POINT); let aB_1 = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR; @@ -1322,14 +1359,14 @@ mod bench { } #[bench] - #[cfg(feature="basepoint_table_creation")] + #[cfg(feature="precomputed_tables")] fn create_basepoint_table(b: &mut Bencher) { let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR; b.iter(|| EdwardsBasepointTable::create(&aB)); } #[bench] - #[cfg(feature="basepoint_table_creation")] + #[cfg(feature="precomputed_tables")] fn ten_fold_scalar_mult(b: &mut Bencher) { let mut csprng: OsRng = OsRng::new().unwrap(); // Create 10 random scalars @@ -1353,7 +1390,7 @@ mod bench { } #[bench] - #[cfg(feature="basepoint_table_creation")] + #[cfg(feature="precomputed_tables")] fn ten_fold_scalar_mult(b: &mut Bencher) { let mut csprng: OsRng = OsRng::new().unwrap(); // Create 10 random scalars diff --git a/src/lib.rs b/src/lib.rs index 51364f0..6c4ed32 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -11,10 +11,9 @@ #![cfg_attr(not(feature = "std"), no_std)] #![cfg_attr(feature = "alloc", feature(alloc))] #![cfg_attr(feature = "nightly", feature(i128_type))] +#![cfg_attr(feature = "nightly", feature(cfg_target_feature))] #![cfg_attr(feature = "bench", feature(test))] -#![cfg_attr(all(feature = "nightly", feature = "std"), feature(zero_one))] -#![allow(unused_features)] #![deny(missing_docs)] // refuse to compile if documentation is missing //! # curve25519-dalek @@ -52,6 +51,9 @@ extern crate clear_on_drop; #[cfg(all(test, feature = "bench"))] extern crate test; +#[cfg(feature = "yolocrypto")] +extern crate stdsimd; + // The `Digest` trait is implemented using `generic_array`, so we need it // too. Hopefully we can eliminate `generic_array` from `Digest` once const // generics land.