c++ - raiz - Newton Raphson con SSE2-¿alguien me puede explicar estas 3 líneas

raices de un polinomio en c++

Question

Estoy leyendo este documento: http://software.intel.com/en-us/articles/interactive-ray-tracing

y tropecé con estas tres líneas de código:

La versión SIMD es bastante más rápida, pero podemos hacerlo mejor. Intel ha agregado una función rápida 1 / sqrt (x) al conjunto de instrucciones SSE2. El único inconveniente es que su precisión es limitada. Necesitamos la precisión, por lo que la refinamos usando Newton-Rhapson:

 __m128 nr = _mm_rsqrt_ps( x );
__m128 muls = _mm_mul_ps( _mm_mul_ps( x, nr ), nr );
result = _mm_mul_ps( _mm_mul_ps( half, nr ), _mm_sub_ps( three, muls ) ); 

Este código asume la existencia de una variable __m128 llamada ''mitad'' (cuatro veces 0.5f) y una variable ''tres'' (cuatro veces 3.0f).

Sé cómo utilizar Newton Raphson para calcular el cero de una función y sé cómo usarla para calcular la raíz cuadrada de un número, pero no puedo ver cómo lo realiza.

¿Alguien me puede explicar por favor?

Answer 1

Dada la iteración de Newton , debería ser bastante sencillo ver esto en el código fuente.

 __m128 nr = _mm_rsqrt_ps( x ); // The initial approximation y_0
__m128 muls = _mm_mul_ps( _mm_mul_ps( x, nr ), nr ); // muls = x*nr*nr == x(y_n)^2
result = _mm_mul_ps(
_mm_sub_ps( three, muls ) // this is 3.0 - mul;
/*multiplied by */ __mm_mul_ps(half,nr) // y_0 / 2 or y_0 * 0.5
);

Y para ser precisos, este algoritmo es para la raíz cuadrada inversa .

Tenga en cuenta que esto todavía no da un resultado totalmente exacto . rsqrtps con una iteración NR da casi 23 bits de precisión, contra los 24 bits de sqrtps con redondeo correcto para el último bit.

La precisión limitada es un problema si desea truncar el resultado a entero . (int)4.99999 es 4 . Además, tenga cuidado con el caso x == 0.0 si usa sqrt(x) ~= x * sqrt(x) , porque 0 * +Inf = NaN .

Answer 2

Para calcular la raíz cuadrada inversa de a , el método de Newton se aplica a la ecuación 0=f(x)=ax^(-2) con la derivada f''(x)=2*x^(-3) y por lo tanto el paso de iteración

N(x) = x - f(x)/f''(x) = x - (a*x^3-x)/2
= x/2 * (3 - a*x^2)

Este método libre de división tiene, en contraste con el método de Heron convergente a nivel mundial, una región de convergencia limitada, por lo que necesita una buena aproximación de la raíz cuadrada inversa para obtener una mejor aproximación.