¿Cómo calcular la raíz cuadrada de un número?
Por ejemplo, √20 es una notación abreviada y no es necesario escribir el número específicamente.
La raíz cuadrada es el valor de la operación cuadrática. Su operación inversa se multiplica por la segunda potencia.
//
// Calcular el recíproco de la raíz cuadrada del parámetro x
//
float InvSqrt (float x)
{
float xhalf = 0.5f*x;
int i = *(int*)amp;x
i = 0x5f3759df - ( i gt; gt ; 1); // Calcula la primera raíz del valor aproximado
x = *(float*)amp; - xhalf*x*x) ; // Método de iteración de Newton
return x;
}
La esencia del algoritmo es en realidad el método de Newton-Raphson. (NR) basado en la serie de Taylor. NR es un método para encontrar raíces aproximadas de una ecuación. Primero estima un valor más cercano a la raíz de la ecuación, luego infiere el siguiente valor más cercano basándose en la fórmula, y así sucesivamente hasta lograr una precisión satisfactoria. La fórmula es la siguiente:
Función: y=f(x)
La primera derivada es: y'=f'(x)
Entonces n 1 de la ecuación Una raíz aproximada: f(x)=0 es
x[n 1] = x[n] - f(x[n]) / f'(x[n])
< La parte más crítica de p > NR es estimar las primeras raíces aproximadas. Si la raíz aproximada está lo suficientemente cerca de la raíz verdadera, entonces sólo se necesitan unas pocas iteraciones para obtener una solución satisfactoria.Ahora volvamos a ver cómo se puede utilizar el método de Newton para resolver nuestro problema. Encontrar la inversa de una raíz cuadrada es en realidad resolver la ecuación 1/(x^2)-a=0. Expanda esta ecuación usando el método de iteración de Newton de la siguiente manera:
x[n 1]=1/2*x[n]*(3-a*x[n]*x[n])
Poner 1/2 entre paréntesis nos da el inverso de la segunda línea de la función anterior.
A continuación, debemos intentar estimar una primera raíz aproximada. Esta es la parte más mágica de la función anterior. De alguna manera encuentra una raíz aproximada que está muy cerca de la raíz verdadera, por lo que usando solo un proceso de iteración, se obtiene una solución más satisfactoria. ¿Cómo hace esto? Justo en esta línea:
i = 0x5f3759df - (i gt; gt; 1); // Calcula la primera raíz de la aproximación
Declaración súper desconcertante, ¿no? Pero puedes entenderlo si lo piensas detenidamente. Sabemos que según el estándar IEEE, los datos de coma flotante se representan en un sistema de 32 bits como este (a grandes rasgos, pero se omiten muchos detalles, si estás interesado, puedes buscarlo en Google):
bits: 31 30...0
31: bit de signo
30-23: ***8 bits, guardar exponente (E)
22-0: ****23 bits, conservando la mantisa (M)
Por lo tanto, un número de coma flotante de 32 bits representado como un número real decimal es M*2^E. Tome sus raíces y luego corresponda, que es M^(-1/2)*2^(-E/2). Está claro ahora. La declaración igt;gt;1 divide el exponente por 2 para lograr la parte 2^(E/2). El propósito de restarlo con una constante anterior es obtener M^(1/2) invirtiendo los signos de todos los exponentes.
En cuanto a ese 0x5f3759df, sólo puedo decir que efectivamente es un número súper mágico.
Este número mágico se puede derivar, pero no voy a discutirlo aquí porque es demasiado tedioso.
En pocas palabras, funciona así: dado que en los números de punto flotante IEEE, la mantisa M omite el 1 inicial, la mantisa real es 1 M. Si no te quedaste dormido en la clase de matemáticas de la universidad, cuando ves un número de la forma (1 M)^(-1/2), inmediatamente piensas en la expansión de la serie de Taylor, donde el primer término es una constante. La siguiente es una derivación simple:
Para un número real Rgt; 0, suponiendo que en la representación de punto flotante IEEE,
el exponente es E y la mantisa es M, entonces:
R^(-1/2)
= (1 M)^(-1/2) * 2^(-E/2)
Utilice la etapa de Taylor Expanda el número (1 M)^(-1/2) y tome el primer término de la expansión. Expande la serie de Taylor y toma el primer término para obtener:
Ecuación original
= (1-M/2) * 2^(-E/2)
= 2^(-E/2) - (M/2) * 2^(-E/2)
Si no se considera el signo del exponente,
(M/ 2) *2^(E/2) es exactamente (Rgt; gt.1),
Y en la notación IEEE, el signo del exponente simplemente se agrega como un desplazamiento,
Y la primera mitad de la ecuación resulta ser una constante, por lo que la ecuación original se puede convertir a:
Ecuación original = C - (M/2)*2^(E/2 ) = C - (Rgt ;gt; 1), donde C es una constante
Así que solo necesitamos resolver la ecuación:
R^(-1/2)
= ( 1 M)^(-1/2) * 2^(-E/2)
= C - (Rgt; gt; 1)
Encuentre el valor C del error relativo mínimo, ya está.
La derivación anterior es solo mi comprensión personal y no ha sido confirmada. Chris Lomont, por otro lado, analiza en detalle cómo resolver la última ecuación y tratar de encontrar la constante óptima C en una máquina real en su artículo; los lectores interesados pueden encontrar un enlace a su artículo al final de este artículo.
Por lo tanto, el llamado "número mágico" no cayó a la Tierra desde una galaxia en el universo N-dimensional debido a la distorsión del espacio y el tiempo, sino que es una teoría matemática que existe desde hace cientos de años. Siempre que estemos familiarizados con las series NR y Taylor, usted y yo también podemos realizar optimizaciones similares.
Hemos probado en GameDev.net y el error relativo de esta función es de aproximadamente 0,177585, que es más de 20 veces más rápido que la biblioteca estándar de C sqrt. Si se agrega una iteración más, el error relativo se puede reducir al orden de e-004, pero la velocidad también se reducirá a aproximadamente la misma que sqrt. Se dice que en "DOOM3", Carmack optimizó aún más el algoritmo utilizando una tabla de búsqueda, con una precisión casi perfecta y una velocidad mejorada en comparación con el algoritmo original (el código fuente está en estudio, si alguien lo tiene, envíeme una copia). ).
Vale la pena señalar que en el algoritmo de Chris Lomont, la mejor constante teórica (la mayor precisión) es 0x5f37642f y, en las pruebas reales, funciona mejor si solo se usa una iteración. Pero, curiosamente, después de dos iteraciones, la precisión de la solución a esta constante disminuye significativamente (¡Dios sabe qué!). Después de pruebas prácticas, Chris Lomon concluyó que la constante óptima era 0x5f375a86. Si cambia a la versión doble de 64 bits, el algoritmo sigue siendo el mismo, siendo la constante óptima 0x5fe6ec85e7de30da (otro número mágico vergonzoso: b).
Este algoritmo se basa en la representación interna y el orden de bytes de números de coma flotante y, por lo tanto, no es portátil. Si se ejecuta en una Mac, se bloqueará. Si desea portabilidad, debe usar sqrt directamente. Pero las ideas algorítmicas son universales. Puede intentar derivar el algoritmo de raíz cuadrada correspondiente.
El algoritmo de raíz cuadrada utilizado por Carmack en QUAKE3 se detalla a continuación. Carmack ha donado todo el código fuente de QUAKE3 a organizaciones de código abierto, para que puedas usarlo sin preocuparte de recibir una carta de un abogado.
//
// Función de Carmack utilizada para calcular raíces cuadradas en QUAKE3
//
float CarmSqrt(float x){
unión{
intPart;
float floatPart
} convertidor
unión{
int intPart;
float floatPart;
} convertor2;
convertor.floatPart = x; > convertor2.floatPart = x;
convertor.intPart = 0x1FBCF800 (convertor.intPart gt; gt; 1);
convertor2.intPart = 0x5f3759df - (convertor2.intPart gt; gt ; 1);
return 0.5f*(convertor.floatPart (x * convertor2.floatPart));
}
Otra implementación basada en el mismo algoritmo. Un método sqrt más rápido es el siguiente. Simplemente divide el exponente por 2 y no tiene en cuenta la raíz cuadrada de la mantisa. Para entender este código, debes darte cuenta de que en el formato de punto flotante IEEE, E se encuentra sumando 127 al exponente real. Por ejemplo, si el número real es 0,1234*2^10, entonces en notación de punto flotante, el valor de E (bits 23-30) es en realidad 10 127=137. Por lo tanto, el siguiente código trata del desplazamiento 127, que es lo que hace la constante 0x3f800000. En realidad no he probado esta función, por lo que no puedo comentar sobre sus ventajas y desventajas, pero supongo que debería ser mucho menos precisa.
float Faster_Sqrtf(float f)
{
resultado flotante;
_asm
{
mov eax, f
sub eax, 0x3f800000
sar eax, 1
agregar eax, 0x3f800000
resultado de mov , eax
}
devuelve resultado;
}
Además de los métodos basados en NR, otros algoritmos rápidos comunes son las aproximaciones polinómicas. . La siguiente función está tomada de Master Tips for 3D Game Programming y usa polinomios para aproximar la ecuación de longitud original, pero no sé cómo se derivó la fórmula que usó el autor (si lo sabe, hágamelo saber, gracias).
//
// Esta función calcula la distancia de (0, 0) a (x, y), el error relativo es 3,5
//
int FastDistance2D(int x, int y)
{
x = abs(x
y = abs( z); ) relativo El error es 8
//
float FastDistance3D(float fx, float fy, float fz)
{
int temp;
int x, y, z
// Asegúrate de que todos los valores sean positivos<
x = int(fabs(fx) * 1024);
y = int(fabs(fy) * 1024
z = int(fabs(fz) * 1024); p>
// Ordenar
if (y lt; x) SWAP(x, y, temp)
if (z lt; y) SWAP(y, z, temp)
if (y lt; p>
return((float) (dist gt; gt; 10));
}
Hay otro método llamado estimación de distancia? Como se muestra en la siguiente figura:
Los puntos en el octágono regular representado por la línea roja son:
octágono(x, y) = min((1/√2) * (|x | |y|), max(|x|, |y|))
Para encontrar las longitudes de los vectores v1 y v2, necesitamos
√(x^ 2 y^2 ) = (|v1| || v2|)/2 * octagon(x, y)
Hasta ahora, hemos discutido el algoritmo de raíz cuadrada de números de punto flotante, y a continuación es el algoritmo de raíz cuadrada de números enteros. Se podría pensar que no tiene sentido realizar la operación de raíz cuadrada en datos enteros porque daría un resultado como 99^(1/2)=9. Esto es generalmente cierto, pero cuando usamos números de punto fijo (y muchos sistemas todavía usan números de punto fijo, como el GBA de Nintendo y otras computadoras portátiles), la raíz cuadrada del número entero se vuelve muy importante. El algoritmo para elevar al cuadrado números enteros es el siguiente. No lo discutiré aquí (de hecho, no lo he estudiado en detalle ya que no se usará pronto - -b), pero puedes encontrar una derivación muy detallada en el artículo de James Ulery al final del artículo.
//
// Para facilitar la lectura, agregué saltos de línea a la definición de macro a continuación
//
# define paso (shift)
if((0x40000000l gt; gt; shift) sqrtVal lt; = val)
{
val -= (0x40000000l gt; gt; cambio) sqrtVal;
sqrtVal = (sqrtVal gt; 1) | (0x40000000l gt; gt; cambio
}
else
<); p>{sqrtVal = sqrtVal gt;
}
//
// Calcular 32 bits Raíz cuadrada de entero
///
int32 xxgluSqrtFx(int32 val)
{
// Nota: esta función rápida de raíz cuadrada
// Sólo aplicable a Q_FACTOR par
int32 sqrtVal = 0;
paso(0);
paso(4);
paso(6);
paso(8);
paso(10);
p>
p>
paso(12);
paso(14);
paso(16); ;
paso(20);
paso(22);
paso(24); /p>
paso(28);
paso(30);
if(sqrtVal lt; val)
{
sqrtVal;
}
sqrtVal lt; = (Q_FACTOR)/2
retorno(sqrtVal)
}