algorithm - usual - La raíz cuadrada de enteros más rápida

tutorial de raiz cuadrada (5)

Estoy en la necesidad de una raíz cuadrada Integer rápida que no implique ninguna división explícita. La arquitectura de RISC de destino puede realizar operaciones como agregar, mul, sub, cambiar en un ciclo (bueno, el resultado de la operación está escrito en el tercer ciclo, en realidad, pero hay intercalado), por lo que cualquier algoritmo Integer que use estas operaciones sea rápido apreciado.

Esto es lo que tengo ahora y estoy pensando que una búsqueda binaria debería ser más rápida, ya que el siguiente bucle se ejecuta 16 veces cada vez (independientemente del valor). Todavía no lo he depurado extensivamente (pero pronto), así que quizás sea posible tener una salida temprana allí:

unsigned short int int_sqrt32(unsigned int x) { unsigned short int res=0; unsigned short int add= 0x8000; int i; for(i=0;i<16;i++) { unsigned short int temp=res | add; unsigned int g2=temp*temp; if (x>=g2) { res=temp; } add>>=1; } return res; }

Parece que el costo de rendimiento actual de lo anterior [en el contexto del RISC objetivo] es un bucle de 5 instrucciones (conjunto de bits, mul, comparar, almacenar, cambiar). Probablemente no haya espacio en el caché para desenrollarse completamente (pero este será el candidato principal para un desenrollamiento parcial [por ejemplo, un bucle de 4 en lugar de 16], por supuesto). Por lo tanto, el costo es de 16 * 5 = 80 instrucciones (más la sobrecarga del bucle, si no se desenrolla). Lo que, si está totalmente intercalado, costaría solo 80 ciclos (+2 para la última instrucción).

¿Puedo obtener alguna otra implementación de sqrt (usando solo add, mul, bitshift, store / cmp) en 82 ciclos?

PREGUNTAS MÁS FRECUENTES:

¿Por qué no confías en el compilador para producir un buen código rápido?

No hay ningún compilador C-> RISC en funcionamiento para la plataforma. Transmitiré el código de referencia C actual al ASM RISC escrito a mano.

¿Perfilaste el código para ver si sqrt es realmente un cuello de botella?

No, no hay necesidad de eso. El chip RISC objetivo es de unos veinte MHz, por lo que cada instrucción cuenta. El bucle central (que calcula el factor de forma de transferencia de energía entre el tirador y el parche del receptor), donde se usa este sqrt , se ejecutará ~ 1,000 veces cada fotograma de renderización (suponiendo que será lo suficientemente rápido, por supuesto), hasta 60,000 por segundo , y aproximadamente 1,000,000 de veces para una demo completa.

¿Ha intentado optimizar el algoritmo para eliminar tal vez el sqrt ?

Sí, ya lo hice. De hecho, me deshice de 2 sqrt s y muchas divisiones (eliminadas o reemplazadas por turnos). Puedo ver un gran aumento de rendimiento (en comparación con la versión flotante de referencia) incluso en mi portátil gigahertz.

¿Cuál es la aplicación?

Es un procesador de radiosidad de refinamiento progresivo en tiempo real para la demostración compo. La idea es tener un ciclo de disparo en cada fotograma, de modo que visiblemente converja y se vea mejor con cada fotograma renderizado (por ejemplo, hasta 60 veces por segundo, aunque el rasterizador SW probablemente no sea tan rápido [pero al menos puede ejecutarse en el otro chip en paralelo con el RISC, por lo que si toma 2-3 cuadros para representar la escena, el RISC habrá trabajado en 2-3 cuadros de datos de radiosidad, en paralelo]).

¿Por qué no trabajas directamente en ASM objetivo?

Debido a que la radiosidad es un algoritmo ligeramente involucrado, necesito la capacidad de edición y continuación de depuración instantánea de Visual Studio. Lo que he hecho durante el fin de semana en VS (un par de cientos de cambios de código para convertir las matemáticas de punto flotante a solo enteros) me llevaría 6 meses en la plataforma de destino con solo imprimir "depuración".

¿Por qué no puedes usar una división?

Porque es 16 veces más lento en el RISC objetivo que cualquiera de los siguientes: mul, add, sub, shift, compare, load / store (que toma solo 1 ciclo). Por lo tanto, se usa solo cuando es absolutamente necesario (un par de veces, desafortunadamente, cuando no se pudo usar el cambio).

¿Se pueden usar tablas de consulta?

El motor ya necesita otros LUT y la copia desde la RAM principal al pequeño caché de RISC es prohibitivamente costosa (y definitivamente no todos y cada uno de los cuadros). Pero, tal vez podría ahorrar 128-256 Bytes si me diera al menos un aumento de 100-200% por sqrt.

¿Cuál es el rango de los valores para sqrt ?

Logré reducirlo a mero int de 32 bits sin firmar (4,294,967,295)

Desde el camino de comentarios, parece que el procesador RISC solo proporciona 32x32-> 32 bits de multiplicación y 16x16-> 32 bits de multiplicación. No se proporciona una multiplicación de ampliación de 32x-32-> 64 bits, o una instrucción MULHI que devuelve los 32 bits superiores de un producto de 64 bits.

Esto parecería excluir los enfoques basados en la iteración de Newton-Raphson, que probablemente serían ineficientes, ya que generalmente requieren una instrucción MULHI o una multiplicación de la ampliación para la aritmética de punto fijo intermedio.

El código C99 a continuación utiliza un enfoque iterativo diferente que requiere solo 16x16-> 32 bits multiplicados, pero converge de forma un tanto lineal, requiriendo hasta seis iteraciones. Este enfoque requiere la funcionalidad CLZ para determinar rápidamente una estimación inicial para las iteraciones. Asker declaró en los comentarios que el procesador RISC utilizado no proporciona la funcionalidad CLZ. Por lo tanto, se requiere la emulación de CLZ, y dado que la emulación se suma al almacenamiento y al conteo de instrucciones, esto puede hacer que este enfoque no sea competitivo. Realicé una búsqueda de fuerza bruta para determinar la tabla de búsqueda deBruijn con el multiplicador más pequeño.

Este algoritmo iterativo entrega resultados brutos bastante cercanos a los resultados deseados, es decir (int)sqrt(x) , pero siempre algo en el lado alto debido a la naturaleza truncadora de la aritmética de enteros. Para llegar al resultado final, el resultado se decrementa condicionalmente hasta que el cuadrado del resultado es menor o igual que el argumento original.

El uso del calificador volatile en el código solo sirve para establecer que todas las variables nombradas pueden asignarse como datos de 16 bits sin afectar la funcionalidad. No sé si esto proporciona alguna ventaja, pero noté que el OP utilizó específicamente variables de 16 bits en su código.

Tenga en cuenta que para la mayoría de los procesadores, los pasos de corrección al final no deben incluir ninguna ramificación. El producto y*y se puede restar de x con carry-out (o préstamo), luego y se corrige mediante una resta con carry-in (o préstamo). Así que cada paso debe ser una secuencia MUL , SUBcc , SUBcc .

Debido a que la implementación de la iteración por un bucle incurre en una sobrecarga considerable, he elegido desenrollar completamente el bucle, pero proporcionar dos comprobaciones anticipadas. Al contar las operaciones manualmente, cuento 46 operaciones para el caso más rápido, 54 operaciones para el caso promedio y 60 operaciones para el caso más desfavorable.

#include <stdio.h> #include <stdlib.h> #include <stdint.h> #include <math.h> static const uint8_t clz_tab[32] = { 31, 22, 30, 21, 18, 10, 29, 2, 20, 17, 15, 13, 9, 6, 28, 1, 23, 19, 11, 3, 16, 14, 7, 24, 12, 4, 8, 25, 5, 26, 27, 0}; uint8_t clz (uint32_t a) { a |= a >> 16; a |= a >> 8; a |= a >> 4; a |= a >> 2; a |= a >> 1; return clz_tab [0x07c4acdd * a >> 27]; } /* 16 x 16 -> 32 bit unsigned multiplication; should be single instruction */ uint32_t umul16w (uint16_t a, uint16_t b) { return (uint32_t)a * b; } /* Reza Hashemian, "Square Rooting Algorithms for Integer and Floating-Point Numbers", IEEE Transactions on Computers, Vol. 39, No. 8, Aug. 1990, p. 1025 */ uint16_t isqrt (uint32_t x) { volatile uint16_t y, z, lsb, mpo, mmo, lz, t; if (x == 0) return x; // early out, code below can''t handle zero lz = clz (x); // #leading zeros, 32-lz = #bits of argument lsb = lz & 1; mpo = 17 - (lz >> 1); // m+1, result has roughly half the #bits of argument mmo = mpo - 2; // m-1 t = 1 << mmo; // power of two for two''s complement of initial guess y = t - (x >> (mpo - lsb)); // initial guess for sqrt t = t + t; // power of two for two''s complement of result z = y; y = (umul16w (y, y) >> mpo) + z; y = (umul16w (y, y) >> mpo) + z; if (x >= 0x40400) { y = (umul16w (y, y) >> mpo) + z; y = (umul16w (y, y) >> mpo) + z; if (x >= 0x1002000) { y = (umul16w (y, y) >> mpo) + z; y = (umul16w (y, y) >> mpo) + z; } } y = t - y; // raw result is 2''s complement of iterated solution y = y - umul16w (lsb, (umul16w (y, 19195) >> 16)); // mult. by sqrt(0.5) if ((int32_t)(x - umul16w (y, y)) < 0) y--; // iteration may overestimate if ((int32_t)(x - umul16w (y, y)) < 0) y--; // result, adjust downward if if ((int32_t)(x - umul16w (y, y)) < 0) y--; // necessary return y; // (int)sqrt(x) } int main (void) { uint32_t x = 0; uint16_t res, ref; do { ref = (uint16_t)sqrt((double)x); res = isqrt (x); if (res != ref) { printf ("!!!! x=%08x res=%08x ref=%08x/n", x, res, ref); return EXIT_FAILURE; } x++; } while (x); return EXIT_SUCCESS; }

Otra posibilidad es usar la iteración de Newton para la raíz cuadrada, a pesar del alto costo de la división. Para entradas pequeñas solo se requerirá una iteración. Aunque el autor de la pregunta no lo dijo, en base al tiempo de ejecución de 16 ciclos para la operación DIV , sospecho que esta es realmente una división de 32/16->16 bits que requiere un código de guarda adicional para evitar el desbordamiento, definido como un cociente que No cabe en 16 bits. He añadido protecciones apropiadas a mi código en base a esta suposición.

Dado que la iteración de Newton duplica el número de bits buenos cada vez que se aplica, solo necesitamos una estimación inicial de baja precisión que se puede recuperar fácilmente de una tabla basada en los cinco bits principales del argumento. Para capturarlos, primero normalizamos el argumento en un formato de punto fijo de 2.30 con un factor de escala implícito adicional de 2 ^{32- (lz y ~ 1)} donde lz es el número de ceros iniciales en el argumento. Como en el enfoque anterior, la iteración no siempre ofrece un resultado preciso, por lo que se debe aplicar una corrección si el resultado preliminar es demasiado grande. Cuento 49 ciclos para la ruta rápida, 70 ciclos para la ruta lenta (promedio de 60 ciclos).

static const uint16_t sqrt_tab[32] = { 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x85ff, 0x8cff, 0x94ff, 0x9aff, 0xa1ff, 0xa7ff, 0xadff, 0xb3ff, 0xb9ff, 0xbeff, 0xc4ff, 0xc9ff, 0xceff, 0xd3ff, 0xd8ff, 0xdcff, 0xe1ff, 0xe6ff, 0xeaff, 0xeeff, 0xf3ff, 0xf7ff, 0xfbff, 0xffff }; /* 32/16->16 bit division. Note: Will overflow if x[31:16] >= y */ uint16_t udiv_32_16 (uint32_t x, uint16_t y) { uint16_t r = x / y; return r; } /* table lookup for initial guess followed by division-based Newton iteration*/ uint16_t isqrt (uint32_t x) { volatile uint16_t q, lz, y, i, xh; if (x == 0) return x; // early out, code below can''t handle zero // initial guess based on leading 5 bits of argument normalized to 2.30 lz = clz (x); i = ((x << (lz & ~1)) >> 27); y = sqrt_tab[i] >> (lz >> 1); xh = (x >> 16); // needed for overflow check on division // first Newton iteration, guard against overflow in division q = 0xffff; if (xh < y) q = udiv_32_16 (x, y); y = (q + y) >> 1; if (lz < 10) { // second Newton iteration, guard against overflow in division q = 0xffff; if (xh < y) q = udiv_32_16 (x, y); y = (q + y) >> 1; } if (umul16w (y, y) > x) y--; // adjust quotient if too large return y; // (int)sqrt(x) }

Echa un vistazo here .

Por ejemplo, en 3 (a) existe este método, que es trivialmente adaptable para hacer una raíz cuadrada de 64-> 32 bits, y también trivialmente transcribible al ensamblador:

/* by Jim Ulery */ static unsigned julery_isqrt(unsigned long val) { unsigned long temp, g=0, b = 0x8000, bshft = 15; do { if (val >= (temp = (((g << 1) + b)<<bshft--))) { g += b; val -= temp; } } while (b >>= 1); return g; }

Sin divisiones, sin multiplicaciones, solo cambios de bit. Sin embargo, el tiempo empleado será algo impredecible, especialmente si utiliza una rama (en ARM RISC, las instrucciones condicionales funcionarían).

En general, esta página enumera formas de calcular raíces cuadradas. Si desea producir una raíz cuadrada inversa rápida (es decir, x**(-0.5) ), o simplemente está interesado en formas asombrosas de optimizar el código, eche un vistazo a this , this y en.wikipedia.org/wiki/… .

Esto es lo mismo que el tuyo, pero con menos operaciones. (Cuento 9 operaciones en el bucle en su código, incluyendo prueba e incremento i en el bucle for y 3 asignaciones, pero ¿quizás algunas de ellas desaparecen cuando se codifican en ASM? Hay 6 operaciones en el código a continuación, si cuenta g*g>n como dos (sin asignación)).

int isqrt(int n) { int g = 0x8000; int c = 0x8000; for (;;) { if (g*g > n) { g ^= c; } c >>= 1; if (c == 0) { return g; } g |= c; } }

Lo tengo here . Es posible que pueda eliminar una comparación si desenrolla el bucle y salta al lugar adecuado en función del bit más alto que no sea cero en la entrada.

Actualizar

He estado pensando más en usar el método de Newton. En teoría, el número de bits de precisión debería duplicarse para cada iteración. Eso significa que el método de Newton es mucho peor que cualquiera de las otras sugerencias cuando hay pocos bits correctos en la respuesta; sin embargo, la situación cambia cuando hay muchos bits correctos en la respuesta. Teniendo en cuenta que la mayoría de las sugerencias parecen tomar 4 ciclos por bit, eso significa que una iteración del método de Newton (16 ciclos para la división + 1 para la adición + 1 para el cambio = 18 ciclos) no vale la pena a menos que proporcione más de 4 bits.

Entonces, mi sugerencia es acumular 8 bits de la respuesta por uno de los métodos sugeridos (8 * 4 = 32 ciclos) y luego realizar una iteración del método de Newton (18 ciclos) para duplicar el número de bits a 16. Eso es un total de 50 ciclos (más quizás un extra de 4 ciclos para obtener 9 bits antes de aplicar el método de Newton, y quizás 2 ciclos para superar el rebasamiento ocasionalmente experimentado por el método de Newton). Eso es un máximo de 56 ciclos que, por lo que puedo ver, rivaliza con cualquiera de las otras sugerencias.

Segunda actualización

Codifiqué la idea del algoritmo híbrido. El método de Newton en sí no tiene sobrecarga; solo aplica y duplica el número de dígitos significativos. El problema es tener un número predecible de dígitos significativos antes de aplicar el método de Newton. Para eso, necesitamos averiguar dónde aparecerá la parte más significativa de la respuesta. Usando una modificación del método de secuencia DeBruijn rápido dado por otro póster, podemos realizar ese cálculo en aproximadamente 12 ciclos en mi estimación. Por otro lado, saber la posición del msb de la respuesta acelera todos los métodos (promedio, no en el peor de los casos), por lo que parece que vale la pena sin importar qué.

Después de calcular el msb de la respuesta, ejecuto varias rondas del algoritmo sugerido anteriormente, luego lo termino con una o dos rondas del método de Newton. Decidimos cuándo ejecutar el método de Newton mediante el siguiente cálculo: un bit de la respuesta toma aproximadamente 8 ciclos según el cálculo de los comentarios; una ronda del método de Newton toma alrededor de 18 ciclos (división, adición y desplazamiento, y quizás asignación), por lo que solo deberíamos ejecutar el método de Newton si vamos a obtener al menos tres bits de él. Entonces, para respuestas de 6 bits, podemos ejecutar el método lineal 3 veces para obtener 3 bits significativos, luego ejecutar el método de Newton 1 vez para obtener otro 3. Para respuestas de 15 bits, ejecutamos el método lineal 4 veces para obtener 4 bits, luego el de Newton Método dos veces para obtener otros 4 y luego otros 7. Y así sucesivamente.

Esos cálculos dependen de saber exactamente cuántos ciclos se requieren para obtener un poco por el método lineal en comparación con cuántos son requeridos por el método de Newton. Si la "economía" cambia, por ejemplo, al descubrir una forma más rápida de acumular bits de forma lineal, la decisión de cuándo invocar el método de Newton cambiará.

Desenrollé los bucles e implementé las decisiones como conmutadores, que espero se traduzcan en búsquedas rápidas de tablas en ensamblador. No estoy absolutamente seguro de tener el número mínimo de ciclos en cada caso, por lo que tal vez sea posible un ajuste adicional. Por ejemplo, para s = 10, puede intentar obtener 5 bits y luego aplicar el método de Newton una vez en lugar de dos veces.

He probado el algoritmo a fondo para la corrección. Algunas aceleraciones menores adicionales son posibles si está dispuesto a aceptar respuestas ligeramente incorrectas en algunos casos. Se usan al menos dos ciclos después de aplicar el método de Newton para corregir un error off-by-one que ocurre con números de la forma m^2-1 . Y se usa un ciclo para probar la entrada 0 al principio, ya que el algoritmo no puede manejar esa entrada. Si sabes que nunca vas a tomar la raíz cuadrada de cero, puedes eliminar esa prueba. Finalmente, si solo necesita 8 bits significativos en la respuesta, puede descartar uno de los cálculos del método de Newton.

#include <inttypes.h> #include <stdint.h> #include <stdbool.h> #include <stdio.h> uint32_t isqrt1(uint32_t n); int main() { uint32_t n; bool it_works = true; for (n = 0; n < UINT32_MAX; ++n) { uint32_t sr = isqrt1(n); if ( sr*sr > n || ( sr < 65535 && (sr+1)*(sr+1) <= n )) { it_works = false; printf("isqrt(%" PRIu32 ") = %" PRIu32 "/n", n, sr); } } if (it_works) { printf("it works/n"); } return 0; } /* table modified to return shift s to move 1 to msb of square root of x */ /* static const uint8_t debruijn32[32] = { 0, 31, 9, 30, 3, 8, 13, 29, 2, 5, 7, 21, 12, 24, 28, 19, 1, 10, 4, 14, 6, 22, 25, 20, 11, 15, 23, 26, 16, 27, 17, 18 }; */ static const uint8_t debruijn32[32] = { 15, 0, 11, 0, 14, 11, 9, 1, 14, 13, 12, 5, 9, 3, 1, 6, 15, 10, 13, 8, 12, 4, 3, 5, 10, 8, 4, 2, 7, 2, 7, 6 }; /* based on CLZ emulation for non-zero arguments, from * http://.com/questions/23856596/counting-leading-zeros-in-a-32-bit-unsigned-integer-with-best-algorithm-in-c-pro */ uint8_t shift_for_msb_of_sqrt(uint32_t x) { x |= x >> 1; x |= x >> 2; x |= x >> 4; x |= x >> 8; x |= x >> 16; x++; return debruijn32 [x * 0x076be629 >> 27]; } uint32_t isqrt1(uint32_t n) { if (n==0) return 0; uint32_t s = shift_for_msb_of_sqrt(n); uint32_t c = 1 << s; uint32_t g = c; switch (s) { case 9: case 5: if (g*g > n) { g ^= c; } c >>= 1; g |= c; case 15: case 14: case 13: case 8: case 7: case 4: if (g*g > n) { g ^= c; } c >>= 1; g |= c; case 12: case 11: case 10: case 6: case 3: if (g*g > n) { g ^= c; } c >>= 1; g |= c; case 2: if (g*g > n) { g ^= c; } c >>= 1; g |= c; case 1: if (g*g > n) { g ^= c; } c >>= 1; g |= c; case 0: if (g*g > n) { g ^= c; } } /* now apply one or two rounds of Newton''s method */ switch (s) { case 15: case 14: case 13: case 12: case 11: case 10: g = (g + n/g) >> 1; case 9: case 8: case 7: case 6: g = (g + n/g) >> 1; } /* correct potential error at m^2-1 for Newton''s method */ return (g==65536 || g*g>n) ? g-1 : g; }

En las pruebas de luz en mi máquina (que no se parecen en nada a las suyas), la nueva rutina isqrt1 ejecuta en promedio un 40% más rápido que la rutina isqrt anterior que proporcioné.

No sé cómo convertirlo en un algoritmo eficiente, pero cuando investigué esto en los años 80 surgió un patrón interesante. Al redondear las raíces cuadradas, hay dos enteros más con esa raíz cuadrada que la anterior (después de cero).

Entonces, un número (cero) tiene una raíz cuadrada de cero, dos tiene una raíz cuadrada de 1 (1 y 2), 4 tiene una raíz cuadrada de dos (3, 4, 5 y 6) y así sucesivamente. Probablemente no sea una respuesta útil pero interesante no obstante.

Si la multiplicación es la misma velocidad (o más rápida que!) De suma y desplazamiento, o si no cuenta con una instrucción rápida de cambio por cantidad contenida en un registro, lo siguiente no será útil. De otra manera:

Estás calculando temp*temp nuevo en cada ciclo de bucle, pero temp = res | add temp = res | add , que es lo mismo que res + add ya que sus bits no se superponen, y (a) ya ha calculado las res*res en un ciclo de bucle anterior, y (b) add tiene una estructura especial (siempre es un solo bit) ). Entonces, aprovechando el hecho de que (a+b)^2 = a^2 + 2ab + b^2 , y ya tiene a^2 , y b^2 es un solo bit desplazado dos veces más hacia la izquierda que el bit único b , y 2ab es solo a desplazada a la izquierda en 1 posición más que la ubicación del bit único en b , puede deshacerse de la multiplicación:

unsigned short int int_sqrt32(unsigned int x) { unsigned short int res = 0; unsigned int res2 = 0; unsigned short int add = 0x8000; unsigned int add2 = 0x80000000; int i; for(i = 0; i < 16; i++) { unsigned int g2 = res2 + (res << i) + add2; if (x >= g2) { res |= add; res2 = g2; } add >>= 1; add2 >>= 2; } return res; }

También creo que es una mejor idea usar el mismo tipo ( unsigned int ) para todas las variables, ya que de acuerdo con el estándar C, toda la aritmética requiere la promoción (conversión) de tipos de enteros más estrechos al tipo más ancho involucrado antes de que se realice la operación aritmética realizado, seguido de posterior conversión posterior si es necesario. (Esto puede, por supuesto, ser optimizado por un compilador suficientemente inteligente, pero ¿por qué correr el riesgo?)