I would like to know what is an optimal initial guess for use with Newton-Raphson method when finding n-th root.
I develop some program which uses GMP C++ library. GMP manual says:
The initial approximation a[1] is generated bitwise by successively powering a trial root with or without new 1 bits, aiming to be just above the true root.
and i would like to know how it works (what "trial root" is used?) and how to calculate it?
Thanks.
Best Answer
Newton's method for root extraction for $x^n=a$ converges for any initial point. So $a[1]=1$ or $a[1]=1+\frac{a-1}n$ are good initial points.
If you start with too large an $a[1]$, then the initial speed of convergence is linear with factor $1-\frac1n$, so it is important to meet the size of the exact root with the initial root.
For that you can use that $\sqrt[n]{a}=2\sqrt[n]{2^{-n}a}$, also with other factors, to reduce the given number to a neighborhood of $1$.
In a bignumber situation it is fast to determine the bitlength $d$ of the number, since the floating point format stores the exponent $d$ as integer and the mantissa $m$ as array of leaves. If each leaf has $b$ bit, then $2^{d-b}\le a=2^d*(0.m)<2^d$. Use $a[1]=2^{[d/n]}$ or use $\sqrt[n]{a}=2^{[d/n]}\sqrt[n]{2^{-n[d/n]}a}$ to reduce the number under the root.
This trial method mentioned in the GMP manual sounds like the bisection method using $u=2^{[(d-b)/n]}$ and $v=2^{[d/n]+1}$ as the end points of the initial bracketing interval, $u^n\le a\le v^n$. Computing the midpoint $(u+v)/2$ or some point close to it amounts to just setting some lower bits to 1 or removing the leading bit. If after some number of recursions the relative error $\tfrac{(v-u)}v$ becomes smaller than $\tfrac1n$, it is time to switch from bisection to Newton.
Halley's method says to apply Newtons method to $$f(x)=x^{(n+1)/2}-ax^{-(n-1)/2}$$ with derivative $$ f'(x)=\tfrac {n+1}2 x^{(n-1)/2}+a\tfrac{n-1}2 x^{-(n+1)/2} $$ resulting in $$ x_+=x-2x\frac{x^n-a}{(n+1)x^n+(n-1)a}=x\frac{(n-1)x^n+(n+1)a}{(n+1)x^n+(n-1)a} $$ promising cubic convergence with just some operations more.