Implementación más rápida de la función exponencial utilizando SSE

optimization vectorization

Question

Estoy buscando una aproximación de la función de exponente que opera en el elemento SSE. A saber: __m128 exp( __m128 x ) .

Tengo una implementación que es rápida pero parece ser muy baja en precisión:

static inline __m128 FastExpSse(__m128 x)
{
__m128 a = _mm_set1_ps(12102203.2f); // (1 << 23) / ln(2)
__m128i b = _mm_set1_epi32(127 * (1 << 23) - 486411);
__m128 m87 = _mm_set1_ps(-87);
// fast exponential function, x should be in [-87, 87]
__m128 mask = _mm_cmpge_ps(x, m87);
__m128i tmp = _mm_add_epi32(_mm_cvtps_epi32(_mm_mul_ps(a, x)), b);
return _mm_and_ps(_mm_castsi128_ps(tmp), mask);
}

¿Alguien podría tener una implementación con mayor precisión pero tan rápida (o más rápida)?

Sería feliz si escribiera en estilo C.

Gracias.

Answer 1

El código C a continuación es una traducción al intrínseco SSE de un algoritmo que utilicé en una respuesta anterior a una pregunta similar.

La idea básica es transformar el cálculo de la función exponencial estándar en el cálculo de una potencia de 2: expf (x) = exp2f (x / logf (2.0f)) = exp2f (x * 1.44269504) . Dividimos t = x * 1.44269504 en un entero i y una fracción f , de modo que t = i + f y 0 <= f <= 1 . Ahora podemos calcular 2 ^f con una aproximación polinómica, luego escalar el resultado en 2 ⁱ sumando i al campo exponente del resultado de punto flotante de precisión simple.

Un problema que existe con una implementación SSE es que queremos calcular i = floorf (t) , pero no hay una forma rápida de calcular la función floor() . Sin embargo, observamos que para números positivos, floor(x) == trunc(x) , y que para números negativos, floor(x) == trunc(x) - 1 , excepto cuando x es un entero negativo. Sin embargo, dado que la aproximación central puede manejar un valor f de 1.0f , usar la aproximación para argumentos negativos es inofensivo. SSE proporciona una instrucción para convertir operandos de punto flotante de precisión simple en enteros con truncamiento, por lo que esta solución es eficiente.

Peter Cordes señala que SSE4.1 admite una función de piso rápida _mm_floor_ps() , por lo que a continuación también se muestra una variante que usa SSE4.1. No todas las cadenas de herramientas predefinen automáticamente la macro __SSE4_1__ cuando la generación de código SSE 4.1 está habilitada, pero gcc sí.

Compiler Explorer (Godbolt) muestra que gcc 7.2 compila el código siguiente en dieciséis instrucciones para SSE simple y doce instrucciones para SSE 4.1.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <emmintrin.h>
#ifdef __SSE4_1__
#include <smmintrin.h>
#endif
/* max. rel. error = 1.72863156e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
__m128 t, f, e, p, r;
__m128i i, j;
__m128 l2e = _mm_set1_ps (1.442695041f); /* log2(e) */
__m128 c0 = _mm_set1_ps (0.3371894346f);
__m128 c1 = _mm_set1_ps (0.657636276f);
__m128 c2 = _mm_set1_ps (1.00172476f);
/* exp(x) = 2^i * 2^f; i = floor (log2(e) * x), 0 <= f <= 1 */
t = _mm_mul_ps (x, l2e); /* t = log2(e) * x */
#ifdef __SSE4_1__
e = _mm_floor_ps (t); /* floor(t) */
i = _mm_cvtps_epi32 (e); /* (int)floor(t) */
#else /* __SSE4_1__*/
i = _mm_cvttps_epi32 (t); /* i = (int)t */
j = _mm_srli_epi32 (_mm_castps_si128 (x), 31); /* signbit(t) */
i = _mm_sub_epi32 (i, j); /* (int)t - signbit(t) */
e = _mm_cvtepi32_ps (i); /* floor(t) ~= (int)t - signbit(t) */
#endif /* __SSE4_1__*/
f = _mm_sub_ps (t, e); /* f = t - floor(t) */
p = c0; /* c0 */
p = _mm_mul_ps (p, f); /* c0 * f */
p = _mm_add_ps (p, c1); /* c0 * f + c1 */
p = _mm_mul_ps (p, f); /* (c0 * f + c1) * f */
p = _mm_add_ps (p, c2); /* p = (c0 * f + c1) * f + c2 ~= 2^f */
j = _mm_slli_epi32 (i, 23); /* i << 23 */
r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
return r;
}
int main (void)
{
union {
float f[4];
unsigned int i[4];
} arg, res;
double relerr, maxrelerr = 0.0;
int i, j;
__m128 x, y;
float start[2] = {-0.0f, 0.0f};
float finish[2] = {-87.33654f, 88.72283f};
for (i = 0; i < 2; i++) {
arg.f[0] = start[i];
arg.i[1] = arg.i[0] + 1;
arg.i[2] = arg.i[0] + 2;
arg.i[3] = arg.i[0] + 3;
do {
memcpy (&x, &arg, sizeof(x));
y = fast_exp_sse (x);
memcpy (&res, &y, sizeof(y));
for (j = 0; j < 4; j++) {
double ref = exp ((double)arg.f[j]);
relerr = fabs ((res.f[j] - ref) / ref);
if (relerr > maxrelerr) {
printf ("arg=% 15.8e res=%15.8e ref=%15.8e err=%15.8e/n",
arg.f[j], res.f[j], ref, relerr);
maxrelerr = relerr;
}
}
arg.i[0] += 4;
arg.i[1] += 4;
arg.i[2] += 4;
arg.i[3] += 4;
} while (fabsf (arg.f[3]) < fabsf (finish[i]));
}
printf ("maximum relative errror = %15.8e/n", maxrelerr);
return EXIT_SUCCESS;
}

Un diseño alternativo para fast_sse_exp() extrae la parte entera del argumento ajustado x / log(2) en modo redondeado al más cercano, utilizando la técnica bien conocida de agregar la constante de conversión "mágica" 1.5 * 2 ²³ para forzar el redondeo en la posición de bit correcta, luego restando el mismo número nuevamente. Esto requiere que el modo de redondeo SSE vigente durante la adición sea "redondear al más cercano o incluso", que es el valor predeterminado. wim señaló en los comentarios que algunos compiladores pueden optimizar la suma y resta de la conversión constante cvt como redundante cuando se utiliza una optimización agresiva, lo que interfiere con la funcionalidad de esta secuencia de código, por lo que se recomienda inspeccionar el código de máquina generado. El intervalo de aproximación para el cálculo de 2 ^f ahora está centrado alrededor de cero, ya que -0.5 <= f <= 0.5 , que requiere una aproximación de núcleo diferente.

/* max. rel. error <= 1.72860465e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
__m128 t, f, p, r;
__m128i i, j;
const __m128 l2e = _mm_set1_ps (1.442695041f); /* log2(e) */
const __m128 cvt = _mm_set1_ps (12582912.0f); /* 1.5 * (1 << 23) */
const __m128 c0 = _mm_set1_ps (0.238428936f);
const __m128 c1 = _mm_set1_ps (0.703448006f);
const __m128 c2 = _mm_set1_ps (1.000443142f);
/* exp(x) = 2^i * 2^f; i = rint (log2(e) * x), -0.5 <= f <= 0.5 */
t = _mm_mul_ps (x, l2e); /* t = log2(e) * x */
r = _mm_sub_ps (_mm_add_ps (t, cvt), cvt); /* r = rint (t) */
f = _mm_sub_ps (t, r); /* f = t - rint (t) */
i = _mm_cvtps_epi32 (t); /* i = (int)t */
p = c0; /* c0 */
p = _mm_mul_ps (p, f); /* c0 * f */
p = _mm_add_ps (p, c1); /* c0 * f + c1 */
p = _mm_mul_ps (p, f); /* (c0 * f + c1) * f */
p = _mm_add_ps (p, c2); /* p = (c0 * f + c1) * f + c2 ~= exp2(f) */
j = _mm_slli_epi32 (i, 23); /* i << 23 */
r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
return r;
}

El algoritmo para el código en la pregunta parece estar tomado del trabajo de Nicol N. Schraudolph, que explota hábilmente la naturaleza semi-logarítmica de los formatos de punto flotante binario IEEE-754:

NN Schraudolph. "Una aproximación rápida y compacta de la función exponencial". Computación neuronal , 11 (4), mayo de 1999, pp.853-862.

Después de eliminar el código de sujeción del argumento, se reduce a solo tres instrucciones SSE. La constante de corrección "mágica" 486411 no es óptima para minimizar el error relativo máximo en todo el dominio de entrada. Basado en una búsqueda binaria simple, el valor 298765 parece ser superior, reduciendo el error relativo máximo para FastExpSse() a 3.56e-2 frente al error relativo máximo de 1.73e-3 para fast_exp_sse() .

/* max. rel. error = 3.55959567e-2 on [-87.33654, 88.72283] */
__m128 FastExpSse (__m128 x)
{
__m128 a = _mm_set1_ps (12102203.0f); /* (1 << 23) / log(2) */
__m128i b = _mm_set1_epi32 (127 * (1 << 23) - 298765);
__m128i t = _mm_add_epi32 (_mm_cvtps_epi32 (_mm_mul_ps (a, x)), b);
return _mm_castsi128_ps (t);
}

El algoritmo de Schraudolph utiliza básicamente la aproximación lineal 2 ^f ~ = 1.0 + f para f en [0,1], y su precisión podría mejorarse agregando un término cuadrático. La parte inteligente del enfoque de Schraudolph es calcular 2 ⁱ * 2 ^f sin separar explícitamente la parte entera i = floor(x * 1.44269504) de la fracción. No veo forma de extender ese truco a una aproximación cuadrática, pero ciertamente se puede combinar el cálculo floor() de Schraudolph con la aproximación cuadrática utilizada anteriormente:

/* max. rel. error <= 1.72886892e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
__m128 f, p, r;
__m128i t, j;
const __m128 a = _mm_set1_ps (12102203.0f); /* (1 << 23) / log(2) */
const __m128i m = _mm_set1_epi32 (0xff800000); /* mask for integer bits */
const __m128 ttm23 = _mm_set1_ps (1.1920929e-7f); /* exp2(-23) */
const __m128 c0 = _mm_set1_ps (0.3371894346f);
const __m128 c1 = _mm_set1_ps (0.657636276f);
const __m128 c2 = _mm_set1_ps (1.00172476f);
t = _mm_cvtps_epi32 (_mm_mul_ps (a, x));
j = _mm_and_si128 (t, m); /* j = (int)(floor (x/log(2))) << 23 */
t = _mm_sub_epi32 (t, j);
f = _mm_mul_ps (ttm23, _mm_cvtepi32_ps (t)); /* f = (x/log(2)) - floor (x/log(2)) */
p = c0; /* c0 */
p = _mm_mul_ps (p, f); /* c0 * f */
p = _mm_add_ps (p, c1); /* c0 * f + c1 */
p = _mm_mul_ps (p, f); /* (c0 * f + c1) * f */
p = _mm_add_ps (p, c2); /* p = (c0 * f + c1) * f + c2 ~= 2^f */
r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
return r;
}

Answer 2

Hay un documento sobre la creación de versiones rápidas de estas ecuaciones (tanh, cosh, artanh, sinh, etc.):

http://ijeais.org/wp-content/uploads/2018/07/IJAER180702.pdf "Creación de una implementación en línea optimizada del compilador de Intel Svml Simd Intrinsics"

su ecuación 6, en la página 9 es muy similar a la respuesta de @NicSchraudolph

Answer 3

Se puede obtener un buen aumento de precisión en mi algoritmo (implementación FastExpSse en la respuesta anterior) a costa de una resta entera y una división de punto flotante utilizando FastExpSse (x / 2) / FastExpSse (-x / 2) en lugar de FastExpSse (X). El truco aquí es establecer el parámetro de desplazamiento (298765 arriba) a cero para que las aproximaciones lineales por partes en el numerador y el denominador se alineen para darle una cancelación sustancial de errores. Enróllelo en una sola función:

__m128 BetterFastExpSse (__m128 x)
{
const __m128 a = _mm_set1_ps ((1 << 22) / float(M_LN2)); // to get exp(x/2)
const __m128i b = _mm_set1_epi32 (127 * (1 << 23)); // NB: zero shift!
__m128i r = _mm_cvtps_epi32 (_mm_mul_ps (a, x));
__m128i s = _mm_add_epi32 (b, r);
__m128i t = _mm_sub_epi32 (b, r);
return _mm_div_ps (_mm_castsi128_ps (s), _mm_castsi128_ps (t));
}

(No soy un experto en hardware, ¿qué tan malo es el rendimiento de la división aquí?)

Si necesita exp (x) solo para obtener y = tanh (x) (por ejemplo, para redes neuronales), use FastExpSse con desplazamiento a cero de la siguiente manera:

a = FastExpSse(x);
b = FastExpSse(-x);
y = (a - b)/(a + b);

para obtener el mismo tipo de beneficio de cancelación de error. La función logística funciona de manera similar, utilizando FastExpSse (x / 2) / (FastExpSse (x / 2) + FastExpSse (-x / 2)) con desplazamiento cero. (Esto es solo para mostrar el principio: obviamente no desea evaluar FastExpSse varias veces aquí, sino que lo convierte en una sola función a lo largo de las líneas de BetterFastExpSse anteriores).

Desarrollé una serie de aproximaciones de orden superior a partir de esto, cada vez más precisas pero también más lentas. Inédito pero feliz de colaborar si alguien quiere darles un giro.

Y finalmente, para divertirse: use en marcha atrás para obtener FastLogSse. Encadenar eso con FastExpSse le brinda al operador y la cancelación de errores, y aparece una función de potencia increíblemente rápida ...

Answer 4

Volviendo a mis notas desde entonces, exploré formas de mejorar la precisión sin usar la división. Utilicé el mismo truco de reinterpretar como flotante, pero apliqué una corrección polinómica a la mantisa que se calculó esencialmente en aritmética de punto fijo de 16 bits (la única forma de hacerlo rápido en ese entonces).

El resp. Cúbico. Las versiones cuárticas le dan 4 resp. 5 dígitos significativos de precisión. No tenía sentido aumentar el orden más allá de eso, ya que el ruido de la aritmética de baja precisión comienza a ahogar el error de la aproximación polinómica. Aquí están las versiones simples de C:

#include <stdint.h>
float fastExp3(register float x) // cubic spline approximation
{
union { float f; int32_t i; } reinterpreter;
reinterpreter.i = (int32_t)(12102203.0f*x) + 127*(1 << 23);
int32_t m = (reinterpreter.i >> 7) & 0xFFFF; // copy mantissa
// empirical values for small maximum relative error (8.34e-5):
reinterpreter.i +=
((((((((1277*m) >> 14) + 14825)*m) >> 14) - 79749)*m) >> 11) - 626;
return reinterpreter.f;
}
float fastExp4(register float x) // quartic spline approximation
{
union { float f; int32_t i; } reinterpreter;
reinterpreter.i = (int32_t)(12102203.0f*x) + 127*(1 << 23);
int32_t m = (reinterpreter.i >> 7) & 0xFFFF; // copy mantissa
// empirical values for small maximum relative error (1.21e-5):
reinterpreter.i += (((((((((((3537*m) >> 16)
+ 13668)*m) >> 18) + 15817)*m) >> 14) - 80470)*m) >> 11);
return reinterpreter.f;
}

El cuarto obedece (fastExp4 (0f) == 1f), que puede ser importante para los algoritmos de iteración de punto fijo.

¿Qué tan eficientes son estas secuencias enteras de multiplicación-cambio-suma en SSE? En arquitecturas donde la aritmética de flotación es igual de rápida, se podría usar eso, reduciendo el ruido aritmético. Esto esencialmente produciría extensiones cúbicas y cuárticas de la respuesta anterior de @ njuffa.