Newton Raphson与SSE2 – 有人可以解释我这3行
我正在阅读这份文件: http : //software.intel.com/en-us/articles/interactive-ray-tracing
我偶然发现了这三行代码:
SIMD版本已经快了很多,但我们可以做得更好。 英特尔为SSE2指令集添加了快速1 / sqrt(x)函数。 唯一的缺点是它的精度有限。 我们需要精度,所以我们使用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 ) );
此代码假定存在名为“half”(四次0.5f)和变量“three”(四次3.0f)的__m128变量。
我知道如何使用Newton Raphson计算函数的零,我知道如何使用它来计算数字的平方根,但我只是看不出这个代码如何执行它。
有人可以向我解释一下吗?
鉴于牛顿迭代 ,在源代码中看到这一点应该非常简单。
__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 );
确切地说,该算法用于反平方根 。
请注意,这仍然无法提供完全准确的结果 。 具有NR迭代的rsqrtps
提供了近23位的精度,而sqrtps
的24位具有正确的最后一位舍入。
如果要将结果截断为整数,则精度有限是一个问题。 (int)4.99999
是4
。 另外,如果使用sqrt(x) ~= x * sqrt(x)
,请注意x == 0.0
情况,因为0 * +Inf = NaN
。
为了计算a
平方根,牛顿方法应用于方程0=f(x)=ax^(-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方法相比,这种无除法方法具有有限的收敛区域,因此您需要已经很好地逼近逆平方根以获得更好的近似。