1use core::arch::x86_64::*;
  2
  3// The following AVX2 implementation has been copied from plonky2:
  4// https://github.com/0xPolygonZero/plonky2/blob/main/plonky2/src/hash/arch/x86_64/poseidon_goldilocks_avx2_bmi2.rs
  5
  6// Preliminary notes:
  7// 1. AVX does not support addition with carry but 128-bit (2-word) addition can be easily
  8//    emulated. The method recognizes that for a + b overflowed iff (a + b) < a:
  9//        i. res_lo = a_lo + b_lo
 10//       ii. carry_mask = res_lo < a_lo
 11//      iii. res_hi = a_hi + b_hi - carry_mask
 12//    Notice that carry_mask is subtracted, not added. This is because AVX comparison instructions
 13//    return -1 (all bits 1) for true and 0 for false.
 14//
 15// 2. AVX does not have unsigned 64-bit comparisons. Those can be emulated with signed comparisons
 16//    by recognizing that a <u b iff a + (1 << 63) <s b + (1 << 63), where the addition wraps around
 17//    and the comparisons are unsigned and signed respectively. The shift function adds/subtracts
 18//    1 << 63 to enable this trick.
 19//      Example: addition with carry.
 20//        i. a_lo_s = shift(a_lo)
 21//       ii. res_lo_s = a_lo_s + b_lo
 22//      iii. carry_mask = res_lo_s <s a_lo_s
 23//       iv. res_lo = shift(res_lo_s)
 24//        v. res_hi = a_hi + b_hi - carry_mask
 25//    The suffix _s denotes a value that has been shifted by 1 << 63. The result of addition is
 26//    shifted if exactly one of the operands is shifted, as is the case on line ii. Line iii.
 27//    performs a signed comparison res_lo_s <s a_lo_s on shifted values to emulate unsigned
 28//    comparison res_lo <u a_lo on unshifted values. Finally, line iv. reverses the shift so the
 29//    result can be returned.
 30//      When performing a chain of calculations, we can often save instructions by letting the shift
 31//    propagate through and only undoing it when necessary. For example, to compute the addition of
 32//    three two-word (128-bit) numbers we can do:
 33//        i. a_lo_s = shift(a_lo)
 34//       ii. tmp_lo_s = a_lo_s + b_lo
 35//      iii. tmp_carry_mask = tmp_lo_s <s a_lo_s
 36//       iv. tmp_hi = a_hi + b_hi - tmp_carry_mask
 37//        v. res_lo_s = tmp_lo_s + c_lo
 38//       vi. res_carry_mask = res_lo_s <s tmp_lo_s
 39//      vii. res_lo = shift(res_lo_s)
 40//     viii. res_hi = tmp_hi + c_hi - res_carry_mask
 41//    Notice that the above 3-value addition still only requires two calls to shift, just like our
 42//    2-value addition.
 43
 44#[inline(always)]
 45pub fn branch_hint() {
 46    // NOTE: These are the currently supported assembly architectures. See the
 47    // [nightly reference](https://doc.rust-lang.org/nightly/reference/inline-assembly.html) for
 48    // the most up-to-date list.
 49    #[cfg(any(
 50        target_arch = "aarch64",
 51        target_arch = "arm",
 52        target_arch = "riscv32",
 53        target_arch = "riscv64",
 54        target_arch = "x86",
 55        target_arch = "x86_64",
 56    ))]
 57    unsafe {
 58        core::arch::asm!("", options(nomem, nostack, preserves_flags));
 59    }
 60}
 61
 62macro_rules! map3 {
 63    ($f:ident::<$l:literal>, $v:ident) => {
 64        ($f::<$l>($v.0), $f::<$l>($v.1), $f::<$l>($v.2))
 65    };
 66    ($f:ident::<$l:literal>, $v1:ident, $v2:ident) => {
 67        ($f::<$l>($v1.0, $v2.0), $f::<$l>($v1.1, $v2.1), $f::<$l>($v1.2, $v2.2))
 68    };
 69    ($f:ident, $v:ident) => {
 70        ($f($v.0), $f($v.1), $f($v.2))
 71    };
 72    ($f:ident, $v0:ident, $v1:ident) => {
 73        ($f($v0.0, $v1.0), $f($v0.1, $v1.1), $f($v0.2, $v1.2))
 74    };
 75    ($f:ident, rep $v0:ident, $v1:ident) => {
 76        ($f($v0, $v1.0), $f($v0, $v1.1), $f($v0, $v1.2))
 77    };
 78
 79    ($f:ident, $v0:ident, rep $v1:ident) => {
 80        ($f($v0.0, $v1), $f($v0.1, $v1), $f($v0.2, $v1))
 81    };
 82}
 83
 84#[inline(always)]
 85unsafe fn square3(
 86    x: (__m256i, __m256i, __m256i),
 87) -> ((__m256i, __m256i, __m256i), (__m256i, __m256i, __m256i)) {
 88    let x_hi = {
 89        // Move high bits to low position. The high bits of x_hi are ignored. Swizzle is faster than
 90        // bitshift. This instruction only has a floating-point flavor, so we cast to/from float.
 91        // This is safe and free.
 92        let x_ps = map3!(_mm256_castsi256_ps, x);
 93        let x_hi_ps = map3!(_mm256_movehdup_ps, x_ps);
 94        map3!(_mm256_castps_si256, x_hi_ps)
 95    };
 96
 97    // All pairwise multiplications.
 98    let mul_ll = map3!(_mm256_mul_epu32, x, x);
 99    let mul_lh = map3!(_mm256_mul_epu32, x, x_hi);
100    let mul_hh = map3!(_mm256_mul_epu32, x_hi, x_hi);
101
102    // Bignum addition, but mul_lh is shifted by 33 bits (not 32).
103    let mul_ll_hi = map3!(_mm256_srli_epi64::<33>, mul_ll);
104    let t0 = map3!(_mm256_add_epi64, mul_lh, mul_ll_hi);
105    let t0_hi = map3!(_mm256_srli_epi64::<31>, t0);
106    let res_hi = map3!(_mm256_add_epi64, mul_hh, t0_hi);
107
108    // Form low result by adding the mul_ll and the low 31 bits of mul_lh (shifted to the high
109    // position).
110    let mul_lh_lo = map3!(_mm256_slli_epi64::<33>, mul_lh);
111    let res_lo = map3!(_mm256_add_epi64, mul_ll, mul_lh_lo);
112
113    (res_lo, res_hi)
114}
115
116#[inline(always)]
117unsafe fn mul3(
118    x: (__m256i, __m256i, __m256i),
119    y: (__m256i, __m256i, __m256i),
120) -> ((__m256i, __m256i, __m256i), (__m256i, __m256i, __m256i)) {
121    let epsilon = _mm256_set1_epi64x(0xffffffff);
122    let x_hi = {
123        // Move high bits to low position. The high bits of x_hi are ignored. Swizzle is faster than
124        // bitshift. This instruction only has a floating-point flavor, so we cast to/from float.
125        // This is safe and free.
126        let x_ps = map3!(_mm256_castsi256_ps, x);
127        let x_hi_ps = map3!(_mm256_movehdup_ps, x_ps);
128        map3!(_mm256_castps_si256, x_hi_ps)
129    };
130    let y_hi = {
131        let y_ps = map3!(_mm256_castsi256_ps, y);
132        let y_hi_ps = map3!(_mm256_movehdup_ps, y_ps);
133        map3!(_mm256_castps_si256, y_hi_ps)
134    };
135
136    // All four pairwise multiplications
137    let mul_ll = map3!(_mm256_mul_epu32, x, y);
138    let mul_lh = map3!(_mm256_mul_epu32, x, y_hi);
139    let mul_hl = map3!(_mm256_mul_epu32, x_hi, y);
140    let mul_hh = map3!(_mm256_mul_epu32, x_hi, y_hi);
141
142    // Bignum addition
143    // Extract high 32 bits of mul_ll and add to mul_hl. This cannot overflow.
144    let mul_ll_hi = map3!(_mm256_srli_epi64::<32>, mul_ll);
145    let t0 = map3!(_mm256_add_epi64, mul_hl, mul_ll_hi);
146    // Extract low 32 bits of t0 and add to mul_lh. Again, this cannot overflow.
147    // Also, extract high 32 bits of t0 and add to mul_hh.
148    let t0_lo = map3!(_mm256_and_si256, t0, rep epsilon);
149    let t0_hi = map3!(_mm256_srli_epi64::<32>, t0);
150    let t1 = map3!(_mm256_add_epi64, mul_lh, t0_lo);
151    let t2 = map3!(_mm256_add_epi64, mul_hh, t0_hi);
152    // Lastly, extract the high 32 bits of t1 and add to t2.
153    let t1_hi = map3!(_mm256_srli_epi64::<32>, t1);
154    let res_hi = map3!(_mm256_add_epi64, t2, t1_hi);
155
156    // Form res_lo by combining the low half of mul_ll with the low half of t1 (shifted into high
157    // position).
158    let t1_lo = {
159        let t1_ps = map3!(_mm256_castsi256_ps, t1);
160        let t1_lo_ps = map3!(_mm256_moveldup_ps, t1_ps);
161        map3!(_mm256_castps_si256, t1_lo_ps)
162    };
163    let res_lo = map3!(_mm256_blend_epi32::<0xaa>, mul_ll, t1_lo);
164
165    (res_lo, res_hi)
166}
167
168/// Addition, where the second operand is `0 <= y < 0xffffffff00000001`.
169#[inline(always)]
170unsafe fn add_small(
171    x_s: (__m256i, __m256i, __m256i),
172    y: (__m256i, __m256i, __m256i),
173) -> (__m256i, __m256i, __m256i) {
174    let res_wrapped_s = map3!(_mm256_add_epi64, x_s, y);
175    let mask = map3!(_mm256_cmpgt_epi32, x_s, res_wrapped_s);
176    let wrapback_amt = map3!(_mm256_srli_epi64::<32>, mask); // EPSILON if overflowed else 0.
177    let res_s = map3!(_mm256_add_epi64, res_wrapped_s, wrapback_amt);
178    res_s
179}
180
181#[inline(always)]
182unsafe fn maybe_adj_sub(res_wrapped_s: __m256i, mask: __m256i) -> __m256i {
183    // The subtraction is very unlikely to overflow so we're best off branching.
184    // The even u32s in `mask` are meaningless, so we want to ignore them. `_mm256_testz_pd`
185    // branches depending on the sign bit of double-precision (64-bit) floats. Bit cast `mask` to
186    // floating-point (this is free).
187    let mask_pd = _mm256_castsi256_pd(mask);
188    // `_mm256_testz_pd(mask_pd, mask_pd) == 1` iff all sign bits are 0, meaning that underflow
189    // did not occur for any of the vector elements.
190    if _mm256_testz_pd(mask_pd, mask_pd) == 1 {
191        res_wrapped_s
192    } else {
193        branch_hint();
194        // Highly unlikely: underflow did occur. Find adjustment per element and apply it.
195        let adj_amount = _mm256_srli_epi64::<32>(mask); // EPSILON if underflow.
196        _mm256_sub_epi64(res_wrapped_s, adj_amount)
197    }
198}
199
200/// Addition, where the second operand is much smaller than `0xffffffff00000001`.
201#[inline(always)]
202unsafe fn sub_tiny(
203    x_s: (__m256i, __m256i, __m256i),
204    y: (__m256i, __m256i, __m256i),
205) -> (__m256i, __m256i, __m256i) {
206    let res_wrapped_s = map3!(_mm256_sub_epi64, x_s, y);
207    let mask = map3!(_mm256_cmpgt_epi32, res_wrapped_s, x_s);
208    let res_s = map3!(maybe_adj_sub, res_wrapped_s, mask);
209    res_s
210}
211
212#[inline(always)]
213unsafe fn reduce3(
214    (lo0, hi0): ((__m256i, __m256i, __m256i), (__m256i, __m256i, __m256i)),
215) -> (__m256i, __m256i, __m256i) {
216    let sign_bit = _mm256_set1_epi64x(i64::MIN);
217    let epsilon = _mm256_set1_epi64x(0xffffffff);
218    let lo0_s = map3!(_mm256_xor_si256, lo0, rep sign_bit);
219    let hi_hi0 = map3!(_mm256_srli_epi64::<32>, hi0);
220    let lo1_s = sub_tiny(lo0_s, hi_hi0);
221    let t1 = map3!(_mm256_mul_epu32, hi0, rep epsilon);
222    let lo2_s = add_small(lo1_s, t1);
223    let lo2 = map3!(_mm256_xor_si256, lo2_s, rep sign_bit);
224    lo2
225}
226
227#[inline(always)]
228unsafe fn mul_reduce(
229    a: (__m256i, __m256i, __m256i),
230    b: (__m256i, __m256i, __m256i),
231) -> (__m256i, __m256i, __m256i) {
232    reduce3(mul3(a, b))
233}
234
235#[inline(always)]
236unsafe fn square_reduce(state: (__m256i, __m256i, __m256i)) -> (__m256i, __m256i, __m256i) {
237    reduce3(square3(state))
238}
239
240#[inline(always)]
241unsafe fn exp_acc(
242    high: (__m256i, __m256i, __m256i),
243    low: (__m256i, __m256i, __m256i),
244    exp: usize,
245) -> (__m256i, __m256i, __m256i) {
246    let mut result = high;
247    for _ in 0..exp {
248        result = square_reduce(result);
249    }
250    mul_reduce(result, low)
251}
252
253#[inline(always)]
254unsafe fn do_apply_sbox(state: (__m256i, __m256i, __m256i)) -> (__m256i, __m256i, __m256i) {
255    let state2 = square_reduce(state);
256    let state4_unreduced = square3(state2);
257    let state3_unreduced = mul3(state2, state);
258    let state4 = reduce3(state4_unreduced);
259    let state3 = reduce3(state3_unreduced);
260    let state7_unreduced = mul3(state3, state4);
261    let state7 = reduce3(state7_unreduced);
262    state7
263}
264
265#[inline(always)]
266unsafe fn do_apply_inv_sbox(state: (__m256i, __m256i, __m256i)) -> (__m256i, __m256i, __m256i) {
267    // compute base^10540996611094048183 using 72 multiplications per array element
268    // 10540996611094048183 = b1001001001001001001001001001000110110110110110110110110110110111
269
270    // compute base^10
271    let t1 = square_reduce(state);
272
273    // compute base^100
274    let t2 = square_reduce(t1);
275
276    // compute base^100100
277    let t3 = exp_acc(t2, t2, 3);
278
279    // compute base^100100100100
280    let t4 = exp_acc(t3, t3, 6);
281
282    // compute base^100100100100100100100100
283    let t5 = exp_acc(t4, t4, 12);
284
285    // compute base^100100100100100100100100100100
286    let t6 = exp_acc(t5, t3, 6);
287
288    // compute base^1001001001001001001001001001000100100100100100100100100100100
289    let t7 = exp_acc(t6, t6, 31);
290
291    // compute base^1001001001001001001001001001000110110110110110110110110110110111
292    let a = square_reduce(square_reduce(mul_reduce(square_reduce(t7), t6)));
293    let b = mul_reduce(t1, mul_reduce(t2, state));
294    mul_reduce(a, b)
295}
296
297#[inline(always)]
298unsafe fn avx2_load(state: &[u64; 12]) -> (__m256i, __m256i, __m256i) {
299    (
300        _mm256_loadu_si256((&state[0..4]).as_ptr().cast::<__m256i>()),
301        _mm256_loadu_si256((&state[4..8]).as_ptr().cast::<__m256i>()),
302        _mm256_loadu_si256((&state[8..12]).as_ptr().cast::<__m256i>()),
303    )
304}
305
306#[inline(always)]
307unsafe fn avx2_store(buf: &mut [u64; 12], state: (__m256i, __m256i, __m256i)) {
308    _mm256_storeu_si256((&mut buf[0..4]).as_mut_ptr().cast::<__m256i>(), state.0);
309    _mm256_storeu_si256((&mut buf[4..8]).as_mut_ptr().cast::<__m256i>(), state.1);
310    _mm256_storeu_si256((&mut buf[8..12]).as_mut_ptr().cast::<__m256i>(), state.2);
311}
312
313#[inline(always)]
314pub unsafe fn apply_sbox(buffer: &mut [u64; 12]) {
315    let mut state = avx2_load(&buffer);
316    state = do_apply_sbox(state);
317    avx2_store(buffer, state);
318}
319
320#[inline(always)]
321pub unsafe fn apply_inv_sbox(buffer: &mut [u64; 12]) {
322    let mut state = avx2_load(&buffer);
323    state = do_apply_inv_sbox(state);
324    avx2_store(buffer, state);
325}