La mayoría de las CPUs modernas proporcionan una operación de raíz cuadrada de forma nativa. Así que con un lenguaje de programación típico en un hardware moderno típico, es casi seguro que esa será la operación en última instancia.
Así que vamos a hablar del método que utilizan casi todas las CPUs modernas de propósito general.
Representación de números
Lo primero que tenemos que entender es cómo un ordenador representa realmente el tipo de números de los que estamos hablando. Se almacenan en una especie de notación científica.
Ahora bien, el tipo de notación científica con la que puede estar familiarizado donde un número como -134,5 se representa como [math]-1,345 times 10^{2}[/math].
Cualquier número que no sea cero puede ser representado de esta forma, que conceptualmente tiene tres partes: el signo (es decir, es el número positivo o negativo), la mantisa (un número entre 1 y 10, incluyendo 1 pero sin incluir 10, en este caso 1.345), y el exponente (la potencia a la que se eleva el radix, en este caso 2; nótese que el exponente puede ser negativo o cero en general).
La representación que utilizan casi todos los ordenadores modernos es muy similar, salvo que el número se almacena utilizando una representación binaria en lugar de una decimal. Así, ignorando el bit de signo y el cero (ya que encontrar la raíz cuadrada del cero es trivial, y la raíz cuadrada de un número negativo no es un número real), los números positivos se almacenan en realidad de la forma [math]m times 2^{e}[/math] donde m es un número entre 1 y 2 (incluyendo 1 pero sin incluir 2).
Esta representación se denomina «punto flotante», por la forma en que se puede mover el punto binario (análogo al punto decimal con el que puede estar familiarizado) ajustando el exponente.
La razón por la que todo esto es importante es porque la CPU utiliza la representación para ayudar a calcular raíces cuadradas. Supongamos que el exponente es par. Entonces:
[math]aqrt{m a2epsilon}} = aqrt{m} N-times 2^{epsilon}[/math]
De manera similar, si el exponente es impar, entonces:
[math]qrt{m times 2^{2epsilon+1} = sqrt{2} N – Cuadrado de m. times 2^{epsilon}[/math]
Recuerda que m está en el rango [1,2). Esto significa que hemos reducido el problema de calcular la raíz cuadrada de un número cualquiera a uno de calcular la raíz cuadrada de un número de ese rango. O, si no queremos precalcular [math]sqrt{2}[/math], un número en el rango [1,4). De cualquier manera, hemos simplificado drásticamente el problema, porque ahora podemos usar métodos que se comportan bien en ese rango de números.
Multiplicación y división
Así que ahora que hemos llegado hasta aquí, podrías pensar (si has mirado las otras respuestas) que la solución es ahora usar un método como la iteración Newton-Raphson para calcular la raíz cuadrada. Para calcular la raíz cuadrada de n, elige una conjetura inicial [math]x_0[/math] e itera:
[math]x_{i+1} = frac{1}{2} left( x_i + frac{n}{x_i} right)[/math]
Esta respuesta es incorrecta. Bueno, no es incorrecta en el sentido de que te dará una respuesta correcta. Pero es incorrecta en el sentido de que ningún diseñador de CPU que se precie (o cualquier escritor de bibliotecas si tuviera que implementarlo en software) implementaría la raíz cuadrada de punto flotante de esta manera.
La división por dos es una operación muy barata en binario (que, recuerde, es de base 2, por lo que sólo es desplazar el punto binario un lugar). Sin embargo, la otra división es una operación muy cara, y en este caso, hay que realizar la operación varias veces.
¡La división es tan cara, de hecho, que las CPUs modernas utilizan un algoritmo iterativo que es similar a la iteración Newton-Raphson (pero no realmente) para realizar la división! Está claro que no queremos hacer esto en hardware, en un bucle interno.
En el hardware de los ordenadores modernos, es mucho más barato realizar una operación de multiplicación que una operación de división. La razón es un poco compleja de explicar; tiene que ver con el hecho de que podemos meter muchos más transistores en un chip de lo que solíamos ser capaces, y la multiplicación es un problema que se puede mover en un montón de transistores de manera eficiente. Busca los árboles de Wallace si te interesan los detalles.
En cualquier caso, la cuestión es que la multiplicación es una operación comparativamente barata. Así que si podemos implementar la operación de raíz cuadrada en términos de multiplicación en lugar de división, esto sería mejor.
Newton-Raphson, toma dos
Ahora viene la primera idea clave: En lugar de calcular la raíz cuadrada, calcular el recíproco de la raíz cuadrada. Es decir, en lugar de [math]qrt{n}[/math], calcular [math]frac{1}{qrt{n}}[/math]. Resulta que este es un número mucho más fácil de calcular, y si necesitas la raíz cuadrada, multiplica este número por n y ya está.
El método Newton-Raphson para la raíz cuadrada recíproca es así. Si [math]x_0[/math] es una estimación inicial a [math]frac{1}{sqrt{n}[/math], iterar:
[math]x_{i+1} = frac{1}{2} x_i left( 3 – n {x_i}^2 right)[/math]
De nuevo, la división por dos es bastante barata, y todo lo demás es multiplicación y suma/resta.
Esta es una gran forma de implementarlo en el software. Además, vale la pena señalar que en muchas situaciones prácticas (por ejemplo, normalizar un vector usando el teorema de Pitágoras), la raíz cuadrada recíproca es en realidad más útil que la raíz cuadrada, ya que dividir por una raíz cuadrada es menos eficiente que multiplicar por una raíz cuadrada recíproca.
Sin embargo, no es así como se suele implementar en hardware.
Algoritmo de Goldschmidt
Ahora podemos ver el algoritmo de Goldschmidt, que calcula la raíz cuadrada y la raíz cuadrada recíproca juntas.
Dado [math]b_0 = n[/math], si podemos encontrar una serie de números [math]Y_i[/math] tal que [math]b_n = b_0 {Y_0}^2 {Y_1}^2 cdots {Y_{n-1}}[/math] se aproxima a 1, entonces [math]y_n = {Y_0} {Y_1} cdots {Y_{n-1}[/math] se acercará a [math]frac{1}{sqrt{b_0}[/math] y [math]x_n = b_0 {Y_0} {Y_1} cdots {Y_{n-1}[/math] se acercará a [math]sqrt{b_0}[/math].
La serie que utilizamos es esencialmente el paso de actualización de Newton-Raphson:
[math]begin{align*}b_i & = b_{i-1} {Y_{i-1}^2 frac{1}{2}(3 – b_i)final{align*}[/math]
Y luego podemos seguir la raíz cuadrada:
[math]begin{align*}x_0 & = n Y_0 \b} & = x_{i-1} Y_iend{align*}[/math]
Y la raíz cuadrada recíproca:
[math]begin{align*}y_0 & = Y_0 \ y_{i} & = y_{i-1} Y_i end{align*}[/math]
Pero si llevamos la cuenta de [math]x_i[/math] y [math]y_i[/math], observa que [math]b_i = x_{i-1} y_{i-1}[/math]. Así que en realidad nunca tenemos que llevar la cuenta de [math]b_i[/math]:
[math]begin{align*} Y_i & = frac{1}{2} left(3 – b_iright) & = 1 + frac{1}{2} Izquierda(1 – b_i\right) & = 1 + frac{1}{2} left(1 – x_{i-1} y_{i-1}right) end{align*}[/math]
Y ahora, tampoco necesitamos seguir la pista de [math]Y_i[/math]:
[math]begin{align*} x_i & = x_{i-1} left( 1 + frac{1}{2} left(1 – x_{i-1} y_{i-1}right) right) \N y = x_{i-1} + x_{i-1} frac{1}{2} left(1 – x_{i-1} y_{i-1}right) y_i & = y_{i-1} left( 1 + frac{1}{2} izquierda(1 – x_{i-1} y_{i-1}decha) derecha) & = y_{i-1} + y_{i-1} frac{1}{2} left(1 – x_{i-1} y_{i-1}right)end{align*}[/math]
Aquí hay algún cálculo redundante, que podemos eliminar, sugiriendo el siguiente algoritmo. Dada [math]Y[/math] como una aproximación a [math]frac{1}{sqrt{n}}[/math], establecer:
[math]begin{align*} x_0 & = n Y y_0 & = Yend{align*}[/math]
Entonces iterar:
[math]begin{align*} r_i & = frac{1}{2}left( 1 – x_i y_i right) x_{i+1} & = x_i + x_i r_i y{i+1} & = y_i + y_i r_iend{align*}[/math]
Aunque la división por dos es barata, podemos evitarlo también llevando la cuenta de [math]h_i = frac{1}{2} y_i[/math] en lugar de [math]y_i[/math]. Este es el algoritmo de Goldschmidt.
Supongamos que [math]Y[/math] es una aproximación a [math]frac{1}{sqrt{n}[/math]. Establecer:
[math]begin{align*} x_0 & = n Y h_0 & = frac{Y}{2}{align*}[/math]
Entonces iterar:
[math]begin{align*} r_i & = frac{1}{2} – x_i h_i \Nde x_{i+1} & = x_i + x_i r_i \Nh_{i+1} & = h_i + h_i r_i\Nfinalizar{{math]
Entonces [math]x_i[/math] converge a [math]qrt{n}[/math] y [math]h_n[/math] converge a [math]frac{1}{2sqrt{n}}[/math].
Implementando esto en hardware
Hasta aquí todo bien. Ciertamente es un algoritmo muy sencillo. Pero, ¿por qué es mejor?
Las CPUs modernas suelen tener un circuito rápido que realiza una operación optimizada de multiplicación-acumulación, a menudo llamada multiplicación-agregación fusionada, o FMA para abreviar. Si ha buscado la referencia a los árboles de Wallace antes, debería estar claro cómo FMA podría ser casi tan eficiente como una operación de multiplicación directa.
FMA es una de las primitivas más útiles para tener alrededor. Si necesitas evaluar un polinomio, por ejemplo:
[math]p(x) = a_0 + a_1 x + a_2 x^2 + cdots a_n x^n[/math]
puedes usar la regla de Horner, que es esencialmente un montón de FMA:
[math]begin{align*}p_{n-1} & = a_{n-1} + x a_n \n p_{n-2} & = a_{n-2} + x p_{n-1} & vdots p_1 & = a_1 + x p_2 \ p_0 & = a_0 + x p_1end{align*}[/math]
El bucle interno del algoritmo de Goldschmidt son tres FMAs y nada más. Por eso es una ventaja: sólo necesitas un tipo de circuito (posiblemente sólo un circuito; ten en cuenta que los dos segundos FMAs son independientes y por tanto pueden beneficiarse del pipelining) más algo de lógica de control para implementarlo todo. Y es un circuito que es útil en muchas otras operaciones, por lo que no está desperdiciando una gran cantidad de hardware sólo en la operación de raíz cuadrada.
La penúltima pieza del rompecabezas es cómo obtener una buena estimación inicial, y la respuesta corta es que el mejor método es utilizar una tabla de búsqueda. Incluso una tabla de tamaño modesto, porque sólo estamos buscando raíces cuadradas en un rango tan pequeño, es bastante eficiente.
La última pieza del rompecabezas es: ¿Cómo sabemos cuándo hemos terminado de iterar? Y la respuesta a eso es que conocemos la precisión de los números que estamos usando, y lo buena que es la estimación inicial. A partir de eso, podemos calcular el número máximo de iteraciones necesarias por adelantado. Así que en realidad no hacemos un bucle como tal, sólo hacemos la iteración un número fijo de veces.