Revert "(work-in-progress) Partially remove unsafe_target_feature"

This reverts commit c67e430cfd.
This commit is contained in:
Jan Bujak 2023-06-05 07:38:55 +00:00
parent c67e430cfd
commit 94247a79d1
No known key found for this signature in database
GPG key ID: 3B438F83D43341D4
3 changed files with 325 additions and 399 deletions

View file

@ -48,6 +48,8 @@ use crate::backend::vector::avx2::constants::{
P_TIMES_16_HI, P_TIMES_16_LO, P_TIMES_2_HI, P_TIMES_2_LO,
};
use unsafe_target_feature::unsafe_target_feature;
/// Unpack 32-bit lanes into 64-bit lanes:
/// ```ascii,no_run
/// (a0, b0, a1, b1, c0, d0, c1, d1)
@ -57,9 +59,9 @@ use crate::backend::vector::avx2::constants::{
/// (a0, 0, b0, 0, c0, 0, d0, 0)
/// (a1, 0, b1, 0, c1, 0, d1, 0)
/// ```
#[target_feature(enable = "avx2")]
#[inline]
unsafe fn unpack_pair(src: u32x8) -> (u32x8, u32x8) {
#[unsafe_target_feature("avx2")]
#[inline(always)]
fn unpack_pair(src: u32x8) -> (u32x8, u32x8) {
let a: u32x8;
let b: u32x8;
let zero = u32x8::splat(0);
@ -81,9 +83,9 @@ unsafe fn unpack_pair(src: u32x8) -> (u32x8, u32x8) {
/// ```ascii,no_run
/// (a0, b0, a1, b1, c0, d0, c1, d1)
/// ```
#[target_feature(enable = "avx2")]
#[inline]
unsafe fn repack_pair(x: u32x8, y: u32x8) -> u32x8 {
#[unsafe_target_feature("avx2")]
#[inline(always)]
fn repack_pair(x: u32x8, y: u32x8) -> u32x8 {
unsafe {
use core::arch::x86_64::_mm256_blend_epi32;
use core::arch::x86_64::_mm256_shuffle_epi32;
@ -153,62 +155,43 @@ pub struct FieldElement2625x4(pub(crate) [u32x8; 5]);
use subtle::Choice;
use subtle::ConditionallySelectable;
#[unsafe_target_feature("avx2")]
impl ConditionallySelectable for FieldElement2625x4 {
#[inline(always)]
fn conditional_select(
a: &FieldElement2625x4,
b: &FieldElement2625x4,
choice: Choice,
) -> FieldElement2625x4 {
#[target_feature(enable = "avx2")]
unsafe fn inner(
a: &FieldElement2625x4,
b: &FieldElement2625x4,
choice: Choice,
) -> FieldElement2625x4 {
let mask = (-(choice.unwrap_u8() as i32)) as u32;
let mask_vec = u32x8::splat(mask);
FieldElement2625x4([
a.0[0] ^ (mask_vec & (a.0[0] ^ b.0[0])),
a.0[1] ^ (mask_vec & (a.0[1] ^ b.0[1])),
a.0[2] ^ (mask_vec & (a.0[2] ^ b.0[2])),
a.0[3] ^ (mask_vec & (a.0[3] ^ b.0[3])),
a.0[4] ^ (mask_vec & (a.0[4] ^ b.0[4])),
])
}
unsafe { inner(a, b, choice) }
let mask = (-(choice.unwrap_u8() as i32)) as u32;
let mask_vec = u32x8::splat(mask);
FieldElement2625x4([
a.0[0] ^ (mask_vec & (a.0[0] ^ b.0[0])),
a.0[1] ^ (mask_vec & (a.0[1] ^ b.0[1])),
a.0[2] ^ (mask_vec & (a.0[2] ^ b.0[2])),
a.0[3] ^ (mask_vec & (a.0[3] ^ b.0[3])),
a.0[4] ^ (mask_vec & (a.0[4] ^ b.0[4])),
])
}
#[inline(always)]
fn conditional_assign(&mut self, other: &FieldElement2625x4, choice: Choice) {
#[target_feature(enable = "avx2")]
unsafe fn inner(
itself: &mut FieldElement2625x4,
other: &FieldElement2625x4,
choice: Choice,
) {
let mask = (-(choice.unwrap_u8() as i32)) as u32;
let mask_vec = u32x8::splat(mask);
itself.0[0] ^= mask_vec & (itself.0[0] ^ other.0[0]);
itself.0[1] ^= mask_vec & (itself.0[1] ^ other.0[1]);
itself.0[2] ^= mask_vec & (itself.0[2] ^ other.0[2]);
itself.0[3] ^= mask_vec & (itself.0[3] ^ other.0[3]);
itself.0[4] ^= mask_vec & (itself.0[4] ^ other.0[4]);
}
unsafe { inner(self, other, choice) }
let mask = (-(choice.unwrap_u8() as i32)) as u32;
let mask_vec = u32x8::splat(mask);
self.0[0] ^= mask_vec & (self.0[0] ^ other.0[0]);
self.0[1] ^= mask_vec & (self.0[1] ^ other.0[1]);
self.0[2] ^= mask_vec & (self.0[2] ^ other.0[2]);
self.0[3] ^= mask_vec & (self.0[3] ^ other.0[3]);
self.0[4] ^= mask_vec & (self.0[4] ^ other.0[4]);
}
}
#[unsafe_target_feature("avx2")]
impl FieldElement2625x4 {
pub const ZERO: FieldElement2625x4 = FieldElement2625x4([u32x8::splat_const::<0>(); 5]);
/// Split this vector into an array of four (serial) field
/// elements.
#[rustfmt::skip] // keep alignment of extracted lanes
#[target_feature(enable = "avx2")]
pub unsafe fn split(&self) -> [FieldElement51; 4] {
pub fn split(&self) -> [FieldElement51; 4] {
let mut out = [FieldElement51::ZERO; 4];
for i in 0..5 {
let a_2i = self.0[i].extract::<0>() as u64; //
@ -235,8 +218,7 @@ impl FieldElement2625x4 {
/// that when this function is inlined, LLVM is able to lower the
/// shuffle using an immediate.
#[inline]
#[target_feature(enable = "avx2")]
pub unsafe fn shuffle(&self, control: Shuffle) -> FieldElement2625x4 {
pub fn shuffle(&self, control: Shuffle) -> FieldElement2625x4 {
#[inline(always)]
fn shuffle_lanes(x: u32x8, control: Shuffle) -> u32x8 {
unsafe {
@ -276,8 +258,7 @@ impl FieldElement2625x4 {
/// that this function can be inlined and LLVM can lower it to a
/// blend instruction using an immediate.
#[inline]
#[target_feature(enable = "avx2")]
pub unsafe fn blend(&self, other: FieldElement2625x4, control: Lanes) -> FieldElement2625x4 {
pub fn blend(&self, other: FieldElement2625x4, control: Lanes) -> FieldElement2625x4 {
#[inline(always)]
fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 {
unsafe {
@ -341,8 +322,7 @@ impl FieldElement2625x4 {
}
/// Convenience wrapper around `new(x,x,x,x)`.
#[target_feature(enable = "avx2")]
pub unsafe fn splat(x: &FieldElement51) -> FieldElement2625x4 {
pub fn splat(x: &FieldElement51) -> FieldElement2625x4 {
FieldElement2625x4::new(x, x, x, x)
}
@ -352,8 +332,7 @@ impl FieldElement2625x4 {
///
/// The resulting `FieldElement2625x4` is bounded with \\( b < 0.0002 \\).
#[rustfmt::skip] // keep alignment of computed lanes
#[target_feature(enable = "avx2")]
pub unsafe fn new(
pub fn new(
x0: &FieldElement51,
x1: &FieldElement51,
x2: &FieldElement51,
@ -392,8 +371,7 @@ impl FieldElement2625x4 {
///
/// The coefficients of the result are bounded with \\( b < 1 \\).
#[inline]
#[target_feature(enable = "avx2")]
pub unsafe fn negate_lazy(&self) -> FieldElement2625x4 {
pub fn negate_lazy(&self) -> FieldElement2625x4 {
// The limbs of self are bounded with b < 0.999, while the
// smallest limb of 2*p is 67108845 > 2^{26+0.9999}, so
// underflows are not possible.
@ -416,8 +394,7 @@ impl FieldElement2625x4 {
///
/// The coefficients of the result are bounded with \\( b < 1.6 \\).
#[inline]
#[target_feature(enable = "avx2")]
pub unsafe fn diff_sum(&self) -> FieldElement2625x4 {
pub fn diff_sum(&self) -> FieldElement2625x4 {
// tmp1 = (B, A, D, C)
let tmp1 = self.shuffle(Shuffle::BADC);
// tmp2 = (-A, B, -C, D)
@ -432,8 +409,7 @@ impl FieldElement2625x4 {
///
/// The coefficients of the result are bounded with \\( b < 0.0002 \\).
#[inline]
#[target_feature(enable = "avx2")]
pub unsafe fn reduce(&self) -> FieldElement2625x4 {
pub fn reduce(&self) -> FieldElement2625x4 {
let shifts = u32x8::new(26, 26, 25, 25, 26, 26, 25, 25);
let masks = u32x8::new(
(1 << 26) - 1,
@ -542,8 +518,7 @@ impl FieldElement2625x4 {
/// The coefficients of the result are bounded with \\( b < 0.007 \\).
#[inline]
#[rustfmt::skip] // keep alignment of carry chain
#[target_feature(enable = "avx2")]
unsafe fn reduce64(mut z: [u64x4; 10]) -> FieldElement2625x4 {
fn reduce64(mut z: [u64x4; 10]) -> FieldElement2625x4 {
// 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);
@ -619,8 +594,7 @@ impl FieldElement2625x4 {
///
/// The coefficients of the result are bounded with \\( b < 0.007 \\).
#[rustfmt::skip] // keep alignment of z* calculations
#[target_feature(enable = "avx2")]
pub unsafe fn square_and_negate_D(&self) -> FieldElement2625x4 {
pub fn square_and_negate_D(&self) -> FieldElement2625x4 {
#[inline(always)]
fn m(x: u32x8, y: u32x8) -> u64x4 {
x.mul32(y)
@ -707,6 +681,7 @@ impl FieldElement2625x4 {
}
}
#[unsafe_target_feature("avx2")]
impl Neg for FieldElement2625x4 {
type Output = FieldElement2625x4;
@ -722,46 +697,36 @@ impl Neg for FieldElement2625x4 {
/// # Postconditions
///
/// The coefficients of the result are bounded with \\( b < 0.0002 \\).
#[inline(always)]
#[inline]
fn neg(self) -> FieldElement2625x4 {
#[inline]
#[target_feature(enable = "avx2")]
unsafe fn inner(itself: FieldElement2625x4) -> FieldElement2625x4 {
FieldElement2625x4([
P_TIMES_16_LO - itself.0[0],
P_TIMES_16_HI - itself.0[1],
P_TIMES_16_HI - itself.0[2],
P_TIMES_16_HI - itself.0[3],
P_TIMES_16_HI - itself.0[4],
])
.reduce()
}
unsafe { inner(self) }
FieldElement2625x4([
P_TIMES_16_LO - self.0[0],
P_TIMES_16_HI - self.0[1],
P_TIMES_16_HI - self.0[2],
P_TIMES_16_HI - self.0[3],
P_TIMES_16_HI - self.0[4],
])
.reduce()
}
}
#[unsafe_target_feature("avx2")]
impl Add<FieldElement2625x4> for FieldElement2625x4 {
type Output = FieldElement2625x4;
/// Add two `FieldElement2625x4`s, without performing a reduction.
#[inline(always)]
#[inline]
fn add(self, rhs: FieldElement2625x4) -> FieldElement2625x4 {
#[inline]
#[target_feature(enable = "avx2")]
unsafe fn inner(itself: FieldElement2625x4, rhs: FieldElement2625x4) -> FieldElement2625x4 {
FieldElement2625x4([
itself.0[0] + rhs.0[0],
itself.0[1] + rhs.0[1],
itself.0[2] + rhs.0[2],
itself.0[3] + rhs.0[3],
itself.0[4] + rhs.0[4],
])
}
unsafe { inner(self, rhs) }
FieldElement2625x4([
self.0[0] + rhs.0[0],
self.0[1] + rhs.0[1],
self.0[2] + rhs.0[2],
self.0[3] + rhs.0[3],
self.0[4] + rhs.0[4],
])
}
}
#[unsafe_target_feature("avx2")]
impl Mul<(u32, u32, u32, u32)> for FieldElement2625x4 {
type Output = FieldElement2625x4;
/// Perform a multiplication by a vector of small constants.
@ -769,40 +734,32 @@ impl Mul<(u32, u32, u32, u32)> for FieldElement2625x4 {
/// # Postconditions
///
/// The coefficients of the result are bounded with \\( b < 0.007 \\).
#[inline(always)]
#[inline]
fn mul(self, scalars: (u32, u32, u32, u32)) -> FieldElement2625x4 {
#[inline]
#[target_feature(enable = "avx2")]
unsafe fn inner(
itself: FieldElement2625x4,
scalars: (u32, u32, u32, u32),
) -> FieldElement2625x4 {
let consts = u32x8::new(scalars.0, 0, scalars.1, 0, scalars.2, 0, scalars.3, 0);
let consts = u32x8::new(scalars.0, 0, scalars.1, 0, scalars.2, 0, scalars.3, 0);
let (b0, b1) = unpack_pair(itself.0[0]);
let (b2, b3) = unpack_pair(itself.0[1]);
let (b4, b5) = unpack_pair(itself.0[2]);
let (b6, b7) = unpack_pair(itself.0[3]);
let (b8, b9) = unpack_pair(itself.0[4]);
let (b0, b1) = unpack_pair(self.0[0]);
let (b2, b3) = unpack_pair(self.0[1]);
let (b4, b5) = unpack_pair(self.0[2]);
let (b6, b7) = unpack_pair(self.0[3]);
let (b8, b9) = unpack_pair(self.0[4]);
FieldElement2625x4::reduce64([
b0.mul32(consts),
b1.mul32(consts),
b2.mul32(consts),
b3.mul32(consts),
b4.mul32(consts),
b5.mul32(consts),
b6.mul32(consts),
b7.mul32(consts),
b8.mul32(consts),
b9.mul32(consts),
])
}
unsafe { inner(self, scalars) }
FieldElement2625x4::reduce64([
b0.mul32(consts),
b1.mul32(consts),
b2.mul32(consts),
b3.mul32(consts),
b4.mul32(consts),
b5.mul32(consts),
b6.mul32(consts),
b7.mul32(consts),
b8.mul32(consts),
b9.mul32(consts),
])
}
}
#[unsafe_target_feature("avx2")]
impl<'a, 'b> Mul<&'b FieldElement2625x4> for &'a FieldElement2625x4 {
type Output = FieldElement2625x4;
/// Multiply `self` by `rhs`.
@ -818,106 +775,98 @@ impl<'a, 'b> Mul<&'b FieldElement2625x4> for &'a FieldElement2625x4 {
/// The coefficients of the result are bounded with \\( b < 0.007 \\).
///
#[rustfmt::skip] // keep alignment of z* calculations
#[inline(always)]
#[inline]
fn mul(self, rhs: &'b FieldElement2625x4) -> FieldElement2625x4 {
#[inline]
#[target_feature(enable = "avx2")]
unsafe fn inner<'a, 'b>(itself: &'a FieldElement2625x4, rhs: &'b FieldElement2625x4) -> FieldElement2625x4 {
#[inline(always)]
fn m(x: u32x8, y: u32x8) -> u64x4 {
x.mul32(y)
}
#[inline(always)]
fn m_lo(x: u32x8, y: u32x8) -> u32x8 {
x.mul32(y).into()
}
let (x0, x1) = unpack_pair(itself.0[0]);
let (x2, x3) = unpack_pair(itself.0[1]);
let (x4, x5) = unpack_pair(itself.0[2]);
let (x6, x7) = unpack_pair(itself.0[3]);
let (x8, x9) = unpack_pair(itself.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);
let y6_19 = m_lo(v19, y6);
let y7_19 = m_lo(v19, y7);
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);
// The bounds on z[i] are the same as in the serial 32-bit code
// and the comment below is copied from there:
// How big is the contribution to z[i+j] from x[i], y[j]?
//
// Using the bounds above, we get:
//
// i even, j even: x[i]*y[j] < 2^(26+b)*2^(26+b) = 2*2^(51+2*b)
// i odd, j even: x[i]*y[j] < 2^(25+b)*2^(26+b) = 1*2^(51+2*b)
// i even, j odd: x[i]*y[j] < 2^(26+b)*2^(25+b) = 1*2^(51+2*b)
// i odd, j odd: 2*x[i]*y[j] < 2*2^(25+b)*2^(25+b) = 1*2^(51+2*b)
//
// We perform inline reduction mod p by replacing 2^255 by 19
// (since 2^255 - 19 = 0 mod p). This adds a factor of 19, so
// we get the bounds (z0 is the biggest one, but calculated for
// posterity here in case finer estimation is needed later):
//
// z0 < ( 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 249*2^(51 + 2*b)
// z1 < ( 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 154*2^(51 + 2*b)
// z2 < ( 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 195*2^(51 + 2*b)
// z3 < ( 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 118*2^(51 + 2*b)
// z4 < ( 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 141*2^(51 + 2*b)
// z5 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 82*2^(51 + 2*b)
// z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 87*2^(51 + 2*b)
// z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 )*2^(51 + 2b) = 46*2^(51 + 2*b)
// z8 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 )*2^(51 + 2b) = 33*2^(51 + 2*b)
// z9 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 )*2^(51 + 2b) = 10*2^(51 + 2*b)
//
// So z[0] fits into a u64 if 51 + 2*b + lg(249) < 64
// if b < 2.5.
// In fact this bound is slightly sloppy, since it treats both
// inputs x and y as being bounded by the same parameter b,
// while they are in fact bounded by b_x and b_y, and we
// already require that b_y < 1.75 in order to fit the
// multiplications by 19 into a u32. The tighter bound on b_y
// means we could get a tighter bound on the outputs, or a
// looser bound on b_x.
FieldElement2625x4::reduce64([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9])
#[inline(always)]
fn m(x: u32x8, y: u32x8) -> u64x4 {
x.mul32(y)
}
unsafe {
inner(self, rhs)
#[inline(always)]
fn m_lo(x: u32x8, y: u32x8) -> u32x8 {
x.mul32(y).into()
}
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);
let y6_19 = m_lo(v19, y6);
let y7_19 = m_lo(v19, y7);
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);
// The bounds on z[i] are the same as in the serial 32-bit code
// and the comment below is copied from there:
// How big is the contribution to z[i+j] from x[i], y[j]?
//
// Using the bounds above, we get:
//
// i even, j even: x[i]*y[j] < 2^(26+b)*2^(26+b) = 2*2^(51+2*b)
// i odd, j even: x[i]*y[j] < 2^(25+b)*2^(26+b) = 1*2^(51+2*b)
// i even, j odd: x[i]*y[j] < 2^(26+b)*2^(25+b) = 1*2^(51+2*b)
// i odd, j odd: 2*x[i]*y[j] < 2*2^(25+b)*2^(25+b) = 1*2^(51+2*b)
//
// We perform inline reduction mod p by replacing 2^255 by 19
// (since 2^255 - 19 = 0 mod p). This adds a factor of 19, so
// we get the bounds (z0 is the biggest one, but calculated for
// posterity here in case finer estimation is needed later):
//
// z0 < ( 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 249*2^(51 + 2*b)
// z1 < ( 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 154*2^(51 + 2*b)
// z2 < ( 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 195*2^(51 + 2*b)
// z3 < ( 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 118*2^(51 + 2*b)
// z4 < ( 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 141*2^(51 + 2*b)
// z5 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 + 1*19 + 1*19 )*2^(51 + 2b) = 82*2^(51 + 2*b)
// z6 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 + 2*19 + 1*19 )*2^(51 + 2b) = 87*2^(51 + 2*b)
// z7 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1*19 + 1*19 )*2^(51 + 2b) = 46*2^(51 + 2*b)
// z8 < ( 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1 + 2 + 1*19 )*2^(51 + 2b) = 33*2^(51 + 2*b)
// z9 < ( 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 )*2^(51 + 2b) = 10*2^(51 + 2*b)
//
// So z[0] fits into a u64 if 51 + 2*b + lg(249) < 64
// if b < 2.5.
// In fact this bound is slightly sloppy, since it treats both
// inputs x and y as being bounded by the same parameter b,
// while they are in fact bounded by b_x and b_y, and we
// already require that b_y < 1.75 in order to fit the
// multiplications by 19 into a u32. The tighter bound on b_y
// means we could get a tighter bound on the outputs, or a
// looser bound on b_x.
FieldElement2625x4::reduce64([z0, z1, z2, z3, z4, z5, z6, z7, z8, z9])
}
}

View file

@ -252,9 +252,9 @@ impl u64x4 {
}
/// Constructs a new instance.
#[target_feature(enable = "avx2")]
#[unsafe_target_feature("avx2")]
#[inline]
pub unsafe fn new(x0: u64, x1: u64, x2: u64, x3: u64) -> u64x4 {
pub fn new(x0: u64, x1: u64, x2: u64, x3: u64) -> u64x4 {
unsafe {
// _mm256_set_epi64 sets the underlying vector in reverse order of the args
u64x4(core::arch::x86_64::_mm256_set_epi64x(
@ -264,9 +264,9 @@ impl u64x4 {
}
/// Constructs a new instance with all of the elements initialized to the given value.
#[target_feature(enable = "avx2")]
#[unsafe_target_feature("avx2")]
#[inline]
pub unsafe fn splat(x: u64) -> u64x4 {
pub fn splat(x: u64) -> u64x4 {
unsafe { u64x4(core::arch::x86_64::_mm256_set1_epi64x(x as i64)) }
}
}
@ -303,18 +303,9 @@ impl u32x8 {
/// Constructs a new instance.
#[allow(clippy::too_many_arguments)]
#[target_feature(enable = "avx2")]
#[unsafe_target_feature("avx2")]
#[inline]
pub unsafe fn new(
x0: u32,
x1: u32,
x2: u32,
x3: u32,
x4: u32,
x5: u32,
x6: u32,
x7: u32,
) -> u32x8 {
pub fn new(x0: u32, x1: u32, x2: u32, x3: u32, x4: u32, x5: u32, x6: u32, x7: u32) -> u32x8 {
unsafe {
// _mm256_set_epi32 sets the underlying vector in reverse order of the args
u32x8(core::arch::x86_64::_mm256_set_epi32(
@ -325,9 +316,9 @@ impl u32x8 {
}
/// Constructs a new instance with all of the elements initialized to the given value.
#[target_feature(enable = "avx2")]
#[unsafe_target_feature("avx2")]
#[inline]
pub unsafe fn splat(x: u32) -> u32x8 {
pub fn splat(x: u32) -> u32x8 {
unsafe { u32x8(core::arch::x86_64::_mm256_set1_epi32(x as i32)) }
}
}

View file

@ -9,183 +9,169 @@
#![allow(non_snake_case)]
macro_rules! implement {
($module:ident, $backend_module:ident, $features:expr) => {
pub mod $module {
use alloc::vec::Vec;
#[unsafe_target_feature::unsafe_target_feature_specialize(
conditional("avx2", feature = "simd_avx2"),
conditional("avx512ifma,avx512vl", all(feature = "simd_avx512", nightly))
)]
pub mod spec {
use core::borrow::Borrow;
use core::cmp::Ordering;
use alloc::vec::Vec;
use crate::backend::vector::$backend_module::{CachedPoint, ExtendedPoint};
use core::borrow::Borrow;
use core::cmp::Ordering;
use crate::edwards::EdwardsPoint;
#[for_target_feature("avx2")]
use crate::backend::vector::avx2::{CachedPoint, ExtendedPoint};
#[for_target_feature("avx512ifma")]
use crate::backend::vector::ifma::{CachedPoint, ExtendedPoint};
use crate::edwards::EdwardsPoint;
use crate::scalar::Scalar;
use crate::traits::{Identity, VartimeMultiscalarMul};
/// Implements a version of Pippenger's algorithm.
///
/// See the documentation in the serial `scalar_mul::pippenger` module for details.
pub struct Pippenger;
impl VartimeMultiscalarMul for Pippenger {
type Point = EdwardsPoint;
fn optional_multiscalar_mul<I, J>(scalars: I, points: J) -> Option<EdwardsPoint>
where
I: IntoIterator,
I::Item: Borrow<Scalar>,
J: IntoIterator<Item = Option<EdwardsPoint>>,
{
let mut scalars = scalars.into_iter();
let size = scalars.by_ref().size_hint().0;
let w = if size < 500 {
6
} else if size < 800 {
7
} else {
8
};
let max_digit: usize = 1 << w;
let digits_count: usize = Scalar::to_radix_2w_size_hint(w);
let buckets_count: usize = max_digit / 2; // digits are signed+centered hence 2^w/2, excluding 0-th bucket
// Collect optimized scalars and points in a buffer for repeated access
// (scanning the whole collection per each digit position).
let scalars = scalars.map(|s| s.borrow().as_radix_2w(w));
let points = points
.into_iter()
.map(|p| p.map(|P| CachedPoint::from(ExtendedPoint::from(P))));
let scalars_points = scalars
.zip(points)
.map(|(s, maybe_p)| maybe_p.map(|p| (s, p)))
.collect::<Option<Vec<_>>>()?;
// Prepare 2^w/2 buckets.
// buckets[i] corresponds to a multiplication factor (i+1).
let mut buckets: Vec<ExtendedPoint> = (0..buckets_count)
.map(|_| ExtendedPoint::identity())
.collect();
let mut columns = (0..digits_count).rev().map(|digit_index| {
// Clear the buckets when processing another digit.
for bucket in &mut buckets {
*bucket = ExtendedPoint::identity();
}
// Iterate over pairs of (point, scalar)
// and add/sub the point to the corresponding bucket.
// Note: if we add support for precomputed lookup tables,
// we'll be adding/subtractiong point premultiplied by `digits[i]` to buckets[0].
for (digits, pt) in scalars_points.iter() {
// Widen digit so that we don't run into edge cases when w=8.
let digit = digits[digit_index] as i16;
match digit.cmp(&0) {
Ordering::Greater => {
let b = (digit - 1) as usize;
buckets[b] = &buckets[b] + pt;
}
Ordering::Less => {
let b = (-digit - 1) as usize;
buckets[b] = &buckets[b] - pt;
}
Ordering::Equal => {}
}
}
// Add the buckets applying the multiplication factor to each bucket.
// The most efficient way to do that is to have a single sum with two running sums:
// an intermediate sum from last bucket to the first, and a sum of intermediate sums.
//
// For example, to add buckets 1*A, 2*B, 3*C we need to add these points:
// C
// C B
// C B A Sum = C + (C+B) + (C+B+A)
let mut buckets_intermediate_sum = buckets[buckets_count - 1];
let mut buckets_sum = buckets[buckets_count - 1];
for i in (0..(buckets_count - 1)).rev() {
buckets_intermediate_sum =
&buckets_intermediate_sum + &CachedPoint::from(buckets[i]);
buckets_sum = &buckets_sum + &CachedPoint::from(buckets_intermediate_sum);
}
buckets_sum
});
// Take the high column as an initial value to avoid wasting time doubling the identity element in `fold()`.
// `unwrap()` always succeeds because we know we have more than zero digits.
let hi_column = columns.next().unwrap();
Some(
columns
.fold(hi_column, |total, p| {
&total.mul_by_pow_2(w as u32) + &CachedPoint::from(p)
})
.into(),
)
}
}
#[cfg(test)]
mod test {
#[test]
fn test_vartime_pippenger() {
use super::*;
use crate::constants;
use crate::scalar::Scalar;
use crate::traits::{Identity, VartimeMultiscalarMul};
/// Implements a version of Pippenger's algorithm.
///
/// See the documentation in the serial `scalar_mul::pippenger` module for details.
pub struct Pippenger;
// Reuse points across different tests
let mut n = 512;
let x = Scalar::from(2128506u64).invert();
let y = Scalar::from(4443282u64).invert();
let points: Vec<_> = (0..n)
.map(|i| constants::ED25519_BASEPOINT_POINT * Scalar::from(1 + i as u64))
.collect();
let scalars: Vec<_> = (0..n)
.map(|i| x + (Scalar::from(i as u64) * y)) // fast way to make ~random but deterministic scalars
.collect();
impl VartimeMultiscalarMul for Pippenger {
type Point = EdwardsPoint;
let premultiplied: Vec<EdwardsPoint> = scalars
.iter()
.zip(points.iter())
.map(|(sc, pt)| sc * pt)
.collect();
#[inline(always)]
fn optional_multiscalar_mul<I, J>(scalars: I, points: J) -> Option<EdwardsPoint>
where
I: IntoIterator,
I::Item: Borrow<Scalar>,
J: IntoIterator<Item = Option<EdwardsPoint>>,
{
#[target_feature(enable = $features)]
unsafe fn inner<I, J>(scalars: I, points: J) -> Option<EdwardsPoint>
where
I: IntoIterator,
I::Item: Borrow<Scalar>,
J: IntoIterator<Item = Option<EdwardsPoint>>,
{
let mut scalars = scalars.into_iter();
let size = scalars.by_ref().size_hint().0;
let w = if size < 500 {
6
} else if size < 800 {
7
} else {
8
};
while n > 0 {
let scalars = &scalars[0..n].to_vec();
let points = &points[0..n].to_vec();
let control: EdwardsPoint = premultiplied[0..n].iter().sum();
let max_digit: usize = 1 << w;
let digits_count: usize = Scalar::to_radix_2w_size_hint(w);
let buckets_count: usize = max_digit / 2; // digits are signed+centered hence 2^w/2, excluding 0-th bucket
let subject = Pippenger::vartime_multiscalar_mul(scalars.clone(), points.clone());
// Collect optimized scalars and points in a buffer for repeated access
// (scanning the whole collection per each digit position).
let scalars = scalars.map(|s| s.borrow().as_radix_2w(w));
assert_eq!(subject.compress(), control.compress());
let points = points
.into_iter()
.map(|p| p.map(|P| CachedPoint::from(ExtendedPoint::from(P))));
let scalars_points = scalars
.zip(points)
.map(|(s, maybe_p)| maybe_p.map(|p| (s, p)))
.collect::<Option<Vec<_>>>()?;
// Prepare 2^w/2 buckets.
// buckets[i] corresponds to a multiplication factor (i+1).
let mut buckets: Vec<ExtendedPoint> = (0..buckets_count)
.map(|_| ExtendedPoint::identity())
.collect();
let mut columns = (0..digits_count).rev().map(|digit_index| {
// Clear the buckets when processing another digit.
for bucket in &mut buckets {
*bucket = ExtendedPoint::identity();
}
// Iterate over pairs of (point, scalar)
// and add/sub the point to the corresponding bucket.
// Note: if we add support for precomputed lookup tables,
// we'll be adding/subtractiong point premultiplied by `digits[i]` to buckets[0].
for (digits, pt) in scalars_points.iter() {
// Widen digit so that we don't run into edge cases when w=8.
let digit = digits[digit_index] as i16;
match digit.cmp(&0) {
Ordering::Greater => {
let b = (digit - 1) as usize;
buckets[b] = &buckets[b] + pt;
}
Ordering::Less => {
let b = (-digit - 1) as usize;
buckets[b] = &buckets[b] - pt;
}
Ordering::Equal => {}
}
}
// Add the buckets applying the multiplication factor to each bucket.
// The most efficient way to do that is to have a single sum with two running sums:
// an intermediate sum from last bucket to the first, and a sum of intermediate sums.
//
// For example, to add buckets 1*A, 2*B, 3*C we need to add these points:
// C
// C B
// C B A Sum = C + (C+B) + (C+B+A)
let mut buckets_intermediate_sum = buckets[buckets_count - 1];
let mut buckets_sum = buckets[buckets_count - 1];
for i in (0..(buckets_count - 1)).rev() {
buckets_intermediate_sum =
&buckets_intermediate_sum + &CachedPoint::from(buckets[i]);
buckets_sum =
&buckets_sum + &CachedPoint::from(buckets_intermediate_sum);
}
buckets_sum
});
// Take the high column as an initial value to avoid wasting time doubling the identity element in `fold()`.
// `unwrap()` always succeeds because we know we have more than zero digits.
let hi_column = columns.next().unwrap();
Some(
columns
.fold(hi_column, |total, p| {
&total.mul_by_pow_2(w as u32) + &CachedPoint::from(p)
})
.into(),
)
}
unsafe { inner(scalars, points) }
}
}
#[cfg(test)]
#[cfg(target_feature = $features)]
mod test {
#[test]
fn test_vartime_pippenger() {
use super::*;
use crate::constants;
use crate::scalar::Scalar;
// Reuse points across different tests
let mut n = 512;
let x = Scalar::from(2128506u64).invert();
let y = Scalar::from(4443282u64).invert();
let points: Vec<_> = (0..n)
.map(|i| constants::ED25519_BASEPOINT_POINT * Scalar::from(1 + i as u64))
.collect();
let scalars: Vec<_> = (0..n)
.map(|i| x + (Scalar::from(i as u64) * y)) // fast way to make ~random but deterministic scalars
.collect();
let premultiplied: Vec<EdwardsPoint> = scalars
.iter()
.zip(points.iter())
.map(|(sc, pt)| sc * pt)
.collect();
while n > 0 {
let scalars = &scalars[0..n].to_vec();
let points = &points[0..n].to_vec();
let control: EdwardsPoint = premultiplied[0..n].iter().sum();
let subject =
Pippenger::vartime_multiscalar_mul(scalars.clone(), points.clone());
assert_eq!(subject.compress(), control.compress());
n = n / 2;
}
}
n = n / 2;
}
}
};
}
}
#[cfg(feature = "simd_avx2")]
implement!(spec_avx2, avx2, "avx2");
#[cfg(all(feature = "simd_avx512", nightly))]
implement!(spec_avx512ifma_avx512vl, ifma, "avx512ifma,avx512vl");