快速反向任意幂根算法实现

Question

许多资料表明，众所周知的fast inverse square root算法可以推广到计算任意功率逆根。不幸的是我没有找到这样的C ++实现，而且我不太擅长数学来自己概括这个方法。你能帮我做这个或者提供一个现成的解决方案吗？我认为这对许多人来说很有用，特别是有很好的解释。

这是原始算法，我不太明白我需要更改以获取例如1 /cbrt(x)：

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the...? 
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

    return y;
}

Answer 1

该算法包括两个步骤 - 粗略的解决方案估计和使用几个Newton method步骤的解决方案改进。

Rough estimation

基本思想是使用float number logarithm log2(x)与其整数表示Ix之间的关系：

（图片来自https://en.wikipedia.org/wiki/Fast_inverse_square_root）

现在使用着名的log​​arithm identity作为root：

.

结合前面获得的身份，我们得到：

用数值代替L * (B - s) = 0x3F7A3BEA，等等

Iy = 0x3F7A3BEA / c * (c + 1) - Ix / c;。

对于一个简单的浮点数表示为整数和后退，使用union类型很方便：

   union 
   { 
      float f; // float representation
      uint32_t i; // integer representation
   } t;

   t.f = x;
   t.i = 0x3F7A3BEA / n * (n + 1) - t.i / n; // Work with integer representation
   float y = t.f; // back to float representation

请注意，对于n=2，表达式被简化为t.i = 0x5f3759df - t.i / 2;，这与原始的i = 0x5f3759df - ( i >> 1 );相同

Newton's solution improvement

转变平等

进入一个应该解决的方程式：

现在构建牛顿步骤：

以编程方式看起来像：y = y * (1 + n - x * pow(y,n)) / n;。作为初始y，我们使用在粗略估计步骤中获得的值。

注意对于平方根（n = 2）的特殊情况，我们得到y = y * (3 - x*y*y) / 2;，其与原始公式y = y * (threehalfs - (x2 * y * y))相同;

最终代码作为模板功能。参数N确定根功率。

template<unsigned N>
float power(float x) {
   if (N % 2 == 0) return power<N / 2>(x * x);
   else if (N % 3 == 0) return power<N / 3>(x * x * x);
   return power<N - 1>(x) * x;
};

template<>
float power<0>(float x){ return 1; }

// fast_inv_nth_root<2>(x) - inverse square root 1/sqrt(x)
// fast_inv_nth_root<3>(x) - inverse cube root 1/cbrt(x)

template <unsigned n>
float fast_inv_nth_root(float x)
{
   union { float f; uint32_t i; } t = { x };

   // Approximate solution
   t.i = 0x3F7A3BEA / n * (n + 1) - t.i / n;
   float y = t.f;

   // Newton's steps. Copy for more accuracy.
   y = y * (n + 1 - x * power<n>(y)) / n;
   y = y * (n + 1 - x * power<n>(y)) / n;
   return y;
}

Testing

测试代码：

int main()
{
   std::cout << "|x          ""|fast2      "" actual2    "
      "|fast3      "" actual3    "
      "|fast4      "" actual4    "
      "|fast5      "" actual5    ""|" << std::endl;

   for (float i = 0.00001; i < 10000; i *= 10)
      std::cout << std::setprecision(5) << std::fixed
      << std::scientific << '|'
      << i << '|'
      << fast_inv_nth_root<2>(i) << " " << 1 / sqrt(i) << "|"
      << fast_inv_nth_root<3>(i) << " " << 1 / cbrt(i) << "|"
      << fast_inv_nth_root<4>(i) << " " << pow(i, -0.25) << "|"
      << fast_inv_nth_root<5>(i) << " " << pow(i, -0.2) << "|"
      << std::endl;
}

结果：

|x          |fast2       actual2    |fast3       actual3    |fast4       actual4    |fast5       actual5    |
|1.00000e-05|3.16226e+02 3.16228e+02|4.64152e+01 4.64159e+01|1.77828e+01 1.77828e+01|9.99985e+00 1.00000e+01|
|1.00000e-04|9.99996e+01 1.00000e+02|2.15441e+01 2.15443e+01|9.99991e+00 1.00000e+01|6.30949e+00 6.30957e+00|
|1.00000e-03|3.16227e+01 3.16228e+01|1.00000e+01 1.00000e+01|5.62339e+00 5.62341e+00|3.98103e+00 3.98107e+00|
|1.00000e-02|9.99995e+00 1.00000e+01|4.64159e+00 4.64159e+00|3.16225e+00 3.16228e+00|2.51185e+00 2.51189e+00|
|1.00000e-01|3.16227e+00 3.16228e+00|2.15443e+00 2.15443e+00|1.77828e+00 1.77828e+00|1.58487e+00 1.58489e+00|
|1.00000e+00|9.99996e-01 1.00000e+00|9.99994e-01 1.00000e+00|9.99991e-01 1.00000e+00|9.99987e-01 1.00000e+00|
|1.00000e+01|3.16226e-01 3.16228e-01|4.64159e-01 4.64159e-01|5.62341e-01 5.62341e-01|6.30948e-01 6.30957e-01|
|1.00000e+02|9.99997e-02 1.00000e-01|2.15443e-01 2.15443e-01|3.16223e-01 3.16228e-01|3.98102e-01 3.98107e-01|
|1.00000e+03|3.16226e-02 3.16228e-02|1.00000e-01 1.00000e-01|1.77827e-01 1.77828e-01|2.51185e-01 2.51189e-01|
|1.00000e+04|9.99996e-03 1.00000e-02|4.64155e-02 4.64159e-02|9.99995e-02 1.00000e-01|1.58487e-01 1.58489e-01|