1#include <immintrin.h>
  2
  3#include "reim_fft_private.h"
  4
  5void reim_from_znx64_bnd50_fma(const REIM_FROM_ZNX64_PRECOMP* precomp, void* r, const int64_t* x) {
  6  static const double EXPO = INT64_C(1) << 52;
  7  const int64_t ADD_CST = INT64_C(1) << 51;
  8  const double SUB_CST = INT64_C(3) << 51;
  9
 10  const __m256d SUB_CST_4 = _mm256_set1_pd(SUB_CST);
 11  const __m256i ADD_CST_4 = _mm256_set1_epi64x(ADD_CST);
 12  const __m256d EXPO_4 = _mm256_set1_pd(EXPO);
 13
 14  double(*out)[4] = (double(*)[4])r;
 15  __m256i* in = (__m256i*)x;
 16  __m256i* inend = (__m256i*)(x + (precomp->m << 1));
 17  do {
 18    // read the next value
 19    __m256i a = _mm256_loadu_si256(in);
 20    a = _mm256_add_epi64(a, ADD_CST_4);
 21    __m256d ad = _mm256_castsi256_pd(a);
 22    ad = _mm256_or_pd(ad, EXPO_4);
 23    ad = _mm256_sub_pd(ad, SUB_CST_4);
 24    // store the next value
 25    _mm256_storeu_pd(out[0], ad);
 26    ++out;
 27    ++in;
 28  } while (in < inend);
 29}
 30
 31// version where the output norm can be as big as 2^63
 32void reim_to_znx64_avx2_bnd63_fma(const REIM_TO_ZNX64_PRECOMP* precomp, int64_t* r, const void* x) {
 33  static const uint64_t SIGN_MASK = 0x8000000000000000UL;
 34  static const uint64_t EXPO_MASK = 0x7FF0000000000000UL;
 35  static const uint64_t MANTISSA_MASK = 0x000FFFFFFFFFFFFFUL;
 36  static const uint64_t MANTISSA_MSB = 0x0010000000000000UL;
 37  const double divisor_bits = precomp->divisor * ((double)(INT64_C(1) << 52));
 38  const double offset = precomp->divisor / 2.;
 39
 40  const __m256d SIGN_MASK_4 = _mm256_castsi256_pd(_mm256_set1_epi64x(SIGN_MASK));
 41  const __m256i EXPO_MASK_4 = _mm256_set1_epi64x(EXPO_MASK);
 42  const __m256i MANTISSA_MASK_4 = _mm256_set1_epi64x(MANTISSA_MASK);
 43  const __m256i MANTISSA_MSB_4 = _mm256_set1_epi64x(MANTISSA_MSB);
 44  const __m256d offset_4 = _mm256_set1_pd(offset);
 45  const __m256i divi_bits_4 = _mm256_castpd_si256(_mm256_set1_pd(divisor_bits));
 46
 47  double(*in)[4] = (double(*)[4])x;
 48  __m256i* out = (__m256i*)r;
 49  __m256i* outend = (__m256i*)(r + (precomp->m << 1));
 50  do {
 51    // read the next value
 52    __m256d a = _mm256_loadu_pd(in[0]);
 53    // a += sign(a) * m/2
 54    __m256d asign = _mm256_and_pd(a, SIGN_MASK_4);
 55    a = _mm256_add_pd(a, _mm256_or_pd(asign, offset_4));
 56    // sign: either 0 or -1
 57    __m256i sign_mask = _mm256_castpd_si256(asign);
 58    sign_mask = _mm256_sub_epi64(_mm256_set1_epi64x(0), _mm256_srli_epi64(sign_mask, 63));
 59    // compute the exponents
 60    __m256i a0exp = _mm256_and_si256(_mm256_castpd_si256(a), EXPO_MASK_4);
 61    __m256i a0lsh = _mm256_sub_epi64(a0exp, divi_bits_4);
 62    __m256i a0rsh = _mm256_sub_epi64(divi_bits_4, a0exp);
 63    a0lsh = _mm256_srli_epi64(a0lsh, 52);
 64    a0rsh = _mm256_srli_epi64(a0rsh, 52);
 65    // compute the new mantissa
 66    __m256i a0pos = _mm256_and_si256(_mm256_castpd_si256(a), MANTISSA_MASK_4);
 67    a0pos = _mm256_or_si256(a0pos, MANTISSA_MSB_4);
 68    a0lsh = _mm256_sllv_epi64(a0pos, a0lsh);
 69    a0rsh = _mm256_srlv_epi64(a0pos, a0rsh);
 70    __m256i final = _mm256_or_si256(a0lsh, a0rsh);
 71    // negate if the sign was negative
 72    final = _mm256_xor_si256(final, sign_mask);
 73    final = _mm256_sub_epi64(final, sign_mask);
 74    // read the next value
 75    _mm256_storeu_si256(out, final);
 76    ++out;
 77    ++in;
 78  } while (out < outend);
 79}
 80
 81// version where the output norm can be as big as 2^50
 82void reim_to_znx64_avx2_bnd50_fma(const REIM_TO_ZNX64_PRECOMP* precomp, int64_t* r, const void* x) {
 83  static const uint64_t MANTISSA_MASK = 0x000FFFFFFFFFFFFFUL;
 84  const int64_t SUB_CST = INT64_C(1) << 51;
 85  const double add_cst = precomp->divisor * ((double)(INT64_C(3) << 51));
 86
 87  const __m256i SUB_CST_4 = _mm256_set1_epi64x(SUB_CST);
 88  const __m256d add_cst_4 = _mm256_set1_pd(add_cst);
 89  const __m256i MANTISSA_MASK_4 = _mm256_set1_epi64x(MANTISSA_MASK);
 90
 91  double(*in)[4] = (double(*)[4])x;
 92  __m256i* out = (__m256i*)r;
 93  __m256i* outend = (__m256i*)(r + (precomp->m << 1));
 94  do {
 95    // read the next value
 96    __m256d a = _mm256_loadu_pd(in[0]);
 97    a = _mm256_add_pd(a, add_cst_4);
 98    __m256i ai = _mm256_castpd_si256(a);
 99    ai = _mm256_and_si256(ai, MANTISSA_MASK_4);
100    ai = _mm256_sub_epi64(ai, SUB_CST_4);
101    // store the next value
102    _mm256_storeu_si256(out, ai);
103    ++out;
104    ++in;
105  } while (out < outend);
106}