5
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ An optimal initial guess would be a root. But since you're using a numerical method, I assume you don't know the roots. That means there is no non-trivial optimal initial guesses. $\endgroup$ – Fly by Night May 8 '14 at 21:05
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ (Incidentally, note that $a_1=1+\frac{a-1}{n}$ comes from the binomial approximation $(1+x)^n = 1+nx$ valid for small $x$.) $\endgroup$ – Steven Stadnicki May 8 '14 at 21:20
  • $\begingroup$ Or from the AGM inequality $$\sqrt[n]a=\sqrt[n]{1^{n-1}\,a}\le\frac{(n-1)\,1+a}{n}$$ which is the stricter, the closer $a$ is to $1$. One can extend this via $$\sqrt[n]a=\sqrt[n]{(2^k)^{n-1}\,(2^{-k(n-1)}a)}\le\frac{(n-1)\,(2^k)+2^{-k(n-1)}a}{n}$$ using $k=[d/n]$ as before. $\endgroup$ – LutzL May 8 '14 at 21:31
  • $\begingroup$ Following the details in this answer, stackoverflow.com/a/7180310/3088138, the bits before the decimal point of an mpf_t A; variable can be read as int d = A->_mp_exp; Then using int k=d/n and an mpf_t X which is initialized to the value 1, you can get the trial power of 2 by setting X->_mp_exp = k. $\endgroup$ – LutzL May 8 '14 at 21:45
  • $\begingroup$ @LutzL Thanks for answer. Can you also provide some reference which advices to use bitlength for bignumber situation? $\endgroup$ – Derp May 9 '14 at 16:56
  • $\begingroup$ There should be something in the standard reference, Knuth's Art of Computer Programming, vol. 2. Also see the introduction of H. Cohen "A Course in Computational Algebraic Number Theory", where the basics of the Pari/GP system are explained. No direct source for root extraction algorithms. $\endgroup$ – LutzL May 9 '14 at 17:25

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.