¿Cómo realizar conversiones dobles / int64 de manera eficiente con SSE / AVX?

SSE2 tiene instrucciones para convertir vectores entre flotantes de precisión simple e 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?

Si está dispuesto a double <-> int64 conversiones double <-> int64 se pueden realizar en solo dos instrucciones:

  • Si no te importa el infinito o NaN .
  • Para double <-> int64_t , solo te importan los valores en el rango [-2^51, 2^51] .
  • Para double <-> uint64_t , solo le importan 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)) ); } 

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 -> 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 generalmente es 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 (redondo hacia cero), en realidad se redondeará hacia el infinito negativo.

¿Como funciona?

A pesar de que este truco tiene solo 2 instrucciones, no es completamente autoexplicativo.

La clave es reconocer que para la coma flotante de precisión doble, 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 , agregue el número mágico M que es el valor de coma flotante de 2^52 . Esto pone a x en el rango “normalizado” de [2^52, 2^53) y redondea convenientemente los bits parciales fraccionarios.

Ahora todo lo que queda es eliminar los 12 bits superiores. Esto se hace fácilmente enmascarando. La manera 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, podemos simplemente restar o XOR por M XOR tiene más rendimiento.

La conversión de uint64_t -> double es simplemente el reverso de este proceso. Vuelve a agregar los bits de exponente de M Luego, anula la normalización del número restando M en coma flotante.

Las conversiones de entero 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

Esta respuesta es de 64 bits de entero a doble conversión, sin cortar esquinas. 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. Lamentablemente, tanto cvtsi2sd como las instrucciones de mezcla 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 alta de 32 bits y una parte alta de 32 bits, al igual que en las respuestas a esta pregunta: ¿Cómo realizar la conversión de uint32 / float con SSE? . Cada una de estas piezas es adecuada para el número entero de Mysticial para la conversión doble. 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 complicida que la conversión sin firma ( uint64_to_double_full_range() ), porque srai_epi64() no existe.

Código:

 #include  #include  #include  /* 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); /* restre 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.1fv =%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 1e9 conversiones (256 bit 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.