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

algorithm floating-point

Question

El __m256d _mm256_log2_pd (__m256d a) no está disponible en otros compiladores que no sean Intel, y dicen que su rendimiento está limitado en los procesadores AMD. ¿Hay algunas implementaciones en Internet referidas en AVX log intrinsics (_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 está la biblioteca de vectores de Agner Fog , sin embargo, es una biblioteca grande que tiene muchas más cosas que solo el vector log2, por lo que, desde la implementación, es difícil descubrir las partes esenciales solo para la operación del vector log2.

Entonces, ¿alguien puede explicar cómo implementar la operación log2() para un vector de 4 números double eficiente? Es decir, lo que hace __m256d _mm256_log2_pd (__m256d a) , pero está disponible para otros compiladores y es razonablemente eficiente tanto para procesadores AMD como 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 suma sobre todo i de P[i]*log(P[i]) . El rango de exponentes de punto flotante para P[i] es grande, por lo que los números pueden estar cerca de 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 proporciona 405 268 490 operaciones por segundo, y parece precisa hasta el octavo dígito. El rendimiento se mide con la siguiente función:

#include <chrono>
#include <cmath>
#include <cstdio>
#include <immintrin.h>

// ... 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<std::chrono::microseconds>(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);
}

En comparación con los resultados de Logarithm en C ++ y ensamblaje , la implementación de vectores 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 utiliza la siguiente fórmula de aproximación:

Answer 1

Finalmente, este es mi mejor resultado, que en Ryzen 1800X @ 3.6GHz proporciona alrededor de 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 : vea al final cómo aumentar el rendimiento a 0.87 mil millones de logaritmos por segundo.

Casos especiales : los números negativos, el infinito negativo y los NaN s con 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 bit de signo positivo dan como resultado un logaritmo de alrededor de 1024. Si no le gusta cómo se tratan los casos especiales, una opción es agregar un código que los verifique 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<<cnLog2TblBits); i++) {
    const uint64_t iZp = (uint64_t(cZeroExp) << 52)
      | (uint64_t(i) << (52 - cnLog2TblBits)) | (1ULL << (52 - cnLog2TblBits - 1));
    const double zp = *reinterpret_cast<const double*>(&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 enfoque de tabla de búsqueda y un polinomio de primer grado, descrito principalmente en Wikipedia (el enlace está en los comentarios del código). Puedo permitirme asignar 8 KB de caché L1 aquí (que es la mitad de 16 KB 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 reduciendo cnLog2TblBits a, por ejemplo, 5 a expensas de disminuir la precisión del cálculo del logaritmo.

O para mantener alta la precisión, puede aumentar el número de términos polinómicos 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 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);

Luego comente // Leading integer part for the logarithm y el resto sin cambios siga.

Normalmente no necesita tantos términos, incluso para una tabla de pocos bits, solo proporcioné los coeficientes y los 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, debes experimentar lo que se adapte a tus necesidades.

Cuantos menos términos polinómicos calcule, obviamente, más rápidos serán los cálculos.

EDITAR : esta pregunta ¿En qué situación el AVX2 recolectaría 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 tanto, el rendimiento aumenta a 0.87 mil millones de logaritmos por segundo. Dejo la implementación actual tal como es porque es más idiomática y debería ser más rápida una vez que las CPU comienzan a realizar correctamente las operaciones de gather (con la fusión como lo hacen las GPU).

EDIT2 : en la CPU Ryzen (pero no en Intel) puede obtener un poco más de velocidad (aproximadamente 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)));

Answer 2

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, 1.0 a 2.0, por lo que un polinomio para log2(mantissa) solo tiene que ajustarse sobre ese rango muy limitado. (O, de manera equivalente, mantisa = 0.5 a 1.0, y cambie la constante de corrección de sesgo de exponente en 1).

Una expansión de la serie Taylor es un buen punto de partida para los coeficientes, pero generalmente desea minimizar el error absoluto máximo (o error relativo) en 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 que el error positivo máximo casi coincida con el error negativo máximo. Entonces puede hacer lo que se llama un ajuste minimax de los coeficientes.

Si es importante que su función evalúe log2(1.0) a exactamente 0.0 , puede organizar que eso suceda usando mantissa-1.0 como su polinomio, y sin coeficiente constante. 0.0 ^ n = 0.0 . Esto también mejora en gran medida el error relativo para las entradas cercanas a 1.0, incluso si el error absoluto sigue siendo pequeño.

¿Qué tan preciso necesita que sea y sobre qué rango de entradas? Como de costumbre, existe 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 polinomial más (y reajustando los coeficientes), o eliminando un poco de evitación de errores de redondeo.

La implementación de VCL de Agner Fog tiene como objetivo una precisión muy alta, utilizando trucos para evitar errores 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 repositorio oficial para VCL, pero ese se actualiza regularmente desde los lanzamientos de tarball aguas arriba).

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

(Ese es el verdadero truco para implementar log() : solo necesita 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 se integra 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 la precisión de float completa, y mucho menos el double , pero podría ajustarse 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 portar la función SSE2 de JRF a AVX2 + FMA (y especialmente AVX512 con _mm512_getexp_ps y _mm512_getmant_ps ), una vez que lo ajusté con cuidado. (Era parte de un proyecto comercial, así que no creo que pueda publicar el código). Una implementación rápida aproximada para float fue exactamente lo que quería.

En mi caso de uso, cada jrf_fastlog() era independiente, por lo que la ejecución de OOO ocultó 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 de ILP más alto que utiliza la función polynomial_5() VCL ( "Estrin''s esquema " , que hace algunas multiplicaciones no FMA antes de las FMA, lo que da como resultado instrucciones más totales).

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

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

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

• extraer 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).

  Adelante ajuste 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 polinómica para log(x) que sea precisa alrededor de x = 1.0. (Para el double , el log_d() VCL usa una relación de dos polinomios de quinto orden. @ Harold dice que a menudo esto es bueno para la precisión . Una división mezclada con una gran cantidad de FMA no suele afectar el rendimiento, pero tiene una latencia más alta que una FMA. Usar vrcpps + una iteración de Newton-Raphson suele ser más lento que simplemente usar vdivps en hardware moderno. Usar una relación también crea más ILP al evaluar dos polinomios de orden inferior en paralelo, en lugar de un polinomio de orden superior, y puede menor latencia general frente a una cadena dep larga para un polinomio de alto orden (que también acumularía un error de redondeo significativo a lo largo de esa cadena larga).

Luego agregue exponent + polynomial_approx_log(mantissa) para obtener el resultado final log () . VCL hace esto en múltiples 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 el material ln2 2 pasos y simplemente usar VM_LN2 si no apunta a una precisión de 0.5 o 1 ulp (o lo que sea que esta función realmente proporcione; IDK).

Supongo que la parte x - 0.5*x2 es realmente un término polinomial adicional. Esto es lo que quise decir con log base e estando integrado: necesitaría un coeficiente en esos términos, o deshacerse de esa línea y volver a ajustar los coeficientes polinómicos para log2. No puede simplemente multiplicar todos los coeficientes polinómicos por una constante.

Después de eso, verifica si hay subflujo, desbordamiento o denormal, y se ramifica si algún elemento en el vector necesita un procesamiento especial para producir un NaN o -Inf adecuado 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 comprobación antes de que la rama tome 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 polinómica, y hacer una corrección mínima de los coeficientes en lugar de simplemente usar una serie de Taylor expansión.

• http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html un enfoque interesante: escribe- uint32_t un float a uint32_t , y convierte ese entero a 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 en 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 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 JRF fastlog 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 ).