带有 SSE2 的 Newton Raphson - 有人可以解释一下这 3 行吗

Posted 2023-03-12

【中文标题】带有 SSE2 的 Newton Raphson - 有人可以解释一下这 3 行吗

【英文标题】：Newton Raphson with SSE2 - can someone explain me these 3 lines

【发布时间】：2013-01-23 00:40:43

【问题描述】：

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

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

谁能给我解释一下？

【参考方案1】：

考虑到牛顿迭代 ，在源代码中应该很容易看到这一点。

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

准确地说，此算法适用于the inverse square root。

注意这个still doesn't give fully a fully accurate result。 rsqrtps 的 NR 迭代提供了几乎 23 位的准确度，而 sqrtps 的 24 位使用正确的最后一位舍入。

如果您想truncate the result to integer，有限的准确性是一个问题。 (int)4.99999 是 4。另外，如果使用sqrt(x) ~= x * sqrt(x)，请注意x == 0.0 的情况，因为0 * +Inf = NaN。

【讨论】：

当截断为整数时，您认为作为最后一步添加一个与结果具有相同指数但只设置最低一位（或两位？）的值是否可行？有效数？这当然是在最不重要的数字总是低于自己的位置的情况下。

取决于应用程序。关键是当使用迭代方法时sqrt(n*n) == n 并不总是成立。这不能随意“修复”——因为sqrt(n*n - epsilon) == n 可能会导致灾难。

【参考方案2】：

为了计算a 的平方根倒数，将牛顿法应用于具有导数f'(x)=2*x^(-3) 的方程0=f(x)=a-x^(-2) 并因此迭代步骤

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

与全局收敛的 Heron's method 相比，这种无除法方法具有有限的收敛区域，因此您需要一个已经很好的平方根倒数近似值才能获得更好的近似值。