Implementazione efficiente di log2 (__ m256d) in AVX2

Lo __m256d _mm256_log2_pd (__m256d a) non è disponibile su altri compilatori oltre a Intel, e dicono che le sue prestazioni sono disabilitate sui processori AMD. Ci sono alcune implementazioni su Internet di AVX Log intrinsec (_mm256_log_ps) mancanti in g ++ – 4.8? e le librerie matematiche SIMD per SSE e AVX , tuttavia sembrano essere più SSE di AVX2. C’è anche la libreria vettoriale di Agner Fog , tuttavia è una grande libreria con molte più cose che si limitano al vettore log2, quindi dall’implementazione in esso è difficile capire le parti essenziali solo per l’operazione vector log2.

Quindi qualcuno può spiegare come implementare l’operazione log2() per un vettore di 4 numeri double modo efficiente? __m256d _mm256_log2_pd (__m256d a) piace ciò che fa __m256d _mm256_log2_pd (__m256d a) , ma disponibile per altri compilatori e abbastanza efficiente per processori AMD e Intel.

EDIT: Nel mio caso specifico corrente, i numeri sono probabilità tra 0 e 1, e il logaritmo è usato per il calcolo entropico: la negazione della sum su tutto i di P[i]*log(P[i]) . La gamma di esponenti in virgola mobile per P[i] è grande, quindi i numeri possono essere vicini a 0. Non sono sicuro dell’accuratezza, quindi considererei qualsiasi soluzione che inizi con 30 bit di mantissa, in particolare una soluzione accordabile è preferibile .

EDIT2: ecco la mia implementazione fino ad ora, basata su “Serie più efficienti” da https://en.wikipedia.org/wiki/Logarithm#Power_series . Come può essere migliorato? (sono desiderati miglioramenti delle prestazioni e della precisione)

 namespace { const __m256i gDoubleExpMask = _mm256_set1_epi64x(0x7ffULL << 52); const __m256i gDoubleExp0 = _mm256_set1_epi64x(1023ULL << 52); const __m256i gTo32bitExp = _mm256_set_epi32(0, 0, 0, 0, 6, 4, 2, 0); const __m128i gExpNormalizer = _mm_set1_epi32(1023); //TODO: some 128-bit variable or two 64-bit variables here? const __m256d gCommMul = _mm256_set1_pd(2.0 / 0.693147180559945309417); // 2.0/ln(2) const __m256d gCoeff1 = _mm256_set1_pd(1.0 / 3); const __m256d gCoeff2 = _mm256_set1_pd(1.0 / 5); const __m256d gCoeff3 = _mm256_set1_pd(1.0 / 7); const __m256d gCoeff4 = _mm256_set1_pd(1.0 / 9); const __m256d gVect1 = _mm256_set1_pd(1.0); } __m256d __vectorcall Log2(__m256d x) { const __m256i exps64 = _mm256_srli_epi64(_mm256_and_si256(gDoubleExpMask, _mm256_castpd_si256(x)), 52); const __m256i exps32_avx = _mm256_permutevar8x32_epi32(exps64, gTo32bitExp); const __m128i exps32_sse = _mm256_castsi256_si128(exps32_avx); const __m128i normExps = _mm_sub_epi32(exps32_sse, gExpNormalizer); const __m256d expsPD = _mm256_cvtepi32_pd(normExps); const __m256d y = _mm256_or_pd(_mm256_castsi256_pd(gDoubleExp0), _mm256_andnot_pd(_mm256_castsi256_pd(gDoubleExpMask), x)); // Calculate t=(y-1)/(y+1) and t**2 const __m256d tNum = _mm256_sub_pd(y, gVect1); const __m256d tDen = _mm256_add_pd(y, gVect1); const __m256d t = _mm256_div_pd(tNum, tDen); const __m256d t2 = _mm256_mul_pd(t, t); // t**2 const __m256d t3 = _mm256_mul_pd(t, t2); // t**3 const __m256d terms01 = _mm256_fmadd_pd(gCoeff1, t3, t); const __m256d t5 = _mm256_mul_pd(t3, t2); // t**5 const __m256d terms012 = _mm256_fmadd_pd(gCoeff2, t5, terms01); const __m256d t7 = _mm256_mul_pd(t5, t2); // t**7 const __m256d terms0123 = _mm256_fmadd_pd(gCoeff3, t7, terms012); const __m256d t9 = _mm256_mul_pd(t7, t2); // t**9 const __m256d terms01234 = _mm256_fmadd_pd(gCoeff4, t9, terms0123); const __m256d log2_y = _mm256_mul_pd(terms01234, gCommMul); const __m256d log2_x = _mm256_add_pd(log2_y, expsPD); return log2_x; } 

Finora la mia implementazione ha dato 405 268 490 operazioni al secondo, e sembra precisa fino all’ottava cifra. La performance è misurata con la seguente funzione:

 #include  #include  #include  #include  // ... Log2() implementation here const int64_t cnLogs = 100 * 1000 * 1000; void BenchmarkLog2Vect() { __m256d sums = _mm256_setzero_pd(); auto start = std::chrono::high_resolution_clock::now(); for (int64_t i = 1; i <= cnLogs; i += 4) { const __m256d x = _mm256_set_pd(double(i+3), double(i+2), double(i+1), double(i)); const __m256d logs = Log2(x); sums = _mm256_add_pd(sums, logs); } auto elapsed = std::chrono::high_resolution_clock::now() - start; double nSec = 1e-6 * std::chrono::duration_cast(elapsed).count(); double sum = sums.m256d_f64[0] + sums.m256d_f64[1] + sums.m256d_f64[2] + sums.m256d_f64[3]; printf("Vect Log2: %.3lf Ops/sec calculated %.3lf\n", cnLogs / nSec, sum); } 

Confrontando i risultati di Logaritmo in C ++ e assemblaggio , l’implementazione del vettore corrente è 4 volte più veloce di std::log2() e 2,5 volte più veloce di std::log() .

In particolare, viene utilizzata la seguente formula di approssimazione: Termini di approssimazione - visivi inserisci la descrizione dell'immagine qui