Newton Raphson con SSE2 - qualcuno mi può spiegare queste 3 righe

Question

sto leggendo questo documento: http://software.intel.com/en-us/articles/interactive-ray-tracing

e sono incappato in queste tre righe di codice:

La versione SIMD è già un bel un po 'più veloce, ma possiamo fare di meglio. Intel ha aggiunto una funzione veloce 1/sqrt (x) al set di istruzioni SSE2. L'unico inconveniente è che la sua precisione è limitata. Abbiamo bisogno della precisione , quindi affiniamo utilizzando 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)); 

Questo codice presuppone l'esistenza di una variabile __m128 denominata 'metà' (quattro volte 0.5f) e una variabile ' tre '(quattro volte 3.0f).

So come utilizzare Newton Raphson per calcolare lo zero di una funzione e so come utilizzarlo per calcolare la radice quadrata di un numero ma non riesco a vedere come questo codice lo esegue.

Qualcuno può spiegarmelo per favore?

Answer 1

Data l'iterazione di Newton , dovrebbe essere abbastanza semplice vederlo nel codice sorgente.

__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 
); 

E per essere precisi, questo algoritmo è per the inverse square root.

Si noti che questo still doesn't give fully a fully accurate result. rsqrtps con un'iterazione NR fornisce quasi 23 bit di precisione, rispetto ai 24 bit di sqrtps con arrotondamento corretto per l'ultimo bit.

La precisione limitata è un problema se si desidera truncate the result to integer. (int)4.99999 è 4. Inoltre, fare attenzione al caso x == 0.0 se si utilizza sqrt(x) ~= x * sqrt(x), perché 0 * +Inf = NaN.

Answer 2

Per calcolare l'inverso radice quadrata di a, metodo di Newton viene applicata l'equazione 0=f(x)=a-x^(-2) con derivata f'(x)=2*x^(-3) e quindi l'iterazione passo

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

Questo metodo privo di divisione ha - in contrasto con la convergenza globale Heron's method - una regione limitata di convergenza, quindi è necessaria una buona approssimazione della radice quadrata inversa per ottenere una migliore approssimazione.