Implementación eficiente de log2 (__ m256d) en AVX2

El archivo __m256d _mm256_log2_pd (__m256d a) no está disponible en otros comstackdores distintos a Intel, y dicen que su rendimiento se ve perjudicado en los procesadores AMD. Hay algunas implementaciones en Internet que se mencionan en los intrínsecos de registro de AVX (_mm256_log_ps) que faltan en g ++ – 4.8? y bibliotecas matemáticas SIMD para SSE y AVX , sin embargo, parecen ser más SSE que AVX2. También hay una biblioteca de vectores de Agner Fog , sin embargo, es una gran biblioteca que tiene muchas más cosas que solo vector log2, así que desde la implementación en ella es difícil descubrir las partes esenciales para la operación vector log2.

Entonces, ¿alguien puede explicar cómo implementar la operación log2() para un vector de 4 números double eficiente? Es decir, como lo hace __m256d _mm256_log2_pd (__m256d a) , pero está disponible para otros comstackdores y es razonablemente eficiente para los procesadores AMD e Intel.

EDITAR: En mi caso específico actual, los números son probabilidades entre 0 y 1, y el logaritmo se usa para el cálculo de entropía: la negación de la sum sobre todo i de P[i]*log(P[i]) . El rango de exponentes de coma flotante para P[i] es grande, por lo que los números pueden ser cercanos a 0. No estoy seguro de la precisión, por lo que consideraría cualquier solución que comience con 30 bits de mantisa, especialmente una solución sintonizable. .

EDIT2: aquí está mi implementación hasta ahora, basada en “Series más eficientes” de https://en.wikipedia.org/wiki/Logarithm#Power_series . ¿Cómo puede ser mejorado? (Se desean mejoras de rendimiento y precisión)

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

Hasta ahora, mi implementación da 405 268 490 operaciones por segundo, y parece ser precisa hasta el octavo dígito. El rendimiento se mide con la siguiente función:

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

Comparando los resultados de Logaritmo en C ++ y ensamblado , la implementación del vector actual es 4 veces más rápida que std::log2() y 2.5 veces más rápida que std::log() .

Específicamente, se usa la siguiente fórmula de aproximación: Términos aproximados - visual enter image description here

La estrategia habitual se basa en el log(a*b) = log(a) + log(b) identidad log(a*b) = log(a) + log(b) , o en este caso log2( 2^exponent * mantissa) ) = log2( 2^exponent ) + log2(mantissa) . O simplificando, exponent + log2(mantissa) . La mantisa tiene un rango muy limitado, de 1.0 a 2.0, por lo que un polinomio para log2(mantissa) solo tiene que caber en un rango muy limitado. (O de forma equivalente, mantisa = 0.5 a 1.0, y cambie la constante de corrección de polarización del exponente por 1).

Una expansión de la serie Taylor es un buen punto de partida para los coeficientes, pero normalmente se desea minimizar el error máximo absoluto (o error relativo) sobre ese rango específico, y los coeficientes de la serie Taylor probablemente tengan un valor atípico menor o mayor en ese rango , en lugar de tener el error positivo máximo que casi coincide con el error negativo máximo. Entonces puedes hacer lo que se llama un ajuste minimax de los coeficientes.

Si es importante que tu función evalúe log2(1.0) a exactamente 0.0 , puedes organizar que eso suceda utilizando realmente mantissa-1.0 como tu polinomio, y sin coeficiente constante. 0.0 ^ n = 0.0 . Esto mejora enormemente el error relativo para las entradas cercanas a 1.0 también, incluso si el error absoluto es aún pequeño.


¿Cuán preciso lo necesita y sobre qué rango de entradas? Como de costumbre, hay una compensación entre precisión y velocidad, pero afortunadamente es bastante fácil moverse a lo largo de esa escala, por ejemplo, agregando un término polinómico más (y reajustando los coeficientes), o descartando algunos errores de redondeo.

La implementación de VCL de Agner Fog apunta a una precisión muy alta, utilizando trucos para evitar el error de redondeo al evitar cosas que podrían resultar en agregar un número pequeño y grande cuando sea posible. Esto oscurece un poco el diseño básico.

(Agner no tiene un repo oficial para VCL, pero ese se actualiza regularmente desde las versiones preliminares de tarball).


Para obtener un log() float más aproximado y más rápido log() , consulte la implementación del polinomio en http://jrfonseca.blogspot.ca/2008/09/fast-sse2-pow-tables-or-polynomials.html . Deja fuera MUCHOS trucos adicionales para ganar precisión que utiliza VCL, por lo que es más fácil de entender. Utiliza una aproximación polinomial para la mantisa en el rango de 1.0 a 2.0.

(Ese es el verdadero truco para las implementaciones de log() : solo necesitas un polinomio que funcione en un rango pequeño).

Ya solo hace log2 lugar de log , a diferencia de los VCL, donde log-base-e está integrado en las constantes y cómo las usa. Leerlo es probablemente un buen punto de partida para comprender las implementaciones exponent + polynomial(mantissa) de log() .

Incluso la versión de mayor precisión no es de precisión de float completa, ni mucho menos double , pero podría adaptarse a un polinomio con más términos. O aparentemente una relación de dos polinomios funciona bien; eso es lo que usa VCL para el double .

Obtuve excelentes resultados al transferir la función SSE2 de JRF a AVX2 + FMA (y especialmente a AVX512 con _mm512_getexp_ps y _mm512_getmant_ps ), una vez que la sintonicé cuidadosamente. (Era parte de un proyecto comercial, así que no creo que pueda publicar el código). Una implementación aproximada rápida para float era exactamente lo que quería.

En mi caso de uso, cada jrf_fastlog() era independiente, por lo que la ejecución de OOO escondía muy bien la latencia de FMA, y ni siquiera valía la pena usar el método de evaluación polinomial de latencia más corta con mayor ILP que utiliza la función polynomial_5() VCL ( “Estrin’s esquema ” , que hace algunas multiplicaciones no FMA antes de las FMA, lo que resulta en más instrucciones totales).


Si su proyecto puede usar código con licencia de GPL y desea una alta precisión, debe usar VCL directamente. Solo son encabezados, por lo que solo las partes que usas realmente serán parte de tus binarios.

Las funciones doble y flotante de log de VCL están en vectormath_exp.h . Hay dos partes principales del algoritmo:

  • extrae los bits de exponente y convierte ese entero de nuevo en un flotante (después de ajustar el sesgo que usa IEEE FP).

  • extrae la mantisa y OR en algunos bits de exponente para obtener un vector de valores double en el rango [0.5, 1.0) . (O (0.5, 1.0] , lo olvido).

    Ajusta esto con if(mantissa <= SQRT2*0.5) { mantissa += mantissa; exponent++;} if(mantissa <= SQRT2*0.5) { mantissa += mantissa; exponent++;} , y luego mantissa -= 1.0 .

    Use una aproximación polinomial para log(x) que sea precisa alrededor de x = 1.0. (Para el double , log_d() VCL usa una relación de dos polinomios de quinto orden. @Harold dice que esto a menudo es bueno para la precisión . Una división mezclada con muchas FMA generalmente no perjudica el rendimiento, pero sí tiene una mayor latencia que una FMA. Usar vrcpps + una iteración de Newton-Raphson es típicamente más lenta que usar vdivps en hardware moderno. Usar una razón también crea más ILP al evaluar dos polinomios de orden inferior en paralelo, en lugar de un polinomio de alto orden, y puede menor latencia global frente a una larga cadena de depósito para un polinomio de orden alto (que también acumularía un error de redondeo significativo a lo largo de esa larga cadena).

A continuación, agregue el exponent + polynomial_approx_log(mantissa) para obtener el resultado final de log () . VCL hace esto en varios pasos para reducir el error de redondeo. ln2_lo + ln2_hi = ln(2) . Se divide en una constante pequeña y grande para reducir el error de redondeo.

 // res is the polynomial(adjusted_mantissa) result // fe is the float exponent // x is the adjusted_mantissa. x2 = x*x; res = mul_add(fe, ln2_lo, res); // res += fe * ln2_lo; res += nmul_add(x2, 0.5, x); // res += x - 0.5 * x2; res = mul_add(fe, ln2_hi, res); // res += fe * ln2_hi; 

Puede soltar las cosas de 2 pasos de ln2 y simplemente usar VM_LN2 si no apunta a una precisión de 0.5 o 1 ulp (o lo que esta función realmente proporcione, IDK).

La parte x - 0.5*x2 es realmente un término polinomial adicional, supongo. Esto es lo que quise decir con logaritmo e echado en el horno: necesitaría un coeficiente en esos términos, o deshacerse de esa línea y volver a ajustar los coeficientes polinomiales para log2. No se puede simplemente multiplicar todos los coeficientes polinomiales por una constante.

Después de eso, comprueba si hay desbordamiento, desbordamiento o desnormal, y twigs si algún elemento en el vector necesita un procesamiento especial para producir un NaN o -Inf apropiado en lugar de cualquier basura que obtengamos del polinomio + exponente. Si se sabe que sus valores son finitos y positivos, puede comentar esta parte y obtener una aceleración significativa (incluso la verificación antes de la bifurcación requiere varias instrucciones).


Otras lecturas:

  • http://gallium.inria.fr/blog/fast-vectorizable-math-approx/ algunas cosas sobre cómo evaluar el error relativo y absoluto en una aproximación polinomial, y hacer una corrección minimax de los coeficientes en lugar de solo usar una serie de Taylor expansión.

  • http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html un enfoque interesante: teclea un float en uint32_t y convierte ese entero para float . Como los flotantes IEEE binary32 almacenan el exponente en bits más altos que la mantisa, el float resultante representa principalmente el valor del exponente, escalado por 1 << 23 , pero también contiene información de la mantisa.

    Luego usa una expresión con un par de coeficientes para arreglar las cosas y obtener una aproximación de log() . Incluye una división por (constant + mantissa) para corregir la contaminación de la mantisa al convertir el patrón de bits float en float . Descubrí que una versión vectorizada de eso era más lenta y menos precisa con AVX2 en HSW y SKL que en el registro rápido de JRF con polinomios de cuarto orden. (Especialmente cuando se usa como parte de un arcsinh rápido que también usa la unidad de división para vsqrtps ).

Finalmente, este es mi mejor resultado que en Ryzen 1800X @ 3.6GHz da aproximadamente 0.8 mil millones de logaritmos por segundo (200 millones de vectores de 4 logaritmos en cada uno) en un solo hilo, y es preciso hasta unos pocos bits en la mantisa. Spoiler : ver al final cómo boost el rendimiento a 0.87 mil millones de logaritmos por segundo.

Casos especiales : los números negativos, el infinito negativo y los NaN con el bit de signo negativo se tratan como si estuvieran muy cerca de 0 (dan como resultado algunos valores de “logaritmo” negativos grandes). El infinito positivo y los NaN con un bit de signo positivo dan como resultado un logaritmo alrededor de 1024. Si no le gusta cómo se tratan los casos especiales, una opción es agregar un código que los revise y haga lo que más le convenga. Esto hará que el cálculo sea más lento.

 namespace { // The limit is 19 because we process only high 32 bits of doubles, and out of // 20 bits of mantissa there, 1 bit is used for rounding. constexpr uint8_t cnLog2TblBits = 10; // 1024 numbers times 8 bytes = 8KB. constexpr uint16_t cZeroExp = 1023; const __m256i gDoubleNotExp = _mm256_set1_epi64x(~(0x7ffULL << 52)); const __m256d gDoubleExp0 = _mm256_castsi256_pd(_mm256_set1_epi64x(1023ULL << 52)); const __m256i cAvxExp2YMask = _mm256_set1_epi64x( ~((1ULL << (52-cnLog2TblBits)) - 1) ); const __m256d cPlusBit = _mm256_castsi256_pd(_mm256_set1_epi64x( 1ULL << (52 - cnLog2TblBits - 1))); const __m256d gCommMul1 = _mm256_set1_pd(2.0 / 0.693147180559945309417); // 2.0/ln(2) const __m256i gHigh32Permute = _mm256_set_epi32(0, 0, 0, 0, 7, 5, 3, 1); const __m128i cSseMantTblMask = _mm_set1_epi32((1 << cnLog2TblBits) - 1); const __m128i gExpNorm0 = _mm_set1_epi32(1023); // plus |cnLog2TblBits|th highest mantissa bit double gPlusLog2Table[1 << cnLog2TblBits]; } // anonymous namespace void InitLog2Table() { for(uint32_t i=0; i<(1<(&iZp); const double l2zp = std::log2(zp); gPlusLog2Table[i] = l2zp; } } __m256d __vectorcall Log2TblPlus(__m256d x) { const __m256d zClearExp = _mm256_and_pd(_mm256_castsi256_pd(gDoubleNotExp), x); const __m256d z = _mm256_or_pd(zClearExp, gDoubleExp0); const __m128i high32 = _mm256_castsi256_si128(_mm256_permutevar8x32_epi32( _mm256_castpd_si256(x), gHigh32Permute)); // This requires that x is non-negative, because the sign bit is not cleared before // computing the exponent. const __m128i exps32 = _mm_srai_epi32(high32, 20); const __m128i normExps = _mm_sub_epi32(exps32, gExpNorm0); // Compute y as approximately equal to log2(z) const __m128i indexes = _mm_and_si128(cSseMantTblMask, _mm_srai_epi32(high32, 20 - cnLog2TblBits)); const __m256d y = _mm256_i32gather_pd(gPlusLog2Table, indexes, /*number of bytes per item*/ 8); // Compute A as z/exp2(y) const __m256d exp2_Y = _mm256_or_pd( cPlusBit, _mm256_and_pd(z, _mm256_castsi256_pd(cAvxExp2YMask))); // Calculate t=(A-1)/(A+1). Both numerator and denominator would be divided by exp2_Y const __m256d tNum = _mm256_sub_pd(z, exp2_Y); const __m256d tDen = _mm256_add_pd(z, exp2_Y); // Compute the first polynomial term from "More efficient series" of https://en.wikipedia.org/wiki/Logarithm#Power_series const __m256d t = _mm256_div_pd(tNum, tDen); const __m256d log2_z = _mm256_fmadd_pd(t, gCommMul1, y); // Leading integer part for the logarithm const __m256d leading = _mm256_cvtepi32_pd(normExps); const __m256d log2_x = _mm256_add_pd(log2_z, leading); return log2_x; } 

Utiliza una combinación de aproximación de tabla de búsqueda y un polinomio de primer grado, descrito principalmente en Wikipedia (el enlace está en los comentarios del código). Me puedo permitir asignar 8KB de caché L1 aquí (que es la mitad de 16KB de caché L1 disponible por núcleo lógico), porque el cálculo del logaritmo es realmente el cuello de botella para mí y no hay mucho más que necesite caché L1.

Sin embargo, si necesita más caché L1 para las otras necesidades, puede disminuir la cantidad de caché utilizada por el algoritmo de logaritmo al reducir cnLog2TblBits a, por ejemplo, 5 a expensas de disminuir la precisión del cálculo del logaritmo.

O para mantener la precisión alta, puede boost la cantidad de términos polinomiales agregando:

 namespace { // ... 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 gCoeff5 = _mm256_set1_pd(1.0 / 11); } 

Y luego cambiando la cola de Log2TblPlus() después de la const __m256d t = _mm256_div_pd(tNum, tDen); línea 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 t11 = _mm256_mul_pd(t9, t2); // t**11 const __m256d terms012345 = _mm256_fmadd_pd(gCoeff5, t11, terms01234); const __m256d log2_z = _mm256_fmadd_pd(terms012345, gCommMul1, y); 

A continuación, comente // Leading integer part for the logarithm Siga la // Leading integer part for the logarithm y el rest sin cambios.

Normalmente no necesitas tantos términos, incluso para una tabla de pocos bits, solo proporcioné los coeficientes y cálculos para referencia. Es probable que si cnLog2TblBits==5 , no necesitará nada más allá de los terms012 . Pero no he hecho tales mediciones, necesita experimentar lo que se adapte a sus necesidades.

Cuantos menos términos polinómicos calcules, obviamente, más rápidos son los cálculos.


EDITAR : esta pregunta ¿En qué situación reuniría AVX2 las instrucciones sería más rápido que cargar los datos individualmente? sugiere que puede obtener una mejora en el rendimiento si

 const __m256d y = _mm256_i32gather_pd(gPlusLog2Table, indexes, /*number of bytes per item*/ 8); 

es reemplazado por

 const __m256d y = _mm256_set_pd(gPlusLog2Table[indexes.m128i_u32[3]], gPlusLog2Table[indexes.m128i_u32[2]], gPlusLog2Table[indexes.m128i_u32[1]], gPlusLog2Table[indexes.m128i_u32[0]]); 

Para mi implementación, ahorra aproximadamente 1.5 ciclos, reduciendo el recuento total de ciclos para calcular 4 logaritmos de 18 a 16.5, por lo que el rendimiento aumenta a 0.87 mil millones de logaritmos por segundo. Dejo la implementación actual tal como está, porque es más idiomática y debería ser más rápida una vez que las CPU comiencen a hacer operaciones de gather correctamente (con coalescencia como lo hacen las GPU).

EDIT2 : en la CPU Ryzen (pero no en Intel) puede obtener un poco más de velocidad (alrededor de 0,5 ciclos) reemplazando

 const __m128i high32 = _mm256_castsi256_si128(_mm256_permutevar8x32_epi32( _mm256_castpd_si256(x), gHigh32Permute)); 

con

  const __m128 hiLane = _mm_castpd_ps(_mm256_extractf128_pd(x, 1)); const __m128 loLane = _mm_castpd_ps(_mm256_castpd256_pd128(x)); const __m128i high32 = _mm_castps_si128(_mm_shuffle_ps(loLane, hiLane, _MM_SHUFFLE(3, 1, 3, 1)));