Find rational power of two using Newton’s method

approximationcalculusexponentiationnewton raphson

I know how to find integer powers of two, $2^x = \prod_{i=1}^x 2$, I have memorized powers of two up to 32nd power of two, and I use bit-shifts to calculate them. For integer powers of other numbers I use exponentiation by squaring.

I want to find rational powers of two, without using built-in functions such as pow and exp2 and such, I find them slow and I want to beat them.

My idea is simple, Double precision floating point is encoded as a 64 bits binary number, the first bit is the sign bit, the next 11 bits are used to encode the exponent, 1023 is added to the exponent to make the exponents nonnegative, and the last 52 bits is used to encode the fraction.

My idea is very simple, extract the exponent and mantissa from the binary, use a simple polynomial to approximate it, use some simple mathematics to get the approximated value, and use Newton's method to refine the approximation.

Because the mantissa $m \in [1, 2)$, I use a mean square error polynomial approximation of 6th order to approximate it:

$2 ^ x \approx 4.37425075 \cdot 10^{-4} \cdot x^6 -1.48068542 \cdot 10^{-4} \cdot x^5 + 1.35513429 \cdot 10^{-2} \cdot x^4 + 4.94846338 \cdot 10^{-2} \cdot x^3 + 2.45619546 \cdot 10^{-1} \cdot x^2 + 6.90511103 \cdot 10^{-1} \cdot x + 1.00054403, x \in [1, 2)$.

The above polynomial gives at least 8 correct decimal digits, with 0.8385 of times 9 digits correct.

After calculating that, I use simple mathematics and bit shifts to assemble the final approximation, which isn't important to this question.

And now I want to do a few iterations of Newton's method to improve the accuracy of the result.

Newton's method:

$x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}$

My idea is simple, convert the exponent to a fraction, then I just have to solve this equation:

$x = 2^{\frac{a}{b}}$

Simple power law:

$x^b = 2^a$

A simple function: $f(x) = x^b – 2^a$. $2^a$ is a constant and can be calculated by 1 << a (a is integer here), b is integer so x^b can be calculated by exponentiation by squaring.

The problem is that I don't know what the derivative of the function is. Don't laugh at me, I am a high school dropout and my calculus is very bad.

I have looked at the Wikipedia article for derivative and I guess it is:

$f'(x) = b \cdot x^{b – 1}$

If so, then the iteration scheme is:

$x_{n + 1} = x_n – \frac{x^b – 2^a} {b \cdot x^{b – 1}}$.

Am I right? Will this work?

Best Answer

It will indeed work.

You seek a solution to the equation

$x=2^{p/q}\implies x^q=2^p; p\in\mathbb{Z}, q\in\mathbb{Z}_{>0}$

Newton's method gives, as you expect,

$x'=x-\dfrac{x^q-2^p}{qx^{q-1}}=x-\dfrac{x(x^{q-1})-2^p}{qx^{q-1}}$

where $x^{q-1}$ with a whole number exponent is rendered with $O(\log q)$ multiplications via the squaring and multiplication method.

Related Question