Как эффективно выполнять двойные /int64 преобразования с SSE/AVX?

SSE2 имеет инструкции для преобразования векторов между числами с плавающей точкой одинарной точности и 32-разрядными целыми числами.

  • _mm_cvtps_epi32()
  • _mm_cvtepi32_ps()

Но нет эквивалентов для двойных и 64-битных целых чисел. Другими словами, они отсутствуют:

  • _mm_cvtpd_epi64()
  • _mm_cvtepi64_pd()

Кажется, у AVX их тоже нет.

Какой самый эффективный способ симулировать эти свойства?

Источник

4 ответа

Решение

Если вы готовы срезать углы, double <-> int64 Преобразование может быть сделано только в двух инструкциях:

  • Если вы не заботитесь о бесконечности или NaN,
  • За double <-> int64_t, вы заботитесь только о значениях в диапазоне [-2^51, 2^51],
  • За double <-> uint64_t, вы заботитесь только о значениях в диапазоне [0, 2^52),

double -> uint64_t

//  Only works for inputs in the range: [0, 2^52)
__m128i double_to_uint64(__m128d x){
    x = _mm_add_pd(x, _mm_set1_pd(0x0010000000000000));
    return _mm_xor_si128(
        _mm_castpd_si128(x),
        _mm_castpd_si128(_mm_set1_pd(0x0010000000000000))
    );
}

double -> int64_t

//  Only works for inputs in the range: [-2^51, 2^51]
__m128i double_to_int64(__m128d x){
    x = _mm_add_pd(x, _mm_set1_pd(0x0018000000000000));
    return _mm_sub_epi64(
        _mm_castpd_si128(x),
        _mm_castpd_si128(_mm_set1_pd(0x0018000000000000))
    );
}

uint64_t -> двойной

//  Only works for inputs in the range: [0, 2^52)
__m128d uint64_to_double(__m128i x){
    x = _mm_or_si128(x, _mm_castpd_si128(_mm_set1_pd(0x0010000000000000)));
    return _mm_sub_pd(_mm_castsi128_pd(x), _mm_set1_pd(0x0010000000000000));
}

int64_t -> double

//  Only works for inputs in the range: [-2^51, 2^51]
__m128d int64_to_double(__m128i x){
    x = _mm_add_epi64(x, _mm_castpd_si128(_mm_set1_pd(0x0018000000000000)));
    return _mm_sub_pd(_mm_castsi128_pd(x), _mm_set1_pd(0x0018000000000000));
}

Поведение при округлении:

  • Для double -> uint64_t преобразование, округление работает правильно, следуя текущему режиму округления. (который обычно округляется до четного)
  • Для double -> int64_t преобразование, округление будет следовать текущему режиму округления для всех режимов, кроме усечения. Если текущим режимом округления является усечение (округление до нуля), оно фактически округляется до отрицательной бесконечности.

Как это работает?

Несмотря на то, что этот трюк состоит всего из двух инструкций, он не совсем понятен

Ключ должен признать, что для чисел с плавающей запятой двойной точности значения в диапазоне [2^52, 2^53) иметь "двоичное место" чуть ниже самого нижнего бита мантиссы. Другими словами, если вы обнуляете биты экспоненты и знака, мантисса становится именно целочисленным представлением.

Преобразовать x от double -> uint64_tдобавляешь магический номер M которая является значением с плавающей запятой 2^52, Это ставит x в "нормализованный" диапазон [2^52, 2^53) и удобно округляет биты дробной части.

Теперь осталось только удалить верхние 12 бит. Это легко сделать, маскируя его. Самый быстрый способ - признать, что эти старшие 12 бит идентичны M, Таким образом, вместо того, чтобы вводить дополнительную константу маски, мы можем просто вычесть или XOR M, XOR имеет большую пропускную способность.

Преобразование из uint64_t -> double просто обратная сторона этого процесса. Вы добавляете обратно экспонентные биты M, Затем ненормализовать число путем вычитания M в плавающей запятой.

Целочисленные преобразования со знаком немного сложнее, так как вам нужно иметь дело с расширением знака дополнения 2. Я оставлю это как упражнение для читателя.

Связанный: объяснен быстрый метод округления двойного к 32-битному целому

Источник

Этот ответ о 64-битном целочисленном значении для двойного преобразования без срезания углов. Мы предполагаем, что как целочисленный, так и двойной выход находятся в регистрах AVX шириной 256 бит. Рассматриваются два подхода:

  1. int64_to_double_based_on_cvtsi2sd(): как предлагается в комментариях к вопросу, используйте cvtsi2sd 4 раза вместе с перетасовкой данных. К сожалению оба cvtsi2sd и для команд перетасовки данных требуется порт 5 выполнения. Это ограничивает производительность этого подхода.

  2. int64_to_double_full_range(): мы можем использовать метод быстрого преобразования Mysticial дважды, чтобы получить точное преобразование для полного 64-битного целочисленного диапазона. 64-разрядное целое число разбивается на 32-битное и 32-битное старшие, как и в ответах на этот вопрос: Как выполнить преобразование uint32/float с SSE?, Каждый из этих кусков подходит для целочисленного двойного преобразования Mysticial. Наконец, верхняя часть умножается на 2^32 и добавляется к нижней части. Преобразование со знаком немного сложнее, чем преобразование без знака (uint64_to_double_full_range()), так как srai_epi64() не существует

Код:

#include <stdio.h>
#include <immintrin.h>
#include <stdint.h>

/* 
gcc -O3 -Wall -m64 -mfma -mavx2 -march=broadwell cvt_int_64_double.c
./a.out A
time ./a.out B
time ./a.out C
etc.
*/


inline __m256d uint64_to_double256(__m256i x){                  /*  Mysticial's fast uint64_to_double. Works for inputs in the range: [0, 2^52)     */
    x = _mm256_or_si256(x, _mm256_castpd_si256(_mm256_set1_pd(0x0010000000000000)));
    return _mm256_sub_pd(_mm256_castsi256_pd(x), _mm256_set1_pd(0x0010000000000000));
}

inline __m256d int64_to_double256(__m256i x){                   /*  Mysticial's fast int64_to_double. Works for inputs in the range: (-2^51, 2^51)  */
    x = _mm256_add_epi64(x, _mm256_castpd_si256(_mm256_set1_pd(0x0018000000000000)));
    return _mm256_sub_pd(_mm256_castsi256_pd(x), _mm256_set1_pd(0x0018000000000000));
}


__m256d int64_to_double_full_range(const __m256i v)
{
    __m256i msk_lo       =_mm256_set1_epi64x(0xFFFFFFFF);
    __m256d cnst2_32_dbl =_mm256_set1_pd(4294967296.0);                 /* 2^32                                                                    */

    __m256i v_lo         = _mm256_and_si256(v,msk_lo);                  /* extract the 32 lowest significant bits of v                             */
    __m256i v_hi         = _mm256_srli_epi64(v,32);                     /* 32 most significant bits of v. srai_epi64 doesn't exist                 */
    __m256i v_sign       = _mm256_srai_epi32(v,32);                     /* broadcast sign bit to the 32 most significant bits                      */
            v_hi         = _mm256_blend_epi32(v_hi,v_sign,0b10101010);  /* restore the correct sign of v_hi                                        */
    __m256d v_lo_dbl     = int64_to_double256(v_lo);                    /* v_lo is within specified range of int64_to_double                       */ 
    __m256d v_hi_dbl     = int64_to_double256(v_hi);                    /* v_hi is within specified range of int64_to_double                       */ 
            v_hi_dbl     = _mm256_mul_pd(cnst2_32_dbl,v_hi_dbl);        /* _mm256_mul_pd and _mm256_add_pd may compile to a single fma instruction */
    return _mm256_add_pd(v_hi_dbl,v_lo_dbl);                            /* rounding occurs if the integer doesn't exist as a double                */   
}


__m256d int64_to_double_based_on_cvtsi2sd(const __m256i v)
{   __m128d zero         = _mm_setzero_pd();                            /* to avoid uninitialized variables in_mm_cvtsi64_sd                       */
    __m128i v_lo         = _mm256_castsi256_si128(v);
    __m128i v_hi         = _mm256_extracti128_si256(v,1);
    __m128d v_0          = _mm_cvtsi64_sd(zero,_mm_cvtsi128_si64(v_lo));
    __m128d v_2          = _mm_cvtsi64_sd(zero,_mm_cvtsi128_si64(v_hi));
    __m128d v_1          = _mm_cvtsi64_sd(zero,_mm_extract_epi64(v_lo,1));
    __m128d v_3          = _mm_cvtsi64_sd(zero,_mm_extract_epi64(v_hi,1));
    __m128d v_01         = _mm_unpacklo_pd(v_0,v_1);
    __m128d v_23         = _mm_unpacklo_pd(v_2,v_3);
    __m256d v_dbl        = _mm256_castpd128_pd256(v_01);
            v_dbl        = _mm256_insertf128_pd(v_dbl,v_23,1);
    return v_dbl;
}


__m256d uint64_to_double_full_range(const __m256i v)                    
{
    __m256i msk_lo       =_mm256_set1_epi64x(0xFFFFFFFF);
    __m256d cnst2_32_dbl =_mm256_set1_pd(4294967296.0);                 /* 2^32                                                                    */

    __m256i v_lo         = _mm256_and_si256(v,msk_lo);                  /* extract the 32 lowest significant bits of v                             */
    __m256i v_hi         = _mm256_srli_epi64(v,32);                     /* 32 most significant bits of v                                           */
    __m256d v_lo_dbl     = uint64_to_double256(v_lo);                   /* v_lo is within specified range of uint64_to_double                      */ 
    __m256d v_hi_dbl     = uint64_to_double256(v_hi);                   /* v_hi is within specified range of uint64_to_double                      */ 
            v_hi_dbl     = _mm256_mul_pd(cnst2_32_dbl,v_hi_dbl);        
    return _mm256_add_pd(v_hi_dbl,v_lo_dbl);                            /* rounding may occur for inputs >2^52                                     */ 
}



int main(int argc, char **argv){
  int i;
  uint64_t j;
  __m256i j_4, j_inc;
  __m256d v, v_acc;
  double x[4];
  char test = argv[1][0];

  if (test=='A'){               /* test the conversions for several integer values                                       */
    j = 1ull;
    printf("\nint64_to_double_full_range\n");
    for (i = 0; i<30; i++){
      j_4= _mm256_set_epi64x(j-3,j+3,-j,j);
      v  = int64_to_double_full_range(j_4);
      _mm256_storeu_pd(x,v);
      printf("j =%21li    v =%23.1f    -v=%23.1f    v+3=%23.1f    v-3=%23.1f  \n",j,x[0],x[1],x[2],x[3]);
      j  = j*7ull;
    }

    j = 1ull;
    printf("\nint64_to_double_based_on_cvtsi2sd\n");
    for (i = 0; i<30; i++){
      j_4= _mm256_set_epi64x(j-3,j+3,-j,j);
      v  = int64_to_double_based_on_cvtsi2sd(j_4);
      _mm256_storeu_pd(x,v);
      printf("j =%21li    v =%23.1f    -v=%23.1f    v+3=%23.1f    v-3=%23.1f  \n",j,x[0],x[1],x[2],x[3]);
      j  = j*7ull;
    }

    j = 1ull;                       
    printf("\nuint64_to_double_full_range\n");
    for (i = 0; i<30; i++){
      j_4= _mm256_set_epi64x(j-3,j+3,j,j);
      v  = uint64_to_double_full_range(j_4);
      _mm256_storeu_pd(x,v);
      printf("j =%21lu    v =%23.1f   v+3=%23.1f    v-3=%23.1f \n",j,x[0],x[2],x[3]);
      j  = j*7ull;    
    }
  }
  else{
    j_4   = _mm256_set_epi64x(-123,-4004,-312313,-23412731);  
    j_inc = _mm256_set_epi64x(1,1,1,1);  
    v_acc = _mm256_setzero_pd();
    switch(test){

      case 'B' :{                  
        printf("\nLatency int64_to_double_cvtsi2sd()\n");      /* simple test to get a rough idea of the latency of int64_to_double_cvtsi2sd()     */
        for (i = 0; i<1000000000; i++){
          v  =int64_to_double_based_on_cvtsi2sd(j_4);
          j_4= _mm256_castpd_si256(v);                         /* cast without conversion, use output as an input in the next step                 */
        }
        _mm256_storeu_pd(x,v);
      }
      break;

      case 'C' :{                  
        printf("\nLatency int64_to_double_full_range()\n");    /* simple test to get a rough idea of the latency of int64_to_double_full_range()    */
        for (i = 0; i<1000000000; i++){
          v  = int64_to_double_full_range(j_4);
          j_4= _mm256_castpd_si256(v);
        }
        _mm256_storeu_pd(x,v);
      }
      break;

      case 'D' :{                  
        printf("\nThroughput int64_to_double_cvtsi2sd()\n");   /* simple test to get a rough idea of the throughput of int64_to_double_cvtsi2sd()   */
        for (i = 0; i<1000000000; i++){
          j_4   = _mm256_add_epi64(j_4,j_inc);                 /* each step a different input                                                       */
          v     = int64_to_double_based_on_cvtsi2sd(j_4);
          v_acc = _mm256_xor_pd(v,v_acc);                      /* use somehow the results                                                           */
        }
        _mm256_storeu_pd(x,v_acc);
      }
      break;

      case 'E' :{                  
        printf("\nThroughput int64_to_double_full_range()\n"); /* simple test to get a rough idea of the throughput of int64_to_double_full_range() */
        for (i = 0; i<1000000000; i++){
          j_4   = _mm256_add_epi64(j_4,j_inc);  
          v     = int64_to_double_full_range(j_4);
          v_acc = _mm256_xor_pd(v,v_acc);           
        }    
        _mm256_storeu_pd(x,v_acc);
      }
      break;

      default : {}
    }  
    printf("v =%23.1f    -v =%23.1f    v =%23.1f    -v =%23.1f  \n",x[0],x[1],x[2],x[3]);
  }

  return 0;
}

Фактическая производительность этих функций может зависеть от окружающего кода и генерации процессора.

Результаты синхронизации для преобразований 1e9 (шириной 256 бит) с простыми тестами B, C, D и E в приведенном выше коде для системы Intel Skylake i5 6500:

Latency experiment int64_to_double_based_on_cvtsi2sd()      (test B)  5.02 sec.
Latency experiment int64_to_double_full_range()             (test C)  3.77 sec.
Throughput experiment int64_to_double_based_on_cvtsi2sd()   (test D)  2.82 sec.
Throughput experiment int64_to_double_full_range()          (test E)  1.07 sec.

Разница в пропускной способности между int64_to_double_full_range() а также int64_to_double_based_on_cvtsi2sd() больше, чем я ожидал.

Источник

Я преобразовал ответ @JanWassenberg во внутреннюю форму для полного диапазона AVX и SSE2.doubleкint64_tконверсия. (Вы можете найти их внизу этого ответа).

Но затем я попробовал инвертировать алгоритм @wim, в котором 64-битное целое число разбивается на две 32-битные части, и похоже, что этот подход будет немного быстрее и намного проще.

AVX

      // Only works for inputs in the range: [-2^51, 2^51]
__m256i double_to_int64(__m256d x){
    x = _mm256_add_pd(x, _mm256_set1_pd(0x0018000000000000));
    return _mm256_sub_epi64(
        _mm256_castpd_si256(x),
        _mm256_castpd_si256(_mm256_set1_pd(0x0018000000000000))
    );
}
// Only works for inputs in the range: [-2^51, 2^51]
__m256d int64_to_double(__m256i x){
    x = _mm256_add_epi64(x, _mm256_castpd_si256(_mm256_set1_pd(0x0018000000000000)));
    return _mm256_sub_pd(_mm256_castsi256_pd(x), _mm256_set1_pd(0x0018000000000000));
}

__m256i my_mm256_cvtpd_epi64(__m256d v) {
    const __m256d k2_32inv_dbl = _mm256_set1_pd(1.0/4294967296.0); // 1 / 2^32
    const __m256d k2_32_dbl = _mm256_set1_pd(4294967296.0); // 2^32
    
    // Multiply by inverse instead of dividing.
    const __m256d v_hi_dbl = _mm256_mul_pd(v, k2_32inv_dbl);
    // Convert to integer.
    const __m256i v_hi = double_to_int64(v_hi_dbl);
    // Convert high32 integer to double and multiply by 2^32.
    const __m256d v_hi_int_dbl = _mm256_mul_pd(int64_to_double(v_hi), k2_32_dbl);
    // Subtract that from the original to get the remainder.
    const __m256d v_lo_dbl = _mm256_sub_pd(v, v_hi_int_dbl);
    // Convert to low32 integer.
    const __m256i v_lo = double_to_int64(v_lo_dbl);
    // Reconstruct integer from shifted high32 and remainder.
    return _mm256_add_epi64(_mm256_slli_epi64(v_hi, 32), v_lo);
}

SSE2

Мне нужна была реализация только для SSE2, поскольку я работаю со старыми машинами, не поддерживающими AVX.

      // Only works for inputs in the range: [-2^51, 2^51]
__m128i double_to_int64(__m128d x){
    x = _mm_add_pd(x, _mm_set1_pd(0x0018000000000000));
    return _mm_sub_epi64(
        _mm_castpd_si128(x),
        _mm_castpd_si128(_mm_set1_pd(0x0018000000000000))
    );
}
// Only works for inputs in the range: [-2^51, 2^51]
__m128d int64_to_double(__m128i x){
    x = _mm_add_epi64(x, _mm_castpd_si128(_mm_set1_pd(0x0018000000000000)));
    return _mm_sub_pd(_mm_castsi128_pd(x), _mm_set1_pd(0x0018000000000000));
}

__m128i my_mm_cvtpd_epi64(__m128d v) {
    const __m128d k2_32inv_dbl = _mm_set1_pd(1.0/4294967296.0); // 1 / 2^32
    const __m128d k2_32_dbl = _mm_set1_pd(4294967296.0); // 2^32
    
    // Multiply by inverse instead of dividing.
    const __m128d v_hi_dbl = _mm_mul_pd(v, k2_32inv_dbl);
    // Convert to integer.
    const __m128i v_hi = double_to_int64(v_hi_dbl);
    // Convert high32 integer to double and multiply by 2^32.
    const __m128d v_hi_int_dbl = _mm_mul_pd(int64_to_double(v_hi), k2_32_dbl);
    // Subtract that from the original to get the remainder.
    const __m128d v_lo_dbl = _mm_sub_pd(v, v_hi_int_dbl);
    // Convert to low32 integer.
    const __m128i v_lo = double_to_int64(v_lo_dbl);
    // Reconstruct integer from shifted high32 and remainder.
    return _mm_add_epi64(_mm_slli_epi64(v_hi, 32), v_lo);
}

Самое замечательное в них то, что их довольно легко понять. Кроме того, вам не нужна отдельная функция для беззнакового случая.

Скаляр

Учитываяcvttsd2siинструкция существует, может ли преобразование каждого элемента по отдельности быть быстрее?

Я прогнал все эти алгоритмы через llvm-mca (анализатор машинного кода LLVM). Он предоставляет число RThroughput, которое может дать хорошее представление о том, как что-то может работать, но вам следует выполнить реальные тесты, если вы хотите добиться максимальной производительности. (Например, тесты @wim не совпадают с некоторыми результатами здесь)

cvtpd_epi64

      __m128i my_mm_cvtpd_epi64(__m128d v) {
    double dd[2];
    _mm_storeu_pd(dd, v);

    int64_t ii[2];
    ii[0] = static_cast<int64_t>(dd[0]);
    ii[1] = static_cast<int64_t>(dd[1]);

    return _mm_loadu_si128(reinterpret_cast<__m128i*>(ii));
}

божий болт SSE

llvm-mca RЦифры пропускной способности (чем меньше, тем лучше)

Скалярные реализации выигрываютcvtepu64_pd.

Переводы ответа @JanWassenberg на необработанные встроенные функции:

Как отметил @MrUnbeliebable92, поведениеstatic_cast<int64_t>(double x)карты переполняютсяINT64_MIN, поэтому, чтобы соответствовать этому поведению, я смог удалить одинxorинструкция.

Мне также удалось реализовать еще несколько оптимизаций, например, отдельно сдвинуть неявный бит.

AVX2

      __m256i my_mm256_cvtpd_epi64(__m256d v) {
    const __m256i k51p3FF = _mm256_set1_epi64x(51 + 0x3FF);

    // Exponent indicates whether the number can be represented as uuint64_t.
    const __m256i biased_exp = _mm256_and_si256(
        _mm256_srli_epi64(_mm256_castpd_si256(v), 52),
        _mm256_set1_epi64x(0x7FF)
    );
    const __m256i mantissa = _mm256_and_si256(
        _mm256_castpd_si256(v),
        _mm256_set1_epi64x((1ULL << 52) - 1)
    );

    // Calculate left and right shifts to move mantissa into place.
    const __m256i shift_right = _mm256_subs_epu16(k51p3FF, biased_exp);
    const __m256i shift_left = _mm256_subs_epu16(biased_exp, k51p3FF);

    // Shift mantissa into place.
    const __m256i shifted = _mm256_srli_epi64(_mm256_srlv_epi64(
        _mm256_sllv_epi64(mantissa, shift_left),
        shift_right
    ), 1);
    // Include implicit 1-bit.
    const __m256i implicit_bit_shifted = _mm256_srlv_epi64(
        _mm256_sllv_epi64(_mm256_set1_epi64x(1ULL << 51), shift_left),
        shift_right
    );
    const __m256i magnitude = _mm256_or_si256(
        shifted,
        implicit_bit_shifted
    );

    // Fill each 64-bit part with sign bits.
    const __m256i sign_mask = _mm256_shuffle_epi32(
        _mm256_srai_epi32(_mm256_castpd_si256(v), 31),
        _MM_SHUFFLE(3, 3, 1, 1)
    );
    // Adjust for negative values.
    const __m256i sign_adjusted = _mm256_sub_epi64(
        _mm256_xor_si256(magnitude, sign_mask),
        sign_mask
    );

    // 0xFF is exp < 64
    const __m256i upper_bound_mask = _mm256_shuffle_epi32(_mm256_cmpgt_epi32(
        _mm256_set1_epi32(64 + 0x3FF),
        biased_exp
    ), _MM_SHUFFLE(2, 2, 0, 0));
    // Saturate overflow values to INT64_MIN.
    const __m256i bounded = _mm256_blendv_epi8(
        _mm256_set1_epi64x(INT64_MAX),
        sign_adjusted,
        upper_bound_mask
    );

    return bounded;
}

__m256i my_mm256_cvtpd_epu64(__m256d v) {
    const __m256i k51p3FF = _mm256_set1_epi64x(51 + 0x3FF);

    // Exponent indicates whether the number can be represented as uuint64_t.
    const __m256i biased_exp = _mm256_and_si256(
        _mm256_srli_epi64(_mm256_castpd_si256(v), 52),
        _mm256_set1_epi64x(0x7FF)
    );
    const __m256i mantissa = _mm256_and_si256(
        _mm256_castpd_si256(v),
        _mm256_set1_epi64x((1ULL << 52) - 1)
    );

    // Calculate left and right shifts to move mantissa into place.
    const __m256i shift_right = _mm256_subs_epu16(k51p3FF, biased_exp);
    const __m256i shift_left = _mm256_subs_epu16(biased_exp, k51p3FF);

    // Shift mantissa into place.
    const __m256i shifted = _mm256_srli_epi64(_mm256_srlv_epi64(
        _mm256_sllv_epi64(mantissa, shift_left),
        shift_right
    ), 1);
    // Include implicit 1-bit.
    const __m256i implicit_bit_shifted = _mm256_srlv_epi64(
        _mm256_sllv_epi64(_mm256_set1_epi64x(1ULL << 51), shift_left),
        shift_right
    );
    const __m256i magnitude = _mm256_or_si256(
        shifted,
        implicit_bit_shifted
    );

    // Fill each 64-bit part with sign bits.
    const __m256i sign_mask = _mm256_shuffle_epi32(
        _mm256_srai_epi32(_mm256_castpd_si256(v), 31),
        _MM_SHUFFLE(3, 3, 1, 1)
    );
    // Mask out negative values to 0.
    const __m256i lower_bounded = _mm256_andnot_si256(sign_mask, magnitude);

    // 0xFF is exp < 64
    const __m256i upper_bound_mask = _mm256_shuffle_epi32(_mm256_cmpgt_epi32(
        _mm256_set1_epi32(64 + 0x3FF),
        biased_exp
    ), _MM_SHUFFLE(2, 2, 0, 0));
    // Mask out overflow values to 0.
    const __m256i fully_bounded = _mm256_and_si256(lower_bounded, upper_bound_mask);
    
    return fully_bounded;
}

SSE2

      // Emulate mm_srlv_epi64 from AVX2
__m128i my_mm_srlv_epi64(__m128i a, __m128i count) {
    __m128i shift_low = _mm_srl_epi64(a, count);            // high 64 is garbage
    __m128i count_high = _mm_unpackhi_epi64(count, count);  // broadcast the high element
    __m128i shift_high = _mm_srl_epi64(a, count_high);      // low 64 is garbage

    // use movsd as a blend.
    return _mm_castpd_si128(_mm_move_sd(
        _mm_castsi128_pd(shift_high),
        _mm_castsi128_pd(shift_low)
    ));
}
// Emulate mm_sllv_epi64 from AVX2
__m128i my_mm_sllv_epi64(__m128i a, __m128i count) {
    __m128i shift_low = _mm_sll_epi64(a, count);            // high 64 is garbage
    __m128i count_high = _mm_unpackhi_epi64(count, count);  // broadcast the high element
    __m128i shift_high = _mm_sll_epi64(a, count_high);      // low 64 is garbage

    // use movsd as a blend.
    return _mm_castpd_si128(_mm_move_sd(
        _mm_castsi128_pd(shift_high),
        _mm_castsi128_pd(shift_low)
    ));
}

__m128i my_mm_cvtpd_epi64(__m128d v) {
    const __m128i k51p3FF = _mm_set1_epi64x(51 + 0x3FF);

    // Exponent indicates whether the number can be represented as uint64_t.
    const __m128i biased_exp = _mm_and_si128(
        _mm_srli_epi64(_mm_castpd_si128(v), 52),
        _mm_set1_epi64x(0x7FF)
    );
    const __m128i mantissa = _mm_and_si128(
        _mm_castpd_si128(v),
        _mm_set1_epi64x((1ULL << 52) - 1)
    );

    // Calculate left and right shifts to move mantissa into place.
    const __m128i shift_right = _mm_subs_epu16(k51p3FF, biased_exp);
    const __m128i shift_left = _mm_subs_epu16(biased_exp, k51p3FF);

    // Shift mantissa into place.
    const __m128i shifted = _mm_srli_epi64(my_mm_srlv_epi64(
        my_mm_sllv_epi64(mantissa, shift_left),
        shift_right
    ), 1);
    // Include implicit 1-bit.
    const __m128i implicit_bit_shifted = my_mm_srlv_epi64(
        my_mm_sllv_epi64(_mm_set1_epi64x(1ULL << 51), shift_left),
        shift_right
    );
    const __m128i magnitude = _mm_or_si128(
        shifted,
        implicit_bit_shifted
    );

    // Fill each 64-bit part with sign bits.
    const __m128i sign_mask = _mm_shuffle_epi32(
        _mm_srai_epi32(_mm_castpd_si128(v), 31),
        _MM_SHUFFLE(3, 3, 1, 1)
    );
    // Adjust for negative values.
    const __m128i sign_adjusted = _mm_sub_epi64(
        _mm_xor_si128(magnitude, sign_mask),
        sign_mask
    );

    // 0xFF is exp < 64
    const __m128i upper_bound_mask = _mm_shuffle_epi32(_mm_cmpgt_epi32(
        _mm_set1_epi32(64 + 0x3FF),
        biased_exp
    ), _MM_SHUFFLE(2, 2, 0, 0));
    // Saturate overflow values to INT64_MIN.
    const __m128i bounded = _mm_or_si128(
        _mm_and_si128(upper_bound_mask, sign_adjusted),
        _mm_andnot_si128(upper_bound_mask, _mm_set1_epi64x(INT64_MAX))
    );
    
    return bounded;
}

__m128i my_mm_cvtpd_epu64(__m128d v) {
    const __m128i k51p3FF = _mm_set1_epi64x(51 + 0x3FF);

    // Exponent indicates whether the number can be represented as uint64_t.
    const __m128i biased_exp = _mm_and_si128(
        _mm_srli_epi64(_mm_castpd_si128(v), 52),
        _mm_set1_epi64x(0x7FF)
    );
    const __m128i mantissa = _mm_and_si128(
        _mm_castpd_si128(v),
        _mm_set1_epi64x((1ULL << 52) - 1)
    );

    // Calculate left and right shifts to move mantissa into place.
    const __m128i shift_right = _mm_subs_epu16(k51p3FF, biased_exp);
    const __m128i shift_left = _mm_subs_epu16(biased_exp, k51p3FF);

    // Shift mantissa into place.
    const __m128i shifted = _mm_srli_epi64(my_mm_srlv_epi64(
        my_mm_sllv_epi64(mantissa, shift_left),
        shift_right
    ), 1);
    // Include implicit 1-bit.
    const __m128i implicit_bit_shifted = my_mm_srlv_epi64(
        my_mm_sllv_epi64(_mm_set1_epi64x(1ULL << 51), shift_left),
        shift_right
    );
    const __m128i magnitude = _mm_or_si128(
        shifted,
        implicit_bit_shifted
    );

    // Fill each 64-bit part with sign bits.
    const __m128i sign_mask = _mm_shuffle_epi32(
        _mm_srai_epi32(_mm_castpd_si128(v), 31),
        _MM_SHUFFLE(3, 3, 1, 1)
    );
    // Mask out negative values to 0.
    const __m128i lower_bounded = _mm_andnot_si128(sign_mask, magnitude);

    // 0xFF is exp < 64
    const __m128i upper_bound_mask = _mm_shuffle_epi32(_mm_cmpgt_epi32(
        _mm_set1_epi32(64 + 0x3FF),
        biased_exp
    ), _MM_SHUFFLE(2, 2, 0, 0));
    // Mask out overflow values to 0.
    const __m128i fully_bounded = _mm_and_si128(lower_bounded, upper_bound_mask);
    
    return fully_bounded;
}
Источник

Спасибо @mysticial и @wim за полнофункциональный i64->f64. Я придумал полное усечение f64->i64 для оболочки Highway SIMD .

Первая версия пыталась изменить режим округления, но Clang сортирует их и игнорирует ASM летучие памяти / куб.см затирает и даже атомный забор. Мне не ясно, как сделать это безопасным; NOINLINE работает, но вызывает много утечек.

Вторая версия (ссылка на Compiler Explorer) эмулирует перенормировку FP и оказывается быстрее в соответствии с llvm-mca (8-10 циклов rthroughput / total).

      // Full-range F64 -> I64 conversion
#include <hwy/highway.h>

namespace hwy {
namespace HWY_NAMESPACE {

HWY_API Vec256<int64_t> I64FromF64(Full256<int64_t> di, const Vec256<double> v) {
  const RebindToFloat<decltype(di)> dd;
  using VD = decltype(v);
  using VI = decltype(Zero(di));

  const VI k0 = Zero(di);
  const VI k1 = Set(di, 1);
  const VI k51 = Set(di, 51);

  // Exponent indicates whether the number can be represented as int64_t.
  const VI biased_exp = ShiftRight<52>(BitCast(di, v)) & Set(di, 0x7FF);
  const VI exp = biased_exp - Set(di, 0x3FF);
  const auto in_range = exp < Set(di, 63);

  // If we were to cap the exponent at 51 and add 2^52, the number would be in
  // [2^52, 2^53) and mantissa bits could be read out directly. We need to
  // round-to-0 (truncate), but changing rounding mode in MXCSR hits a
  // compiler reordering bug: https://gcc.godbolt.org/z/4hKj6c6qc . We instead
  // manually shift the mantissa into place (we already have many of the
  // inputs anyway).
  const VI shift_mnt = Max(k51 - exp, k0);
  const VI shift_int = Max(exp - k51, k0);
  const VI mantissa = BitCast(di, v) & Set(di, (1ULL << 52) - 1);
  // Include implicit 1-bit; shift by one more to ensure it's in the mantissa.
  const VI int52 = (mantissa | Set(di, 1ULL << 52)) >> (shift_mnt + k1);
  // For inputs larger than 2^52, insert zeros at the bottom.
  const VI shifted = int52 << shift_int;
  // Restore the one bit lost when shifting in the implicit 1-bit.
  const VI restored = shifted | ((mantissa & k1) << (shift_int - k1));

  // Saturate to LimitsMin (unchanged when negating below) or LimitsMax.
  const VI sign_mask = BroadcastSignBit(BitCast(di, v));
  const VI limit = Set(di, LimitsMax<int64_t>()) - sign_mask;
  const VI magnitude = IfThenElse(in_range, restored, limit);

  // If the input was negative, negate the integer (two's complement).
  return (magnitude ^ sign_mask) - sign_mask;
}

void Test(const double* pd, int64_t* pi) {
    Full256<int64_t> di;
    Full256<double> dd;
    for (size_t i = 0; i < 256; i += Lanes(di)) {
      Store(I64FromF64(di, Load(dd, pd + i)), di, pi + i);
    }
}

}
}

Если кто-то видит потенциал для упрощения алгоритма, оставьте, пожалуйста, комментарий.

Источник
Другие вопросы по тегам