Calcular la parte fraccionaria de la raíz cuadrada sin tomar la raíz cuadrada

Question

Digamos que tengo un número de $x>0$ y tengo que calcular la parte fraccionaria de su raíz cuadrada:

$$f(x) = \sqrt x-\lfloor\sqrt x\rfloor$$

Si he a $\lfloor\sqrt x\rfloor$ disponible, hay una manera que puedo lograr (o aproximado) sin tener que calcular cualquier otro raíces cuadradas? Me gustaría evitar raíces cuadradas por motivos de rendimiento, ya que suelen ser bastante lenta la operación de calcular.

Como un ejemplo de lo que estoy pensando, aquí es una opción:

$$f(x)\approx\frac{x-\lfloor\sqrt x\rfloor^2}{(\lfloor\sqrt x\rfloor+1)^2-\lfloor\sqrt x\rfloor^2}$$

Como $x$ se hace grande, esta aproximación se vuelve mejor y mejor (menos de 1% de error si $\sqrt x\ge40$), pero para los pequeños $x$, es bastante terrible (sin límites de error como $x$$0$).

Puedo hacer mejor? Preferiblemente un error de menos del 1% para cualquier $\sqrt x<10000$.

El santo grial sería encontrar una fórmula que no necesita ningún tipo de división.

A un lado: En caso de que alguien se interese por qué tengo esto: estoy tratando de ver si puedo acelerar Xiaolin de Wu anti-alias círculo funciones de dibujo mediante el cálculo de las variables necesarias de forma incremental (Bresenham de estilo)

Answer 1

Si usted está usando IEEE 754 compatible con los números de punto flotante, que puede ser que la raíz cuadrada de las operaciones son más rápidas que usted podría suponer, como la raíz cuadrada de la operación es directamente compatible con el estándar y compatible con cualquier hardware está garantizado para redondear correctamente (a diferencia de seno y coseno) y es a menudo el uso de una implementación de hardware de la de Newton-Raphson método de Alijah descrito.

Dicho esto, el algoritmo se ha vinculado a sólo utiliza la parte fraccionaria de la raíz cuadrada para calcular los píxeles de la opacidad, y, en consecuencia, el final de opacidad valor varía de 0 a 255. Debido a la pequeña de la gama, números de punto flotante pueden ser una exageración y un punto fijo entero representación podría ser mejor. Si el rango es realmente sólo un byte y el tamaño máximo de sus radios no son muy grandes, puede utilizar una tabla con un punto fijo de entrada para saltar a la cara de la raíz cuadrada de cálculo. Una de 16 bits número de punto fijo daría una 64KB look-up table, que no es demasiado malo.

Usted también puede ser capaz de evitar una división y raíz cuadrada de la operación para el cálculo de la $45^\circ$ valor (ffd en el algoritmo) mediante la Rápida Inversa de la Raíz Cuadrada hack.

Ahora la pregunta de si existe un método para calcular la parte fraccionaria de una raíz cuadrada a sabiendas de que sólo la parte entera, hay un enfoque que de forma iterativa calcula una raíz cuadrada y minimiza divisiones: Fracciones continuas. El enfoque más simple que conozco para lo que están queriendo hacer (convergencia rápida) es el que se detalla aquí y funciona de la siguiente manera:

$a_0 = x - \lfloor\sqrt x\rfloor^2$

$b_0 = 2\lfloor\sqrt x\rfloor$

$a_n = a_0b_{n-1}$

$b_n = b_0b_{n-1} + a_{n-1}$

que le da cuadráticamente mejor y mejores aproximaciones de la parte fraccional de $\sqrt x$, y se divide $\frac{a_n}{b_n}$ cuando se ha hecho lo suficiente de iteraciones para asegurar la exactitud. Si usted está en última instancia, que sólo necesitan byte de tamaño valores de opacidad, sólo debe tomar 3 o 4 iteraciones o así, y nos ahorramos la división hasta el final, que es la única diferencia significativa entre este y el de Newton-Raphson método que no sea el hecho de que te da la parte fraccionaria directamente.

Si realmente quieren seguir de forma incremental el cálculo de variables tan lejos como puede ir, puede utilizar Gosper continuo de la fracción de algoritmos (véase especialmente la sección sobre las raíces cuadradas) y calcular todas las variables que intervienen como fracciones continuas un término en un tiempo. Esto le permite evitar raíz cuadrada cálculos y divisiones de bits turnos, así como abortar tan pronto como usted sabe que el píxel correcto de opacidad (o cualquier otra cosa que usted está calculando) sin perder el tiempo en el cálculo de los dígitos que usted no necesita, pero implica una seria revisión para el algoritmo que se vincula, y si me fui en detalle esta respuesta se convertiría en un libro.

Así que, esencialmente, si tienes memoria de sobra y la longitud máxima de sus radios no es enorme y su tamaño de salida es pequeño, ir con una tabla con un punto fijo de números. Si quieres sencillez de la aplicación, vaya con punto flotante o fija de números de punto. Si desea evitar absolutamente la raíz cuadrada sin una look-up table ir con Newton-Raphson o la continuación de la fracción variante. Si desea absolutamente minimizar el desperdicio de cómputo a expensas de algunos frente a la sobrecarga de ir con Gosper continuo de la fracción de los algoritmos.

Answer 2

Usted podría tratar de usar el de Newton-Raphson método, donde se calcula la raíz cuadrada de $x$ en su totalidad, utilizando la siguiente actualización (para la iteración $n+1$):- $$y(n+1)=y(n)-\frac{y(n)^2-x}{2y(n)}=0.5\left(y(n)+\frac{x}{y(n)}\right)$$ y establecimiento $y(1)=\lfloor\sqrt x\rfloor$ como valor de partida. Por lo tanto, para la iteración $n+1$, la parte fraccionaria $f(x)$$y(n+1)-\lfloor\sqrt x\rfloor$.

Las ventajas de este enfoque es que el algoritmo converge muy rápidamente (es decir, en alrededor de $4$ $5$iteraciones), y no hay necesidad de calcular cualquier otra de las raíces cuadradas. Desafortunadamente, usted necesita para llevar a cabo una división por iteración.

Answer 3

Después de un poco más de lectura en Ted de Matemáticas del Mundo (primera señalado por @hatch22 en su respuesta), me encontré con una muy bonita manera de llevar a cabo este cálculo utilizando el teorema de Pitágoras:

Dado que algunos estiman $y$ de la raíz cuadrada (en el caso que aquí tenemos$\lfloor\sqrt{x}\rfloor$), si dejamos que la hipotenusa de un triángulo $x+y^2$ y uno de sus otros dos lados $x-y^2$, luego el resto de la cara tendrá la longitud de la $2y\sqrt{x}$, es decir, que contendrá la respuesta deseada.

El ángulo de $\alpha$ formado por el lado de los intereses y la hipotenusa puede ser calculado como:

$$\alpha=\sin^{-1}\frac{x-y^2}{x+y^2}\approx\frac{x-y^2}{x+y^2}$$

donde la aproximación es el primer término de la Serie de Maclaurin para $\sin^{-1}$.

El lateral de interés puede ser calculada a partir de:

$$2y\sqrt{x}=(x+y^2)\cos\alpha\approx(x+y^2)\cos\frac{x-y^2}{x+y^2}\approx(x+y^2)\left(1-\frac{1}{2}\left(\frac{x-y^2}{x+y^2}\right)^2\right)$$

Donde la segunda aproximación son los dos primeros términos de la Serie de Maclaurin para $\cos$.

A partir de esto, ahora podemos obtener:

$$\sqrt{x}\approx\frac{x+y^2}{2y}\left(1-\frac{1}{2}\left(\frac{x-y^2}{x+y^2}\right)^2\right)=\frac{x^2+6xy^2+y^4}{4y(x+y^2)}$$

Para obtener la parte fraccionaria de $\sqrt{x}$ en el rango $0..255$, este puede ser optimizado para:

$$y_{\,\text{square}}=y\times y$$ $$s=x+y_{\,\text{square}}$$ $$r=\frac{(s\times s\ll6) + (x\times y_{\,\text{square}}\ll8)}{s\times y}\,\,\&\,\,255$$

donde $\ll$ significa que un bit a bit de desplazamiento a la izquierda (es decir, $\ll6$ $\ll8$ son equivalentes a $\times\,64$ $\times\,256$ respectivamente) y $\&$ significa que un bit a bit (es decir, $\&\,255$ es equivalente a $\%\,256$ donde $\%$ representa el operador de módulo).

La parte sorprendente es que a pesar de los mínimos de la Serie de Maclaurin utilizado, si podemos usar más cerca de $\lfloor\sqrt{x}\rfloor$ $\lceil\sqrt{x}\rceil$ como la estimación de $y$ (tengo ambos disponibles), la respuesta en el rango de $0..255$ es realmente EXACTO!!! para todos los valores de $x\ge1$ que no conducen a un error por desbordamiento durante el cálculo (es decir, $x<134\,223\,232$ para las de 64 bits enteros y $x<2\,071$ 32 bits enteros).

Es posible ampliar el alcance de utilización de la aproximación a $x<2\,147\,441\,941 $ for 64-bit signed integers and $x<41\,324 dólares para la de 32 bits enteros por el cambio de la fórmula:

$$r=\left(\frac{s\ll6}{y} + \frac{x\times y\ll8}{s}\right)\,\&\,\,255$$

Pero debido a la anterior redondeo, esto conduce a una reducción en la precisión de tal forma que el valor es $1$ en muchos casos.

Ahora el problema: Un poco de benchmarking y la lectura indica que en muchos de los procesadores de una operación de división es en realidad no es mucho más rápido que una raíz cuadrada. Por lo menos que se puede encontrar una manera de deshacerse de la división así, este enfoque no es que realmente me va a ayudar mucho. :(

Actualización:

Si una precisión de $\pm 1$ es aceptable, el rango puede ser aumentado de manera significativa con este cálculo: $$k = \frac{(x + y \times y) \ll 5}{y}$$ $$r = \left(\left(k + \frac{x \ll 12}{k}\right)\ll 1\right)\,\&\,\,255$$

Para 32 bits enteros, esto funciona para cualquier $x<524\,288$, es decir, no se rompe tan pronto como $x \ll 12$ desbordamientos. Así, puede ser utilizado para los círculos hasta el radio 723 píxeles. Nota, que $y$ no cambia en cada paso del algoritmo de Bresenham, por lo $1/y$ puede ser pre-calculado y por lo tanto no agrega un segundo completo de la división para el algoritmo.