c++ - ¿Cómo realizar eficientemente conversiones dobles/int64 con SSE/AVX?

performance floating-point

Question

SSE2 tiene instrucciones para convertir vectores entre flotantes de precisión simple y enteros de 32 bits.

• _mm_cvtps_epi32()
• _mm_cvtepi32_ps()

Pero no hay equivalentes para enteros de doble precisión y de 64 bits. En otras palabras, faltan:

• _mm_cvtpd_epi64()
• _mm_cvtepi64_pd()

Parece que AVX tampoco los tiene.

¿Cuál es la forma más eficiente de simular estos intrínsecos?

Answer 1

Esta respuesta es aproximadamente un entero de 64 bits a doble conversión, sin cortar esquinas. En una versión anterior de esta respuesta (ver párrafo Conversión rápida y precisa al dividir ... , a continuación), se demostró que es bastante eficiente dividir los enteros de 64 bits en un nivel bajo de 32 bits y uno de 32 bits. parte alta, convierta estas partes a doble y calcule low + high * 2^32 .

Los recuentos de instrucciones de estas conversiones fueron:

• int64_to_double_full_range 9 instrucciones (con mul y add como una fma )
• uint64_to_double_full_range 7 instrucciones (con mul y add como una fma )

Inspirado por la respuesta actualizada de Mysticial, con conversiones precisas mejor optimizadas, int64_t aún más la conversión int64_t a doble:

• int64_to_double_fast_precise : 5 instrucciones.
• uint64_to_double_fast_precise : 5 instrucciones.

La conversión int64_to_double_fast_precise requiere una instrucción menos que la solución de Mysticial. El código uint64_to_double_fast_precise es esencialmente idéntico a la solución de Mysticial (pero con un vpblendd lugar de vpblendw ). Se incluye aquí debido a sus similitudes con la conversión int64_to_double_fast_precise : las instrucciones son idénticas, solo las constantes son diferentes:

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

__m256d int64_to_double_fast_precise(const __m256i v)
/* Optimized full range int64_t to double conversion           */
/* Emulate _mm256_cvtepi64_pd()                                */
{
    __m256i magic_i_lo   = _mm256_set1_epi64x(0x4330000000000000);                /* 2^52               encoded as floating-point  */
    __m256i magic_i_hi32 = _mm256_set1_epi64x(0x4530000080000000);                /* 2^84 + 2^63        encoded as floating-point  */
    __m256i magic_i_all  = _mm256_set1_epi64x(0x4530000080100000);                /* 2^84 + 2^63 + 2^52 encoded as floating-point  */
    __m256d magic_d_all  = _mm256_castsi256_pd(magic_i_all);

    __m256i v_lo         = _mm256_blend_epi32(magic_i_lo, v, 0b01010101);         /* Blend the 32 lowest significant bits of v with magic_int_lo                                                   */
    __m256i v_hi         = _mm256_srli_epi64(v, 32);                              /* Extract the 32 most significant bits of v                                                                     */
            v_hi         = _mm256_xor_si256(v_hi, magic_i_hi32);                  /* Flip the msb of v_hi and blend with 0x45300000                                                                */
    __m256d v_hi_dbl     = _mm256_sub_pd(_mm256_castsi256_pd(v_hi), magic_d_all); /* Compute in double precision:                                                                                  */
    __m256d result       = _mm256_add_pd(v_hi_dbl, _mm256_castsi256_pd(v_lo));    /* (v_hi - magic_d_all) + v_lo  Do not assume associativity of floating point addition !!                        */
            return result;                                                        /* With gcc use -O3, then -fno-associative-math is default. Do not use -Ofast, which enables -fassociative-math! */
                                                                                  /* With icc use -fp-model precise                                                                                */
}


__m256d uint64_to_double_fast_precise(const __m256i v)                    
/* Optimized full range uint64_t to double conversion          */
/* This code is essentially identical to Mysticial''s solution. */
/* Emulate _mm256_cvtepu64_pd()                                */
{
    __m256i magic_i_lo   = _mm256_set1_epi64x(0x4330000000000000);                /* 2^52        encoded as floating-point  */
    __m256i magic_i_hi32 = _mm256_set1_epi64x(0x4530000000000000);                /* 2^84        encoded as floating-point  */
    __m256i magic_i_all  = _mm256_set1_epi64x(0x4530000000100000);                /* 2^84 + 2^52 encoded as floating-point  */
    __m256d magic_d_all  = _mm256_castsi256_pd(magic_i_all);

    __m256i v_lo         = _mm256_blend_epi32(magic_i_lo, v, 0b01010101);         /* Blend the 32 lowest significant bits of v with magic_int_lo                                                   */
    __m256i v_hi         = _mm256_srli_epi64(v, 32);                              /* Extract the 32 most significant bits of v                                                                     */
            v_hi         = _mm256_xor_si256(v_hi, magic_i_hi32);                  /* Blend v_hi with 0x45300000                                                                                    */
    __m256d v_hi_dbl     = _mm256_sub_pd(_mm256_castsi256_pd(v_hi), magic_d_all); /* Compute in double precision:                                                                                  */
    __m256d result       = _mm256_add_pd(v_hi_dbl, _mm256_castsi256_pd(v_lo));    /* (v_hi - magic_d_all) + v_lo  Do not assume associativity of floating point addition !!                        */
            return result;                                                        /* With gcc use -O3, then -fno-associative-math is default. Do not use -Ofast, which enables -fassociative-math! */
                                                                                  /* With icc use -fp-model precise                                                                                */
}


int main(){
    int i;
    uint64_t j;
    __m256i j_4;
    __m256d v;
    double x[4];
    double x0, x1, a0, a1;

    j = 0ull;
    printf("/nAccurate int64_to_double/n");
    for (i = 0; i < 260; i++){
        j_4= _mm256_set_epi64x(0, 0, -j, j);

        v  = int64_to_double_fast_precise(j_4);
        _mm256_storeu_pd(x,v);
        x0 = x[0];
        x1 = x[1];
        a0 = _mm_cvtsd_f64(_mm_cvtsi64_sd(_mm_setzero_pd(),j));
        a1 = _mm_cvtsd_f64(_mm_cvtsi64_sd(_mm_setzero_pd(),-j));
        printf(" j =%21li   v =%23.1f   v=%23.1f   -v=%23.1f   -v=%23.1f   d=%.1f   d=%.1f/n", j, x0, a0, x1, a1, x0-a0, x1-a1);
        j  = j+(j>>2)-(j>>5)+1ull;
    }

    j = 0ull;
    printf("/nAccurate uint64_to_double/n");
    for (i = 0; i < 260; i++){
        if (i==258){j=-1;}
        if (i==259){j=-2;}
        j_4= _mm256_set_epi64x(0, 0, -j, j);

        v  = uint64_to_double_fast_precise(j_4);
        _mm256_storeu_pd(x,v);
        x0 = x[0];
        x1 = x[1];
        a0 = (double)((uint64_t)j);
        a1 = (double)((uint64_t)-j);
        printf(" j =%21li   v =%23.1f   v=%23.1f   -v=%23.1f   -v=%23.1f   d=%.1f   d=%.1f/n", j, x0, a0, x1, a1, x0-a0, x1-a1);
        j  = j+(j>>2)-(j>>5)+1ull;
    }
    return 0;
}

Las conversiones pueden fallar si las opciones de optimización matemática inseguras están habilitadas. Con gcc, -O3 es seguro, pero -Ofast puede conducir a resultados incorrectos, porque aquí no podemos asumir la asociatividad de la adición de punto flotante (lo mismo vale para las conversiones de Mysticial). Con icc use -fp-model precise .

Conversión rápida y precisa al dividir los enteros de 64 bits en una parte baja de 32 bits y una parte alta de 32 bits.

Suponemos que tanto la entrada entera como la salida doble están en registros AVX de 256 bits de ancho. Se consideran dos enfoques:

1. int64_to_double_based_on_cvtsi2sd() : como se sugiere en los comentarios sobre la pregunta, use cvtsi2sd 4 veces junto con algunos cvtsi2sd de datos. Desafortunadamente, tanto cvtsi2sd como las instrucciones de barajado de datos necesitan el puerto de ejecución 5. Esto limita el rendimiento de este enfoque.

2. int64_to_double_full_range() : podemos usar el método de conversión rápida de Mysticial dos veces para obtener una conversión precisa para el rango entero de 64 bits. El entero de 64 bits se divide en una parte baja de 32 bits y una parte alta de 32 bits, de manera similar a las respuestas a esta pregunta: ¿Cómo realizar la conversión uint32 / float con SSE? . Cada una de estas piezas es adecuada para el entero de Mysticial para duplicar la conversión. Finalmente, la parte alta se multiplica por 2 ^ 32 y se agrega a la parte baja. La conversión firmada es un poco más complicada que la conversión sin firmar ( uint64_to_double_full_range() ), porque srai_epi64() no existe.

Código:

#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;
}

El rendimiento real de estas funciones puede depender del código circundante y la generación de la CPU.

Resultados de tiempo para conversiones 1e9 (256 bits de ancho) con pruebas simples B, C, D y E en el código anterior, en un sistema 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.

La diferencia en el rendimiento entre int64_to_double_full_range() y int64_to_double_based_on_cvtsi2sd() es mayor de lo que esperaba.

Answer 2

Si está dispuesto a cortar esquinas, las conversiones double <-> int64 se pueden hacer en solo dos instrucciones:

• Si no te importa el infinito o el NaN .
• Para double <-> int64_t , solo le importan los valores en el rango [-2^51, 2^51] .
• Para double <-> uint64_t , solo le interesan los valores en el rango [0, 2^52) .

doble -> 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))
    );
}

doble -> 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 -> doble

//  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 -> doble

//  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));
}

Comportamiento de redondeo:

• Para la conversión double -> uint64_t , el redondeo funciona correctamente siguiendo el modo de redondeo actual. (que suele ser redondo a par)
• Para la conversión double -> int64_t , el redondeo seguirá el modo de redondeo actual para todos los modos, excepto el truncamiento. Si el modo de redondeo actual es el truncamiento (redondeado hacia cero), en realidad se redondeará hacia el infinito negativo.

¿Como funciona?

A pesar de que este truco tiene solo 2 instrucciones, no se explica por sí solo.

La clave es reconocer que para punto flotante de doble precisión, los valores en el rango [2^52, 2^53) tienen el "lugar binario" justo debajo del bit más bajo de la mantisa. En otras palabras, si pone a cero el exponente y los bits de signo, la mantisa se convierte precisamente en la representación entera.

Para convertir x de double -> uint64_t , agrega el número mágico M que es el valor de coma flotante de 2^52 . Esto pone x en el rango "normalizado" de [2^52, 2^53) y redondea convenientemente los bits de parte fraccionaria.

Ahora todo lo que queda es eliminar los 12 bits superiores. Esto se hace fácilmente enmascarándolo. La forma más rápida es reconocer que esos 12 bits superiores son idénticos a los de M Entonces, en lugar de introducir una constante de máscara adicional, simplemente podemos restar o XOR por M XOR tiene más rendimiento.

La conversión de uint64_t -> double es simplemente lo contrario de este proceso. Vuelve a agregar los bits exponentes de M Luego, normalice el número restando M en coma flotante.

Las conversiones de enteros con signo son un poco más complicadas ya que necesita lidiar con la extensión de signo del complemento de 2. Los dejaré como ejercicio para el lector.

Relacionado: Un método rápido para redondear un doble a un int de 32 bits explicado

Gama completa int64 -> doble:

Después de muchos años, finalmente tuve una necesidad de esto.

• 5 instrucciones para uint64_t -> double
• 6 instrucciones para int64_t -> double

uint64_t -> doble

__m128d uint64_to_double_full(__m128i x){
    __m128i xH = _mm_srli_epi64(x, 32);
    xH = _mm_or_si128(xH, _mm_castpd_si128(_mm_set1_pd(19342813113834066795298816.)));          //  2^84
    __m128i xL = _mm_blend_epi16(x, _mm_castpd_si128(_mm_set1_pd(0x0010000000000000)), 0xcc);   //  2^52
    __m128d f = _mm_sub_pd(_mm_castsi128_pd(xH), _mm_set1_pd(19342813118337666422669312.));     //  2^84 + 2^52
    return _mm_add_pd(f, _mm_castsi128_pd(xL));
}

int64_t -> doble

__m128d int64_to_double_full(__m128i x){
    __m128i xH = _mm_srai_epi32(x, 16);
    xH = _mm_blend_epi16(xH, _mm_setzero_si128(), 0x33);
    xH = _mm_add_epi64(xH, _mm_castpd_si128(_mm_set1_pd(442721857769029238784.)));              //  3*2^67
    __m128i xL = _mm_blend_epi16(x, _mm_castpd_si128(_mm_set1_pd(0x0010000000000000)), 0x88);   //  2^52
    __m128d f = _mm_sub_pd(_mm_castsi128_pd(xH), _mm_set1_pd(442726361368656609280.));          //  3*2^67 + 2^52
    return _mm_add_pd(f, _mm_castsi128_pd(xL));
}

Funcionan para todo el rango de 64 bits y se redondean correctamente al comportamiento de redondeo actual.

Estas son las respuestas de wim similares a continuación, pero con optimizaciones más abusivas. Como tal, descifrarlos también se dejará como un ejercicio para el lector.