Hay una pregunta existente "Promedio de 3 enteros largos" que se ocupa específicamente del cálculo eficiente del promedio de tres enteros con signo .
Sin embargo, el uso de números enteros sin signo permite optimizaciones adicionales que no se aplican al escenario cubierto en la pregunta anterior. Esta pregunta trata sobre el cálculo eficiente del promedio de tres enteros sin signo , donde el promedio se redondea hacia cero, es decir, en términos matemáticos, quiero calcular ⌊ (a + b + c) / 3 ⌋.
Una forma sencilla de calcular este promedio es
avg = a / 3 + b / 3 + c / 3 + (a % 3 + b % 3 + c % 3) / 3; En primer orden, los compiladores de optimización modernos transformarán las divisiones en multiplicaciones con un recíproco más un cambio, y las operaciones de módulo en una multiplicación hacia atrás y una resta, donde la multiplicación hacia atrás puede usar un idioma scale_add disponible en muchas arquitecturas, por ejemplo, lea en x86_64, add con lsl #n en ARM, iscadd en GPU NVIDIA.
Al tratar de optimizar lo anterior de una manera genérica adecuada para muchas plataformas comunes, observo que, por lo general, el costo de las operaciones con enteros está en la relación lógico ≤ ( agregar | sub ) ≤ cambio ≤ escala_agregar ≤ mul . El costo aquí se refiere a toda la latencia, las limitaciones de rendimiento y el consumo de energía. Cualquier diferencia de este tipo se vuelve más pronunciada cuando el tipo entero procesado es más ancho que el ancho del registro nativo, por ejemplo, cuando se procesan datos uint64_t en un procesador de 32 bits.
Por lo tanto, mi estrategia de optimización fue minimizar el número de instrucciones y reemplazar las operaciones "costosas" por "baratas" cuando fuera posible, sin aumentar la presión del registro y manteniendo el paralelismo aprovechable para procesadores fuera de servicio amplios.
La primera observación es que podemos reducir una suma de tres operandos en una suma de dos operandos aplicando primero un CSA (carry save sumador) que produce un valor de suma y un valor de acarreo, donde el valor de acarreo tiene el doble del peso de la suma valor. El costo de un CSA basado en software es de cinco s lógicos en la mayoría de los procesadores. Algunos procesadores, como las GPU NVIDIA, tienen una instrucción LOP3 que puede calcular una expresión lógica arbitraria de tres operandos de una sola vez, en cuyo caso CSA se condensa en dos LOP3 (nota: todavía he convencido al compilador CUDA para que emita esos dos LOP3 s; actualmente produce cuatro LOP3 s!).
La segunda observación es que debido a que estamos calculando el módulo de la división por 3, no necesitamos una multiplicación inversa para calcularlo. En su lugar, podemos usar dividend % 3 = ((dividend / 3) + dividend) & 3 , reduciendo el módulo a una suma más una lógica ya que ya tenemos el resultado de la división. Esta es una instancia del algoritmo general: dividendo % (2 n -1) = ((dividendo / (2 n -1) + dividendo) & (2 n -1).
Finalmente para la división por 3 en el término de corrección (a % 3 + b % 3 + c % 3) / 3 no necesitamos el código para la división genérica por 3. Como el dividendo es muy pequeño, en [0, 6 ], podemos simplificar x / 3 en (3 * x) / 8 , lo que requiere solo un scale_add más un cambio .
El siguiente código muestra mi trabajo actual en progreso. El uso de Compiler Explorer para verificar el código generado para varias plataformas muestra el código ajustado que esperaría (cuando se compila con -O3 ).
Sin embargo, al cronometrar el código en mi máquina Ivy Bridge x86_64 usando el compilador Intel 13.x, se hizo evidente una falla: mientras que mi código mejora la latencia (de 18 ciclos a 15 ciclos para datos uint64_t ) en comparación con la versión simple, el rendimiento empeora ( de un resultado cada 6,8 ciclos a un resultado cada 8,5 ciclos para datos uint64_t ). Mirando el código ensamblador más de cerca, es bastante evidente por qué es así: básicamente logré reducir el código de un paralelismo de tres vías a un paralelismo de dos vías.
¿Existe una técnica de optimización de aplicación genérica, beneficiosa para los procesadores comunes, en particular todos los tipos de x86 y ARM, así como para las GPU, que conserve más paralelismo? Alternativamente, ¿existe una técnica de optimización que reduzca aún más el recuento total de operaciones para compensar el paralelismo reducido? El cálculo del término de corrección ( tail en el código a continuación) parece un buen objetivo. La simplificación (carry_mod_3 + sum_mod_3) / 2 parecía tentadora pero ofrece un resultado incorrecto para una de las nueve combinaciones posibles.
#include <stdio.h> #include <stdlib.h> #include <stdint.h> #define BENCHMARK (1) #define SIMPLE_COMPUTATION (0) #if BENCHMARK #define T uint64_t #else // !BENCHMARK #define T uint8_t #endif // BENCHMARK T average_of_3 (T a, T b, T c) { T avg; #if SIMPLE_COMPUTATION avg = a / 3 + b / 3 + c / 3 + (a % 3 + b % 3 + c % 3) / 3; #else // !SIMPLE_COMPUTATION /* carry save adder */ T a_xor_b = a ^ b; T sum = a_xor_b ^ c; T carry = (a_xor_b & c) | (a & b); /* here 2 * carry + sum = a + b + c */ T sum_div_3 = (sum / 3); // {MUL|MULHI}, SHR T sum_mod_3 = (sum + sum_div_3) & 3; // ADD, AND if (sizeof (size_t) == sizeof (T)) { // "native precision" (well, not always) T two_carry_div_3 = (carry / 3) * 2; // MULHI, ANDN T two_carry_mod_3 = (2 * carry + two_carry_div_3) & 6; // SCALE_ADD, AND T head = two_carry_div_3 + sum_div_3; // ADD T tail = (3 * (two_carry_mod_3 + sum_mod_3)) / 8; // ADD, SCALE_ADD, SHR avg = head + tail; // ADD } else { T carry_div_3 = (carry / 3); // MUL, SHR T carry_mod_3 = (carry + carry_div_3) & 3; // ADD, AND T head = (2 * carry_div_3 + sum_div_3); // SCALE_ADD T tail = (3 * (2 * carry_mod_3 + sum_mod_3)) / 8; // SCALE_ADD, SCALE_ADD, SHR avg = head + tail; // ADD } #endif // SIMPLE_COMPUTATION return avg; } #if !BENCHMARK /* Test correctness on 8-bit data exhaustively. Should catch most errors */ int main (void) { T a, b, c, res, ref; a = 0; do { b = 0; do { c = 0; do { res = average_of_3 (a, b, c); ref = ((uint64_t)a + (uint64_t)b + (uint64_t)c) / 3; if (res != ref) { printf ("a=%08x b=%08x c=%08x res=%08x ref=%08x\n", a, b, c, res, ref); return EXIT_FAILURE; } c++; } while (c); b++; } while (b); a++; } while (a); return EXIT_SUCCESS; } #else // BENCHMARK #include <math.h> // A routine to give access to a high precision timer on most systems. #if defined(_WIN32) #if !defined(WIN32_LEAN_AND_MEAN) #define WIN32_LEAN_AND_MEAN #endif #include <windows.h> double second (void) { LARGE_INTEGER t; static double oofreq; static int checkedForHighResTimer; static BOOL hasHighResTimer; if (!checkedForHighResTimer) { hasHighResTimer = QueryPerformanceFrequency (&t); oofreq = 1.0 / (double)t.QuadPart; checkedForHighResTimer = 1; } if (hasHighResTimer) { QueryPerformanceCounter (&t); return (double)t.QuadPart * oofreq; } else { return (double)GetTickCount() * 1.0e-3; } } #elif defined(__linux__) || defined(__APPLE__) #include <stddef.h> #include <sys/time.h> double second (void) { struct timeval tv; gettimeofday(&tv, NULL); return (double)tv.tv_sec + (double)tv.tv_usec * 1.0e-6; } #else #error unsupported platform #endif #define N (3000000) int main (void) { double start, stop, elapsed = INFINITY; int i, k; T a, b; T avg0 = 0xffffffff, avg1 = 0xfffffffe; T avg2 = 0xfffffffd, avg3 = 0xfffffffc; T avg4 = 0xfffffffb, avg5 = 0xfffffffa; T avg6 = 0xfffffff9, avg7 = 0xfffffff8; T avg8 = 0xfffffff7, avg9 = 0xfffffff6; T avg10 = 0xfffffff5, avg11 = 0xfffffff4; T avg12 = 0xfffffff2, avg13 = 0xfffffff2; T avg14 = 0xfffffff1, avg15 = 0xfffffff0; a = 0x31415926; b = 0x27182818; avg0 = average_of_3 (a, b, avg0); for (k = 0; k < 5; k++) { start = second(); for (i = 0; i < N; i++) { avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); avg0 = average_of_3 (a, b, avg0); b = (b + avg0) ^ a; a = (a ^ b) + avg0; } stop = second(); elapsed = fmin (stop - start, elapsed); } printf ("a=%016llx b=%016llx avg=%016llx", (uint64_t)a, (uint64_t)b, (uint64_t)avg0); printf ("\rlatency: each average_of_3() took %.6e seconds\n", elapsed / 16 / N); a = 0x31415926; b = 0x27182818; avg0 = average_of_3 (a, b, avg0); for (k = 0; k < 5; k++) { start = second(); for (i = 0; i < N; i++) { avg0 = average_of_3 (a, b, avg0); avg1 = average_of_3 (a, b, avg1); avg2 = average_of_3 (a, b, avg2); avg3 = average_of_3 (a, b, avg3); avg4 = average_of_3 (a, b, avg4); avg5 = average_of_3 (a, b, avg5); avg6 = average_of_3 (a, b, avg6); avg7 = average_of_3 (a, b, avg7); avg8 = average_of_3 (a, b, avg8); avg9 = average_of_3 (a, b, avg9); avg10 = average_of_3 (a, b, avg10); avg11 = average_of_3 (a, b, avg11); avg12 = average_of_3 (a, b, avg12); avg13 = average_of_3 (a, b, avg13); avg14 = average_of_3 (a, b, avg14); avg15 = average_of_3 (a, b, avg15); b = (b + avg0) ^ a; a = (a ^ b) + avg0; } stop = second(); elapsed = fmin (stop - start, elapsed); } printf ("a=%016llx b=%016llx avg=%016llx", (uint64_t)a, (uint64_t)b, (uint64_t)(avg0 + avg1 + avg2 + avg3 + avg4 + avg5 + avg6 + avg7 + avg8 + avg9 +avg10 +avg11 +avg12 +avg13 +avg14 +avg15)); printf ("\rthroughput: each average_of_3() took %.6e seconds\n", elapsed / 16 / N); return EXIT_SUCCESS; } #endif // BENCHMARKSospecho que SIMPLE está superando el punto de referencia de rendimiento al hacer CSE y sacar a/3+b/3 y a%3+b%3 del ciclo, reutilizando esos resultados para los 16 resultados avg0..15 .
(La versión SIMPLE puede levantar mucho más trabajo que la versión complicada; realmente solo a ^ b y a & b en esa versión).
Obligar a la función a que no esté en línea introduce más sobrecarga de front-end, pero hace que su versión gane, como esperamos que sea en una CPU con búferes de ejecución fuera de orden profundos para superponer el trabajo independiente. Hay muchos ILP para encontrar en las iteraciones, para el punto de referencia de rendimiento. (No miré de cerca el asm para la versión no en línea).
https://godbolt.org/z/j95qn3 (usando __attribute__((noinline)) con clang -O3 -march=skylake en las CPU SKX de Godbolt) muestra un rendimiento de 2,58 nanosegundos para la forma simple, un rendimiento de 2,48 nanosegundos para su forma. frente a un rendimiento de 1,17 nanosegundos con inserción para la versión simple.
-march=skylake permite mulx para una multiplicación completa más flexible, pero por lo demás no se beneficia de BMI2. andn no se usa; la línea que comentó con mulhi / andn es mulx en RCX / and rcx, -2 que solo requiere un signo inmediato extendido.
Otra forma de hacer esto sin forzar la sobrecarga de call/ret sería asm en línea como en Prevención de optimizaciones del compilador durante la evaluación comparativa (la charla CppCon de Chandler Carruth tiene un ejemplo de cómo usa un par de contenedores), o el punto de referencia de Google Benchmark benchmark::DoNotOptimize .
Específicamente, GNU C asm("" : "+r"(a), "+r"(b)) entre cada avgX = average_of_3 (a, b, avgX); hará que el compilador olvide todo lo que sabe sobre los valores de a y b , mientras los mantiene en registros.
Mi respuesta sobre No entiendo la definición de DoNotOptimizeAway entra en más detalles sobre el uso de una restricción de registro "r" de solo lectura para obligar al compilador a materializar un resultado en un registro, frente a "+r" para que asuma el el valor ha sido modificado.
Si entiende bien el asm en línea GNU C, puede ser más fácil implementar el suyo propio de manera que sepa exactamente lo que hacen.
No estoy seguro de si se ajusta a sus requisitos, pero tal vez funcione simplemente para calcular el resultado y luego corregir el error del desbordamiento:
T average_of_3 (T a, T b, T c) { T r = ((T) (a + b + c)) / 3; T o = (a > (T) ~b) + ((T) (a + b) > (T) (~c)); if (o) r += ((T) 0x5555555555555555) << (o - 1); T rem = ((T) (a + b + c)) % 3; if (rem >= (3 - o)) ++r; return r; } [EDITAR] Esta es la mejor versión sin ramas y sin comparación que se me ocurre. En mi máquina, esta versión en realidad tiene un rendimiento ligeramente mayor que el código de njuffa. __builtin_add_overflow(x, y, r) es compatible con gcc y clang y devuelve 1 si la suma x + y desborda el tipo de *r y 0 de lo contrario, por lo que el cálculo de o es equivalente al código portátil en la primera versión, pero al menos gcc produce mejor código con el incorporado.
T average_of_3 (T a, T b, T c) { T r = ((T) (a + b + c)) / 3; T rem = ((T) (a + b + c)) % 3; T dummy; T o = __builtin_add_overflow(a, b, &dummy) + __builtin_add_overflow((T) (a + b), c, &dummy); r += -((o - 1) & 0xaaaaaaaaaaaaaaab) ^ 0x5555555555555555; r += (rem + o + 1) >> 2; return r; }[Falk Hüffner señala en los comentarios que esta respuesta tiene similitudes con su respuesta . Mirando su código más de cerca con retraso, encuentro algunas similitudes. Sin embargo, lo que publiqué aquí es producto de un proceso de pensamiento independiente, una continuación de mi idea original "reducir tres elementos a dos antes de div-mod". Entendí que el enfoque de Hüffner era diferente: "cálculo ingenuo seguido de correcciones".]
He encontrado una mejor manera que la técnica CSA en mi pregunta para reducir el trabajo de división y módulo de tres operandos a dos operandos. Primero, forme la suma completa de dos palabras, luego aplique la división y el módulo por 3 a cada una de las mitades por separado, finalmente combine los resultados. Dado que la mitad más significativa solo puede tomar los valores 0, 1 o 2, calcular el cociente y el resto de la división por tres es trivial. Además, la combinación en el resultado final se vuelve más simple.
En comparación con la variante de código no simple de la pregunta, esto logra una aceleración en todas las plataformas que examiné. La calidad del código generado por los compiladores para la adición simulada de palabras dobles varía, pero en general es satisfactoria. No obstante, puede valer la pena codificar esta parte de forma no portátil, por ejemplo, con ensamblado en línea.
T average_of_3_hilo (T a, T b, T c) { const T fives = (((T)(~(T)0)) / 3); // 0x5555... T avg, hi, lo, lo_div_3, lo_mod_3, hi_div_3, hi_mod_3; /* compute the full sum a + b + c into the operand pair hi:lo */ lo = a + b; hi = lo < a; lo = c + lo; hi = hi + (lo < c); /* determine quotient and remainder of each half separately */ lo_div_3 = lo / 3; lo_mod_3 = (lo + lo_div_3) & 3; hi_div_3 = hi * fives; hi_mod_3 = hi; /* combine partial results into the division result for the full sum */ avg = lo_div_3 + hi_div_3 + ((lo_mod_3 + hi_mod_3 + 1) / 4); return avg; }Déjame tirar mi sombrero en el ring. No estoy haciendo nada demasiado complicado aquí, creo.
#include <stdint.h> uint64_t average_of_three(uint64_t a, uint64_t b, uint64_t c) { uint64_t hi = (a >> 32) + (b >> 32) + (c >> 32); uint64_t lo = hi + (a & 0xffffffff) + (b & 0xffffffff) + (c & 0xffffffff); return 0x55555555 * hi + lo / 3; }Siguiendo la discusión a continuación sobre las diferentes divisiones, aquí hay una versión que ahorra una multiplicación a expensas de tres AND bit a bit:
T hi = (a >> 2) + (b >> 2) + (c >> 2); T lo = (a & 3) + (b & 3) + (c & 3); avg = hi + (hi + lo) / 3;Nueva respuesta, nueva idea. Este está basado en la identidad matemática.
floor((a+b+c)/3) = floor(x + (a+b+c - 3x)/3) ¿Cuándo funciona esto con enteros de máquina y división sin signo?
Cuando la diferencia no se ajusta, es decir, 0 ≤ a+b+c - 3x ≤ T_MAX .
Esta definición de x es rápida y hace el trabajo.
T avg3(T a, T b, T c) { T x = (a >> 2) + (b >> 2) + (c >> 2); return x + (a + b + c - 3 * x) / 3; }Extrañamente, ICC inserta un neg adicional a menos que haga esto:
T avg3(T a, T b, T c) { T x = (a >> 2) + (b >> 2) + (c >> 2); return x + (a + b + c - (x + x * 2)) / 3; } Tenga en cuenta que T debe tener al menos cinco bits de ancho.
Si T tiene dos palabras de plataforma, puede ahorrar algunas operaciones de palabra doble omitiendo la palabra baja de x .
¿Versión alternativa con peor latencia pero quizás un rendimiento ligeramente mayor?
T lo = a + b; T hi = lo < b; lo += c; hi += lo < c; T x = (hi << (sizeof(T) * CHAR_BIT - 2)) + (lo >> 2); avg = x + (T)(lo - 3 * x) / 3;Ya respondí la pregunta a la que se vinculó, así que solo respondo la parte que es diferente de esta: rendimiento.
Si realmente te importaba el rendimiento, entonces la respuesta es:
( a + b + c ) / 3Dado que le importaba el rendimiento, debe tener una intuición sobre el tamaño de los datos con los que está trabajando. No debería haberse preocupado por el desbordamiento en la suma (la multiplicación es otra cuestión) de solo 3 valores, porque si sus datos ya son lo suficientemente grandes como para usar los bits altos del tipo de datos elegido, corre peligro de desbordamiento de todos modos y debería haber usado un tipo entero más grande. Si está desbordado en uint64_t, entonces realmente debería preguntarse por qué exactamente necesita contar con precisión hasta 18 quintillones, y tal vez considerar usar float o double.
Ahora, habiendo dicho todo eso, les daré mi respuesta real: No importa. La pregunta no surge en la vida real y cuando lo hace, el rendimiento no importa.
Podría ser una pregunta de rendimiento real si lo está haciendo un millón de veces en SIMD, porque allí está realmente incentivado para usar números enteros de menor ancho y es posible que necesite ese último espacio libre, pero esa no era su pregunta.
Una compilación experimental de GCC-11 compila la función ingenua obvia en algo como:
uint32_t avg3t (uint32_t a, uint32_t b, uint32_t c) { a += b; b = a < b; a += c; b += a < c; b = b + a; b += b < a; return (a - (b % 3)) * 0xaaaaaaab; }Lo cual es similar a algunas de las otras respuestas publicadas aquí. Cualquier explicación de cómo funcionan estas soluciones sería bienvenida (no estoy seguro de la etiqueta aquí).