28

我正在阅读这份文件： 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）的 __m128 变量和一个变量“three”（四倍 3.0f）。

我知道如何使用 Newton Raphson 来计算函数的零，并且我知道如何使用它来计算数字的平方根，但我只是看不出这段代码是如何执行它的。

有人可以向我解释一下吗？

4

2 回答 2

35

鉴于牛顿迭代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 位对最后一位进行正确舍入。

如果要将结果截断为 integer ，则有限的准确性是一个问题。 (int)4.99999是4。另外，x == 0.0如果使用 , 请注意这种情况sqrt(x) ~= x * sqrt(x)，因为0 * +Inf = NaN.

于 2013-02-07T13:59:12.493 回答
3

为了计算 的反平方根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 方法相比，这种无除法方法具有有限的收敛区域，因此您需要一个已经很好的逆平方根近似值才能获得更好的近似值。

于 2014-03-12T16:37:48.767 回答