algorithm - mundo - ¿Cuál es el algoritmo más rápido de división en lenguaje ensamblador x86-64 para grandes números?
nasm ensamblador (5)
Estoy escribiendo una biblioteca de códigos en x86-64 en lenguaje ensamblador para proporcionar todas las funciones convencionales bitwise, shift, logic, compare, arithmetic y math para s0128 , s0256 , s0512 , s1024 tipos de enteros con signo y f0128 , f0256 , f0512 , f1024 flotante tipos de puntos Hasta ahora estoy trabajando en los tipos enteros con signo, porque las funciones de punto flotante probablemente llamarán a algunas rutinas internas escritas para los tipos enteros.
Hasta ahora he escrito y probado funciones para realizar los diversos operadores unarios, comparar operadores y sumar, restar y multiplicar operadores.
Ahora estoy tratando de decidir cómo implementar funciones para los operadores de división.
Mi primer pensamiento fue: "Newton-Raphson debe ser el mejor enfoque". ¿Por qué? Debido a que converge muy rápidamente dado un buen valor de inicio (conjetura inicial), y me imagino que debería ser capaz de descubrir cómo ejecutar la instrucción de división de 64 bits nativa en los operandos para obtener un excelente valor de inicio. De hecho, si el valor de inicialización es preciso a 64 bits, para obtener las respuestas correctas solo debe tomar:
`s0128` : 1~2 iterations : (or 1 iteration plus 1~2 "test subtracts")
`s0256` : 2~3 iterations : (or 2 iterations plus 1~2 "test subtracts")
`s0512` : 3~4 iterations : (or 3 iterations plus 1~2 "test subtracts")
`s1024` : 4~5 iterations : (or 4 iterations plus 1~2 "test subtracts")
Sin embargo, un poco más de reflexión sobre esta pregunta me hace pensar. Por ejemplo, recuerdo la rutina central que escribí que realiza la operación de multiplicación para todos los tipos de enteros grandes:
s0128 : 4 iterations == 4 (128-bit = 64-bit * 64-bit) multiplies + 12 adds
s0256 : 16 iterations == 16 (128-bit = 64-bit * 64-bit) multiplies + 48 adds
s0512 : 64 iterations == 64 (128-bit = 64-bit * 64-bit) multiplies + 192 adds
s1024 : 256 iterations == 256 (128-bit = 64-bit * 64-bit) multiplies + 768 adds
El crecimiento en las operaciones para los tipos de datos más amplios es bastante sustancial, aunque el bucle es bastante corto y eficiente (incluido el almacenamiento en caché). Este bucle escribe cada porción de 64 bits del resultado solo una vez, y nunca lee ninguna porción del resultado para un procesamiento posterior.
Esto me hizo pensar si los algoritmos de división de tipo de cambio y resta más convencionales podrían ser más rápidos, especialmente para los tipos más grandes.
Mi primer pensamiento fue este:
result = dividend / divisor // if I remember my terminology
remainder = dividend - (result * divisor) // or something along these lines
# 1: para calcular cada bit, generalmente el divisor se resta del dividendo SI el divisor es menor o igual al dividendo. Bueno, generalmente podemos determinar que el divisor es definitivamente menor o mayor que el dividendo solo inspeccionando sus porciones más significativas de 64 bits. Solo cuando esas porciones de ms64 bits son iguales, la rutina debe verificar las siguientes porciones inferiores de 64 bits, y solo cuando sean iguales debemos verificar incluso más abajo, y así sucesivamente. Por lo tanto, en casi todas las iteraciones (calculando cada bit del resultado), podemos reducir en gran medida las instrucciones ejecutadas para calcular esta prueba.
# 2: Sin embargo ... en promedio, alrededor del 50% del tiempo encontraremos que debemos restar el divisor del dividendo, por lo que tendremos que restar todo su ancho de todos modos. En este caso, en realidad ejecutamos más instrucciones de las que tendríamos en el enfoque convencional (donde primero las restamos y luego probamos los indicadores para determinar si el divisor <= dividendo). Por lo tanto, la mitad del tiempo nos damos cuenta de un ahorro y la mitad del tiempo nos damos cuenta de una pérdida. En los tipos grandes, como s1024 (que contiene componentes de -16- 64 bits), los ahorros son sustanciales y las pérdidas son pequeñas, por lo que este enfoque debería lograr un gran ahorro neto. En tipos pequeños como s0128 (que contiene componentes de 2 a 64 bits), los ahorros son pequeños y las pérdidas significativas pero no enormes.
Entonces, mi pregunta es, "cuáles son los algoritmos de división más eficientes", dado:
#1: modern x86-64 CPUs like FX-8350
#2: executing in 64-bit mode only (no 32-bit)
#3: implementation entirely in assembly-language
#4: 128-bit to 1024-bit integer operands (nominally signed, but...)
NOTA: Mi conjetura es que la implementación real operará solo en enteros sin signo. En el caso de la multiplicación, resultó ser más fácil y más eficiente (quizás) convertir los operandos negativos en positivos, luego realizar la multiplicación sin signo, luego negar el resultado si exactamente un operando original era negativo. Sin embargo, consideraré un algoritmo entero con signo si es eficiente.
NOTA: Si las mejores respuestas son diferentes para mis tipos de punto flotante ( f0128 , f0256 , f0512 , f1024 ), explique por qué.
NOTA: Mi rutina interna sin signo de multiplicación múltiple, que realiza la operación de multiplicación para todos estos tipos de datos enteros, produce un resultado de doble ancho. En otras palabras:
u0256 = u0128 * u0128 // cannot overflow
u0512 = u0256 * u0256 // cannot overflow
u1024 = u0512 * u0512 // cannot overflow
u2048 = u1024 * u1024 // cannot overflow
Mi biblioteca de códigos ofrece dos versiones de multiplicación para cada tipo de datos de entero con signo:
s0128 = s0128 * s0128 // can overflow (result not fit in s0128)
s0256 = s0256 * s0256 // can overflow (result not fit in s0256)
s0512 = s0512 * s0512 // can overflow (result not fit in s0512)
s1024 = s1024 * s1024 // can overflow (result not fit in s1024)
s0256 = s0128 * s0128 // cannot overflow
s0512 = s0256 * s0256 // cannot overflow
s1024 = s0512 * s0512 // cannot overflow
s2048 = s1024 * s1024 // cannot overflow
Esto es consistente con la política de mi biblioteca de código para "nunca perder precisión" y "nunca desbordar" (los errores se devuelven cuando la respuesta no es válida debido a pérdida de precisión o debido a desbordamiento / desbordamiento). Sin embargo, cuando se invocan funciones de valor de retorno de doble ancho, no se pueden producir tales errores.
¿Seguro que conoce los paquetes de precisión arbitrarios existentes (por ejemplo, http://gmplib.org/ ) y cómo funcionan? Por lo general, están diseñados para ejecutarse "lo más rápido posible" para precisiones arbitrarias.
Si los especializase en tamaños fijos (p. Ej., Aplicó técnicas de evaluación parcial [manualmente] para plegar constantes y desenrollar bucles) esperaría que obtenga rutinas bastante buenas para precisiones específicas de tamaño fijo del tipo que desea.
Además, si no lo has visto, echa un vistazo a los Algoritmos Seminuméricos de D. Knuth, y antiguo pero realmente bueno, que proporciona algoritmos clave para aritmética de precisión múltiple. (La mayoría de los paquetes se basan en estas ideas, pero Knuth tiene grandes explicaciones y muchísimo derecho).
La idea clave es tratar los números de precisión múltiple como si fueran números de radix muy grande (por ejemplo, radix 2 ^ 64) y aplicar aritmética estándar de 3er grado a los "dígitos" (por ejemplo, palabras de 64 bits). La división consiste en el "dígito del cociente estimado (radix grande), la estimación multiplicada por el divisor, la resta del dividendo, el desplazamiento hacia la izquierda un dígito, la repetición" hasta que obtenga suficientes dígitos para satisfacerlo. Para la división, sí, todo está sin firmar (haciendo el manejo de letreros en envoltorios). El truco básico es estimar bien un dígito de cociente (utilizando las instrucciones de precisión simple que le brinda el procesador), y realizar multiplicaciones rápidas de precisión múltiple por un solo dígito. Ver Knuth para más detalles. Consulte los documentos de investigación técnica sobre aritmética de precisión múltiple (puede hacer algunas investigaciones) para obtener mejoras exóticas ("lo más rápido posible").
La alternativa es la fuerza bruta. Podría tomar los 128 bits más altos de x, dividir por los 64 bits más altos de y, y obtener los 64 bits r más altos del cociente, luego restar rxy de x. Y repita según sea necesario, comprobando cuidadosamente qué tan grandes son los errores.
Las divisiones son bajas. Entonces calculas 2 ^ 127 / (los 64 bits más altos de y) una vez. Luego, para estimar los siguientes 64 bits, multiplique los 64 bits más altos de x con este número y coloque todo en el lugar correcto. La multiplicación es mucho más rápida que la división.
A continuación encontrará que todas estas operaciones tienen largas latencias. Por ejemplo, 5 ciclos para obtener un resultado, pero puedes hacer una multiplicación en cada ciclo. Entonces: Estimar 64 bit del resultado. Comience a restar r * y en el extremo superior de x, para que pueda estimar los siguientes 64 bits lo más rápido posible. Luego restas dos o más productos de x simultáneamente, para evitar la penalización por latencia. Implementar esto es difícil . Algunas cosas pueden no valer la pena incluso para 1024 bits (que son solo dieciséis enteros de 64 bits).
Los enfoques de "radix grande" son más eficientes para los tipos de datos enormes que menciona, especialmente si puede ejecutar 128 bits dividido por instrucciones de 64 bits en lenguaje ensamblador.
Si bien la iteración de Newton-Raphson converge rápidamente, cada iteración requiere un número demasiado grande de multiplicar y acumular pasos para cada iteración.
Para la multiplicación, echa un vistazo aquí:
http://www.math.niu.edu/~rusin/known-math/99/karatsuba
Básicamente, permite realizar una multiplicación de 1024 x 1024 utilizando tres (en lugar de cuatro) multiplicaciones de 512 x 512 bits. O nueve 256 x 256 bits, o veintisiete 128 x 128 bits. La complejidad agregada podría no vencer a la fuerza bruta incluso para 1024 x 1024, pero probablemente para productos más grandes. Ese es el más simple de los algoritmos "rápidos", usando n ^ (log 3 / log 2) = n ^ 1.585 multiplicaciones.
Yo aconsejo no usar ensamblador. Implemente 64 x 64 -> 128 bits de multiplicación con el ensamblador en línea, lo mismo que con add-with-carry (creo que gcc y clang podrían tener operaciones incorporadas para esto hoy en día); luego, por ejemplo, puede multiplicar n bits x 256 bits (cualquier número de palabras por 4 palabras) en paralelo, evitando toda la latencia de la multiplicación, sin volverse loco con el ensamblador.
Para un gran número de bits, aprendí que el algoritmo más rápido es el siguiente: en lugar de dividir x / y, se calcula 1 / y se multiplica por x. Para calcular 1 / y:
1 / y is the solution t of (1 / ty) - 1 = 0.
Newton iteration: t'' = t - f (t) / f'' (t)
= t - (1 / ty - 1) / (-1 / t^2 / y)
= t + (t - t^2 y)
= 2t - t^2 y
La iteración de Newton converge cuadráticamente. Ahora el truco: si quieres una precisión de 1024 bits, comienzas con 32 bits, un paso de iteración da 64 bits, el siguiente paso de iteración da 128 bits, luego 256, luego 512, luego 1024. Así que haces muchas iteraciones, pero solo la última uno usa precisión completa Entonces, en general, haces un 512 x 512-> 1024 productos (t ^ 2), un 1024 x 1024 -> 1024 productos (t ^ 2 y = 1 / y), y otro 1024 x 1024 productos (x * ( 1 / y)).
Por supuesto, tiene que averiguar con mucha precisión cuál es el error después de cada iteración; Probablemente tengas que comenzar con, digamos, 40 bits, perder un poco de precisión en cada paso para que tengas suficiente al final.
No tengo idea en qué momento esto se ejecutaría más rápido que una división de fuerza bruta sencilla como lo aprendiste en la escuela. Y y puede tener menos de la cantidad total de bits.