Ньютон Рафсон с SSE2 - может кто-нибудь объяснить мне эти 3 строки

Я читаю этот документ: http://software.intel.com/en-us/articles/interactive-ray-tracing

и я наткнулся на эти три строки кода:

SIMD-версия уже немного быстрее, но мы можем сделать лучше. Intel добавила быструю 1/sqrt (x) функцию в набор инструкций SSE2. Единственным недостатком является то, что его точность ограничена. Нам нужно точность, поэтому мы уточняем его с помощью 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 ) );

Этот код предполагает существование переменной __m128 с именем 'half' (четыре раза 0,5f) и переменную "три" (четыре раза 3.0f).

Я знаю, как использовать Ньютона Рафсона для вычисления нулевой функции, и я знаю, как использовать его для вычисления квадратного корня числа, но я просто не вижу, как этот код выполняет его.

Может кто-нибудь объяснить это мне, пожалуйста?

Ответы

Ответ 1

Учитывая итерацию Newton y_n+1=y_n(3-x(y_n)^2)/2, это должно быть совершенно прямо, чтобы увидеть это в исходном коде.

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

И если быть точным, этот алгоритм для обратного квадратного корня.

Обратите внимание, что этот по-прежнему не дает полностью точного результата. rsqrtps с итерацией NR дает почти 23 бит точности, против sqrtps 24 бит с правильным округлением для последнего бит.

Ограниченная точность является проблемой, если вы хотите обрезать результат до целого числа. (int)4.99999 - 4. Также обратите внимание на случай x == 0.0 при использовании sqrt(x) ~= x * sqrt(x), потому что 0 * +Inf = NaN.

Ответ 2

Чтобы вычислить обратный квадратный корень из a, метод Ньютона применяется к уравнению 0=f(x)=a-x^(-2) с производной f'(x)=2*x^(-3) и, следовательно, шаг итерации

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

Этот метод без разделения имеет - в отличие от глобально сходящегося метод Heron - ограниченный регион конвергенции, поэтому вам нужно уже хорошее приближение обратного квадратного корня, чтобы получить лучшее приближение.