Eliminate a carry pass through tighter bounds checks

This commit is contained in:
Henry de Valence 2017-12-08 17:49:35 -08:00
parent d62fc7caf1
commit 640888198e
2 changed files with 83 additions and 31 deletions

View file

@ -26,7 +26,7 @@ use scalar::Scalar;
use traits::Identity;
use backend::avx2::field::FieldElement32x4;
use backend::avx2::field::P_TIMES_2;
use backend::avx2::field::P_TIMES_2_MASKED;
use backend::avx2;
@ -119,9 +119,9 @@ impl ExtendedPoint {
t0.0[3] = _mm256_blend_epi32(t0.0[3].into(), P.0[3].into(), 0b01011111).into();
t0.0[4] = _mm256_blend_epi32(t0.0[4].into(), P.0[4].into(), 0b01011111).into();
t1 = t0.square();
t1 = t0.square(0b11_00_00_00);
// Now t1 = (S1 S2 S3 S4)
// Now t1 = (S1 S2 S3 -S4)
let c0 = u32x8::new(0,0,2,2,0,0,2,2); // (ABCD) -> (AAAA)
let c1 = u32x8::new(1,1,3,3,1,1,3,3); // (ABCD) -> (BBBB)
@ -134,9 +134,9 @@ impl ExtendedPoint {
// + | S2 | | | S2 |
// + | | | S3 | |
// + | | | S3 | |
// + | | 2p | 2p | 2p |
// + | | | |-S4 |
// + | | 2p | 2p | |
// - | | S2 | S2 | |
// - | | | | S4 |
// =======================
// S5 S6 S8 S9
//
@ -146,17 +146,17 @@ impl ExtendedPoint {
// + | 2^26 | | | 2^26 | + | 2^25 | | | 2^25 |
// + | | | 2^26 | | + | | | 2^25 | |
// + | | | 2^26 | | + | | | 2^25 | |
// + | | 2^27 | 2^27 | 2^27 | + | | 2^26 | 2^26 | 2^26 |
// + | | | | 2^26 | + | | | | 2^25 |
// + | | 2^27 | 2^27 | | + | | 2^26 | 2^26 | |
// - | | 0 | 0 | | - | | 0 | 0 | |
// - | | | | 0 | - | | | | 0 |
// =================================== ===================================
// < 2^27 2^27.59 2^28.33 2^28 2^26 2^26.59 2^27.33 2^27
// < 2^27 2^27.59 2^28.33 2^27.59 2^26 2^26.59 2^27.33 2^27.59
//
// So, the bit-excess for (S5 S6 S8 S9) is (1, 1.59, 2.33, 2).
// So, the bit-excess for (S5 S6 S8 S9) is (1, 1.59, 2.33, 1.59).
//
// However the multiplication routine only allows (1.75, 1.75, 1.75, 1.75).
//
// This is because we need to have 19*y[i] < 2^32. Otherwise I think we could get b < 2.5.
// This is because we need to have 19*y[i] < 2^32. Otherwise the bound would be b < 2.5.
//
// Can we tighten these bounds to avoid a reduction? Alternately, can we do better than
// the 64-bit reduction that reduce32() calls internally?
@ -169,17 +169,14 @@ impl ExtendedPoint {
let zero = i32x8::splat(0);
let S1 = _mm256_permutevar8x32_epi32(t1.0[i], c0);
let S2 = _mm256_permutevar8x32_epi32(t1.0[i], c1);
let S3_2 = _mm256_blend_epi32(zero, (t1.0[i] + t1.0[i]).into(), 0b01010000).into();
t0.0[i] = (P_TIMES_2.0[i] + S3_2) + S1;
let S3_2: u32x8 = _mm256_blend_epi32(zero, (t1.0[i] + t1.0[i]).into(), 0b01010000).into();
// tmp0 = (0 0 2*S3 -S4)
let tmp0: u32x8 = _mm256_blend_epi32(S3_2.into(), t1.0[i].into(), 0b10100000).into();
t0.0[i] = (P_TIMES_2_MASKED.0[i] + tmp0) + S1;
t0.0[i] = t0.0[i] + _mm256_blend_epi32(zero, S2.into(), 0b10100101).into();
let S4 = _mm256_blend_epi32(zero, t1.0[i].into(), 0b10100000);
let sub = _mm256_blend_epi32(S2.into(), S4, 0b10100101).into();
t0.0[i] = t0.0[i] - sub;
t0.0[i] = t0.0[i] - _mm256_blend_epi32(S2.into(), zero, 0b10100101).into();
}
// This is really sad, see above
t0.reduce32();
let c0 = u32x8::new(4,0,6,2,4,0,6,2); // (ABCD) -> (CACA)
let c1 = u32x8::new(5,1,7,3,1,5,3,7); // (ABCD) -> (DBBD)

View file

@ -27,6 +27,14 @@ pub(crate) static P_TIMES_2: FieldElement32x4 = FieldElement32x4([
u32x8::new(134217726, 134217726, 67108862, 67108862, 134217726, 134217726, 67108862, 67108862)
]);
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)
]);
/// A vector of four `FieldElements`, implemented using AVX2.
#[derive(Clone, Copy, Debug)]
pub(crate) struct FieldElement32x4(pub(crate) [u32x8; 5]);
@ -400,7 +408,13 @@ pub fn repack_pair(x: u32x8, y: u32x8) -> u32x8 {
}
impl FieldElement32x4 {
pub fn square(&self) -> 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;
@ -436,16 +450,55 @@ impl FieldElement32x4 {
let x8_19 = m_lo(v19, x8);
let x9_19 = m_lo(v19, x9);
let 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 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 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 z3 = m(x0_2,x3) + m(x1_2,x2) + m(x5_2,x8_19) + ((m(x4,x9_19) + m(x6,x7_19)) << 1);
let 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 z5 = m(x0_2,x5) + m(x1_2,x4) + m(x2_2,x3) + m(x7_2,x8_19) + ((m(x6,x9_19)) << 1);
let 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 z7 = m(x0_2,x7) + m(x1_2,x6) + m(x2_2,x5) + m(x3_2,x4) + ((m(x8,x9_19)) << 1);
let 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 z9 = m(x0_2,x9) + m(x1_2,x8) + m(x2_2,x7) + m(x3_2,x6) + m(x4_2,x5);
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])
}
@ -556,12 +609,14 @@ mod test {
let vec = FieldElement32x4::new(&x0, &x1, &x2, &x3);
let result = vec.square().split();
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);
assert_eq!(result[3], -&(&x3 * &x3));
}