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}$.

Numebr of Iterations and Error

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?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.