saymrwulf/risc0-curve25519-dalek-source
1
0
Fork 0
mirror of https://github.com/saymrwulf/risc0-curve25519-dalek-source.git synced 2026-09-04 20:03:40 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'release/0.14.2'

Browse source
This commit is contained in:
Isis Lovecruft 2018-01-20 00:02:02 +00:00
commit 1fafc2c41e
Failed to extract signature
14 changed files with 2587 additions and 238 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

15
.travis.yml
View file

@ -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

16
Cargo.toml
View file

@ -1,6 +1,6 @@
[package]
name = "curve25519-dalek"
version = "0.14.1"
version = "0.14.2"
authors = ["Isis Lovecruft <isis@patternsinthevoid.net>",
"Henry de Valence <hdevalence@hdevalence.ca>"]
readme = "README.md"
@ -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"]

5
README.md
View file

@ -1,6 +1,11 @@
# curve25519-dalek [![](https://img.shields.io/crates/v/curve25519-dalek.svg)](https://crates.io/crates/curve25519-dalek) [![](https://docs.rs/curve25519-dalek/badge.svg)](https://docs.rs/curve25519-dalek) [![](https://travis-ci.org/isislovecruft/curve25519-dalek.svg?branch=master)](https://travis-ci.org/isislovecruft/curve25519-dalek)
<img
width="50%"
align="right"
src="https://user-images.githubusercontent.com/797/34898472-83686016-f7f3-11e7-967b-24b2aadd623a.png"/>
**A low-level cryptographic library for point, group, field, and scalar
operations on a curve isomorphic to the twisted Edwards curve defined by -x²+y²
= 1 - 121665/121666 x²y² over GF(2²⁵⁵ - 19).**

3
build.rs
View file

@ -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

BIN
dalek-logo.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 107 KiB

1
dalek-logo.svg Normal file
View file

File diff suppressed because one or more lines are too long
Side by side

After

Width:  |  Height:  |  Size: 59 KiB

112
src/backend/avx2/constants.rs Normal file
View file

@ -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 <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca>
//! 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),
])),
];

1057
src/backend/avx2/edwards.rs Normal file
View file

File diff suppressed because it is too large Load diff

689
src/backend/avx2/field.rs Normal file
View file

@ -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 <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca>
//! 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 );
}
}

487
src/backend/avx2/mod.rs Normal file
View file

@ -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 <isis@patternsinthevoid.net>
// - Henry de Valence <hdevalence@hdevalence.ca>
//! 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;

4
src/backend/mod.rs
View file

@ -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;

305
src/edwards.rs
View file

@ -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::<ProjectiveNielsPoint>::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::<ProjectiveNielsPoint>::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<Item = &'a Scalar>,
J: IntoIterator<Item = &'b ExtendedPoint>
{
//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::<ProjectiveNielsPoint>::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::<ProjectiveNielsPoint>::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<Item = &'a Scalar>,
J: IntoIterator<Item = &'b ExtendedPoint>
{
//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

6
src/lib.rs
View file

@ -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.

125
src/montgomery.rs
View file

@ -31,6 +31,7 @@
use core::ops::{Mul, MulAssign};
use constants;
use constants::APLUS2_OVER_FOUR;
use field::FieldElement;
use edwards::{ExtendedPoint, CompressedEdwardsY};
use scalar::Scalar;
@ -309,85 +310,52 @@ impl MontgomeryPoint {
CompressedMontgomeryU(u_affine.to_bytes())
}
}
/// Differential addition for single-coordinate Montgomery points.
///
/// Montgomery coordinates in projective 𝗣¹ space are odd in that 𝗣¹
/// inherits none of the group structure from E_(A,B). Hence, the mapping
/// of the group operation, `⊕`, is undefined for the pair `(x(P), x(Q))`;
/// that is, given `x(P)` and `x(Q)`, we cannot derive `x(P ⊕ Q)`. This is
/// due to the fact that, in Montgomery coordinates, `x(P)` determines `P`
/// only up to a sign, and thus we cannot differentiate `x(P ⊕ Q)` from
/// `x(P ⊖ Q)`. However, via differential addition, any three of the values
/// `{x(P), x(Q), x(P ⊕ Q), x(P ⊖ Q)}` determines the forth, so we can
/// define *pseudo-addition* for a singular coordinate.
///
/// # Warning
///
/// If the `difference` is the identity point, or a two torsion point, the
/// results of this method are not correct, but instead result in `(0:0)`
/// (an invalid projective point in the Montgomery model).
///
/// The doubling case is degenerate, in that `P ⦵ Q ∉ {O,T}`, where `T` is
/// the two torsion point.
fn differential_add(&self, that: &MontgomeryPoint,
difference: &MontgomeryPoint) -> MontgomeryPoint {
// XXX Do we want these debug assertions? We would need to implement
// XXX is_two_torsion_point(). —isis
// debug_assert!(!difference.is_identity()); // P ⦵ Q ∉ {O,T}
// debug_assert!(!difference.is_two_torsion_point());
/// DOCDOC
fn differential_add_and_double(P: &mut MontgomeryPoint, Q: &mut MontgomeryPoint,
difference: &MontgomeryPoint) {
let t0 = &P.U + &P.W;
let t1 = &P.U - &P.W;
let t2 = &Q.U + &Q.W;
let t3 = &Q.U - &Q.W;
let v1: FieldElement = &(&self.U + &self.W) * &(&that.U - &that.W);
let v2: FieldElement = &(&self.U - &self.W) * &(&that.U + &that.W);
let t4 = t0.square(); // (U_P + W_P)^2 = U_P^2 + 2 U_P W_P + W_P^2
let t5 = t1.square(); // (U_P - W_P)^2 = U_P^2 - 2 U_P W_P + W_P^2
MontgomeryPoint {
U: &difference.W * &(&v1 + &v2).square(), // does reduction on square()
W: &difference.U * &(&v1 - &v2).square(), // does reduction on square()
}
}
let t6 = &t4 - &t5; // 4 U_P W_P
/// Pseudo-doubling for single-coordinate Montgomery points.
///
/// Given a Montgomery U-coordinate of a point `P`, compute the
/// U-coordinate given by
///
///     differential_double: x(P) ⟼ x([2]P)
///
/// # Returns
///
/// A Montgomery point equal to doubling this one.
///
// XXX It seems possible that combining the differential_add() and
// XXX differential_double() methods would save a non-trivial amount of
// XXX computation in the ladder. —isis
fn differential_double(&self) -> MontgomeryPoint {
let mut v1: FieldElement;
let v2: FieldElement;
let v3: FieldElement;
let t7 = &t0 * &t3; // (U_P + W_P) (U_Q - W_Q) = U_P U_Q + W_P U_Q - U_P W_Q - W_P W_Q
let t8 = &t1 * &t2; // (U_P - W_P) (U_Q + W_Q) = U_P U_Q - W_P U_Q + U_P W_Q - W_P W_Q
v1 = (&self.U + &self.W).square();
v2 = (&self.U - &self.W).square();
let t9 = &t7 + &t8; // 2 (U_P U_Q - W_P W_Q)
let t10 = &t7 - &t8; // 2 (W_P U_Q - U_P W_Q)
let U: FieldElement = &v1 * &v2;
let t11 = t9.square(); // 4 (U_P U_Q - W_P W_Q)^2
let t12 = t10.square(); // 4 (W_P U_Q - U_P W_Q)^2
v1 -= &v2;
v3 = &(&constants::APLUS2_OVER_FOUR * &v1) + &v2;
let t13 = &APLUS2_OVER_FOUR * &t6; // (A + 2) U_P U_Q
let W: FieldElement = &v1 * &v3;
let t14 = &t4 * &t5; // ((U_P + W_P)(U_P - W_P))^2 = (U_P^2 - W_P^2)^2
let t15 = &t13 + &t5; // (U_P - W_P)^2 + (A + 2) U_P W_P
MontgomeryPoint{ U: U, W: W }
}
let t16 = &t6 * &t15; // 4 (U_P W_P) ((U_P - W_P)^2 + (A + 2) U_P W_P)
let t17 = &difference.U * &t12; // U_D * 4 (W_P U_Q - U_P W_Q)^2
let t18 = &difference.W * &t11; // W_D * 4 (U_P U_Q - W_P W_Q)^2
P.U = t14; // U_{P'} = (U_P + W_P)^2 (U_P - W_P)^2
P.W = t16; // W_{P'} = (4 U_P W_P) ((U_P - W_P)^2 + ((A + 2)/4) 4 U_P W_P)
Q.U = t18; // U_{Q'} = D_W * 4 (U_P U_Q - W_P W_Q)^2
Q.W = t17; // W_{Q'} = U_D * 4 (W_P U_Q - U_P W_Q)^2
}
/// Multiply this `MontgomeryPoint` by a `Scalar`.
///
/// The reader is refered to §5.3 of ["Montgomery Curves and Their Arithmetic"
/// by Craig Costello and Benjamin Smith](https://eprint.iacr.org/2017/212.pdf)
/// for an overview of side-channel-free Montgomery laddering algorithms.
impl<'a, 'b> Mul<&'b Scalar> for &'a MontgomeryPoint {
type Output = MontgomeryPoint;
fn mul(self, scalar: &'b Scalar) -> MontgomeryPoint {
// Algorithm 8 of Costello-Smith 2017
let mut x0: MontgomeryPoint = MontgomeryPoint::identity();
let mut x1: MontgomeryPoint = *self;
@ -399,8 +367,7 @@ impl<'a, 'b> Mul<&'b Scalar> for &'a MontgomeryPoint {
debug_assert!(mask == 0 || mask == 1);
x0.conditional_swap(&mut x1, mask);
x1 = x0.differential_add(&x1, &self);
x0 = x0.differential_double();
differential_add_and_double(&mut x0, &mut x1, &self);
}
x0.conditional_swap(&mut x1, bits[0] as u8);
x0
@ -476,14 +443,6 @@ mod test {
}
#[test]
fn differential_double_matches_double() {
let p: ExtendedPoint = constants::ED25519_BASEPOINT_POINT.double();
let q: MontgomeryPoint = BASE_COMPRESSED_MONTGOMERY.decompress().differential_double();
assert_eq!(p.to_montgomery().compress(), q.compress());
}
#[test]
#[cfg(feature="precomputed_tables")]
fn montgomery_ct_eq_ne() {
@ -506,26 +465,6 @@ mod test {
assert_eq!(p1.ct_eq(&p1), 1);
}
#[test]
#[cfg(feature="precomputed_tables")]
fn differential_add_matches_edwards_model() {
let mut csprng: OsRng = OsRng::new().unwrap();
let s1: Scalar = Scalar::random(&mut csprng);
let s2: Scalar = Scalar::random(&mut csprng);
let p1: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s1;
let p2: ExtendedPoint = &constants::ED25519_BASEPOINT_TABLE * &s2;
let diff: ExtendedPoint = &p1 - &p2;
let p1m: MontgomeryPoint = p1.to_montgomery();
let p2m: MontgomeryPoint = p2.to_montgomery();
let diffm: MontgomeryPoint = diff.to_montgomery();
let result = p1m.differential_add(&p2m, &diffm);
assert_eq!(result.compress(), (&p1 + &p2).to_montgomery().compress());
}
#[test]
#[cfg(feature="precomputed_tables")]
fn ladder_matches_scalarmult() {