# Numerical Stability of Halley's Recurrence for Integer $n^{\mathrm{th}}$-Root

tl;dr? See last paragraph.

If I use the initial value $2^{\left(\big\lfloor\lfloor\log_2 x \rfloor/n\big\rfloor + 1\right)}$ with Halley's recurrence in the compact form

$x_{k+1} = \frac{x_k\Big[A\left(n+1\right) + x_k^n\left(n-1\right) \Big] }{A\left(n-1\right) + x_k^n\left(n+1\right)}$

to evaluate $\lfloor x^{1/n}\rfloor$ with $x,n \in \mathbb{N}$, $x \gg 1$ and $n > 2$ it seems (empirical tests only!) to work. Slowly.

As is the case with all of these methods: the closer the initial value $x_0$ to the actual root, the smaller the amount of iterations needed. Many papers have been written about it, although not many for the integer versions, but a simple refinement can be implemented by observing that the root lies between $2^{\lfloor \lfloor\log_2 x \rfloor / n \rfloor + 1}$ and $2^{\lfloor \lfloor\log_2 x \rfloor / n \rfloor}$ so a simpel arithmetic average of these limits should give a significant decrease in the number of iterations needed and, lo! and behold, it does. Small problem: the algorithm is unstable with that seed. Visible in the following pretty picture with a highly abused y-axis (#iterations, bitsize of root, and absolute error) and the index along the x-axis. The radicand used was $5987^{797}$.

The range of indices where the error occurs is outside the range where I would use Halley's recurrence and change to bisection. The cut-off point I have choosen is the intersection between the bisection which is linear $ax^{-1}$ and the approximately linear part of the Halley iterations at the beginning $bx$ which puts the intersection at $\sqrt{a/b}$. Some runs with up to $3\,321\,960$ bits (ca. one million decimal digits) showed $\Big\lfloor\sqrt{A_b\left(\left\lfloor\log_2\left(\left\lfloor\log_2 \left(A_b\right) \right\rfloor\right) \right\rfloor + 1\right)}\Big\rfloor$ with $A_b = \left\lfloor\log_2 \left(A\right) \right\rfloor$ to be a good estimate for the big-integer library in use.

Hence my question: is Halley's recurrence, implemented as described above, numerically stable in the range $(3,\Big\lfloor\sqrt{A_b\left(\left\lfloor\log_2\left(\left\lfloor\log_2 \left(A_b\right) \right\rfloor\right) \right\rfloor + 1\right)}\Big\rfloor)$ with the initial value the arithmetic average of $2^{\lfloor \lfloor\log_2 x \rfloor / n \rfloor + 1}$ and $2^{\lfloor \lfloor\log_2 x \rfloor / n \rfloor}$ or not and, much more interesting, why?