Divisionless iterative approximation method of cube root

It seems that the way cbrt is computed in SLEEF is not very common.

Applying Newton's method to the equation (1/x^3) - a = 0 produces a method that converges to the cube root of 1/a:

x_{n+1} = \frac{1}{3} \left( 4x_n - ax_n^4 \right).

Then, an approximation of the cube root of a can be computed as ax_n^2.

The whole code is like the following.

double xcbrt(double a) {
  int e = ilogb(fabs(a));
  e -= e % 3;
  a = ldexp(a, -e);

  double x = ((-0.0014755156 * a + 0.029122174) * a - 0.21788529) * a + 1.1285767;
  for(int i=0;i<4;i++) {
    x = (4 * x - (x * x) * (x * x) * a) * (1.0 / 3.0);
  }

  return ldexp(a * x * x, e / 3);
}

In this code, the argument is first reduced to [1, 8), and the reciprocal of cube root of the reduced argument is approximated with a polynomial. Then, the method described above is applied to obtain a better approximation.