Come eseguire efficientemente conversioni double / int64 con SSE / AVX?

SSE2 ha istruzioni per convertire i vettori tra float a precisione singola e interi a 32 bit.

  • _mm_cvtps_epi32()
  • _mm_cvtepi32_ps()

Ma non ci sono equivalenti per la precisione doppia e gli interi a 64 bit. In altre parole, mancano:

  • _mm_cvtpd_epi64()
  • _mm_cvtepi64_pd()

Sembra che AVX non li abbia neanche.

Qual è il modo più efficace per simulare questi elementi intrinseci?

Se sei disposto a tagliare gli angoli, le conversioni double <-> int64 possono essere eseguite in due sole istruzioni:

  • Se non ti interessa infinito o NaN .
  • Per double <-> int64_t , ti preoccupi solo dei valori nell’intervallo [-2^51, 2^51] .
  • Per il double <-> uint64_t , ti interessano solo i valori nell’intervallo [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 -> double

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

Comportamento di arrotondamento:

  • Per la conversione double -> uint64_t , l’arrotondamento funziona correttamente seguendo la modalità di arrotondamento corrente. (che di solito è rotondo-pari)
  • Per la conversione double -> int64_t , l’arrotondamento seguirà la modalità di arrotondamento corrente per tutte le modalità tranne il troncamento. Se la modalità di arrotondamento corrente è troncata (arrotondata a zero), verrà effettivamente arrotondata all’infinito negativo.

Come funziona?

Nonostante questo trucco sia solo 2 istruzioni, non è completamente auto-esplicativo.

La chiave è riconoscere che per i valori a virgola mobile a precisione doppia, i valori nell’intervallo [2^52, 2^53) hanno il “posto binario” appena sotto il bit più basso della mantissa. In altre parole, se azzerate l’esponente e firmate i bit, la mantissa diventa esattamente la rappresentazione intera.

Per convertire x da double -> uint64_t , aggiungi il numero magico M che è il valore a virgola mobile di 2^52 . Questo inserisce x nell’intervallo “normalizzato” di [2^52, 2^53) e arrotonda convenientemente i bit della parte frazionaria.

Ora tutto ciò che rimane è rimuovere i 12 bit superiori. Questo è facile da mascherare. Il modo più veloce è riconoscere che quei 12 bit superiori sono identici a quelli di M Quindi, piuttosto che introdurre una maschera aggiuntiva costante, possiamo semplicemente sottrarre o XOR per M XOR ha più throughput.

La conversione da uint64_t -> double è semplicemente il contrario di questo processo. Aggiungete i bit esponenti di M Quindi non normalizzare il numero sottraendo M in virgola mobile.

Le conversioni intere con segno sono leggermente più complesse poiché è necessario gestire l’estensione di segno del complemento a 2. Lascerò quelli come esercizio per il lettore.

Correlato: Un metodo veloce per arrotondare un doppio a un int a 32 bit spiegato

Questa risposta è circa 64 bit interi per la doppia conversione, senza tagliare gli angoli. Supponiamo che sia l’intero sia l’uscita doppia siano in registri AVX larghi 256 bit. Sono stati considerati due approcci:

  1. int64_to_double_based_on_cvtsi2sd() : come suggerito nei commenti sulla domanda, usa cvtsi2sd 4 volte insieme ad alcuni dati shuffling. Sfortunatamente entrambe le istruzioni di cvtsi2sd e di shuffling dei dati richiedono la porta di esecuzione 5. Ciò limita le prestazioni di questo approccio.

  2. int64_to_double_full_range() : possiamo utilizzare il metodo di conversione veloce di Mysticial due volte per ottenere una conversione accurata dell’intero intervallo intero a 64 bit. Il numero intero a 64 bit è diviso in una parte a 32 bit e una parte a 32 bit, analogamente alle risposte a questa domanda: Come eseguire la conversione uint32 / float con SSE? . Ognuno di questi pezzi è adatto per l’intero di Mysticial alla doppia conversione. Infine, la parte alta viene moltiplicata per 2 ^ 32 e aggiunta alla parte bassa. La conversione firmata è un po ‘più complicata della conversione senza segno ( uint64_to_double_full_range() ), perché srai_epi64() non esiste.

Codice:

 #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); /* 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.1fv =%23.1f -v =%23.1f \n",x[0],x[1],x[2],x[3]); } return 0; } 

Le prestazioni effettive di queste funzioni possono dipendere dal codice circostante e dalla generazione della CPU.

Risultati di cronometraggio per 1e9 conversioni (larghezza 256 bit) con semplici test B, C, D ed E nel codice sopra, su un sistema i5 6500 intel skylake:

 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 differenza di throughput tra int64_to_double_full_range() e int64_to_double_based_on_cvtsi2sd() è maggiore di quanto mi aspettassi.