# How to derive this curious approximation to the cube root of $a + bi$?

In this Wikipedia article in Portuguese is given the following approximation for the cube root of a complex number $c = a + bi$:

$$\sqrt[3]{c} \approx k\left ( \frac{29z^3 + 261z^2 + 255z + 22}{7z^3 + 165z^2 +324z+71} \right )$$

where $\sqrt[3]{c} = \sqrt[3]{k^3z},\quad$ $\forall k, z\in \mathbb{C}, z \neq 0$

This gives an approximation, say $p_1$, and then you can use this approximation as a new value of $k$ to get another better approximation $p_2$, and so on.

The question is: how this approximation can be derived? I think that it's an application of Newton's method or some truncated series, but I don't know the details.

Also, what precise is the approximation? Does converges to the cube root for any initial values?

• Looks to me like Padé approximant. – Wojowu Nov 11 '16 at 18:17
• It appears that that formula is original research from 2015 and thus not suitable for Wikipedia. -- Furthermore, it is not clear in what respects that formula is better than Newton's or Halley's methods (Newton for $f(k)=c/k-k^2$, $f'(x)=-c/k^2-2k$, $k_+=k-f/f'=k\frac{2c/k^3+1}{c/k^3+2}=k\frac{2z+1}{z+2}$) in respect of gain in accuracy vs. computational effort. – LutzL Nov 12 '16 at 13:18
• @LutzL Thanks you for your detailed analysis. Then this comes from Padè Approximant? Is the factorization of $c$ which makes this approximation an iterative procedure for getting more and more decimals? – Fractional Inquirer Nov 12 '16 at 14:24
• @FractionalInquirer: There is no report available how that fraction was derived, but no, it is not a Pade approximant. There are infinitely many fractions of degree 3 polynomials that will have error order $O((z-1)^5)$. -- Yes, essentially the method presented in the WP article guesses (badly in some examples) an initial $k_0$ and then iterates $k_{n+1}=k_n·F(c/k_n^3)$ where $F(z)$ is the fraction used. $F(z)=\frac{2z+1}{z+2}$ gives the order 3 Halley method, the other Pade approximants will give correspondingly higher order methods. – LutzL Nov 12 '16 at 14:27

## An iterative approximation of the cube root

The intention of the construction of the fraction $F(z)$ is that if $k=k_0$ is chosen such that $z=c/k^3$ has $|z|\approx 1$ and $|\arg(z)|<\frac\pi3$, then $k_+=k·F(z)$ is a good approximation of $\sqrt[3]c$. One obtains successively more accurate root approximations by iterating $$k_{n+1}=k_n·F(c/k_n^3)$$

In view of that iteration it is sensible to express the error in terms of powers in $x=z-1$.

## On the order of approximation

Let's reverse engineer this. Use, for instance, the Magma CAS online calculator with the commands

PS<x> := PowerSeriesRing(Rationals());
P<z> := PolynomialRing(Rationals());
num := Evaluate(29*z^3+261*z^2+255*z+22 , 1+x); num;
// 567 + 864*x + 348*x^2 + 29*x^3
den := Evaluate(7*z^3+165*z^2+324*z+71, 1+x); den;
// 567 + 675*x + 186*x^2 + 7*x^3
f := num/den + O(x^8); f;
// 1 + 1/3*x - 1/9*x^2 + 5/81*x^3 - 10/243*x^4 + 461/15309*x^5 - 7430/321489*x^6 + 367466/20253807*x^7 + O(x^8)
f^3;
// 1 + x - 1/5103*x^5 + 34/35721*x^6 - 12361/6751269*x^7 + O(x^8)

which shows that the expression is correct to order $O((z-1)^5)$. With degree $3$ in both numerator and denominator and thus $2·4-1$ free coefficients, one could find an expression that is $O((z-1)^7)$, i.e., the Taylor series expressions of both sides match in the first $7$ terms.

However, the loss in error order around the origin might have been used to reduce the maximal error in the disk around $1$, or at least around the segment of the unit circle that is relevant according to the initial considerations.

## Some related balanced Padé approximants

The 2/2 order 5 Padé approximant is $$\sqrt[3]z=\frac{14z^2 + 35z + 5}{5z^2 + 35z + 14} + O((z-1)^5)$$ and the 3/3 order 7 Padé approximant is $$\sqrt[3]z=\frac{7z^3 + 42z^2 + 30z + 2}{2z^3 + 30z^2 + 42z + 7} + O((z-1)^7)$$

## Root iterations comparing the 3 fractions

From WP take the moderately interesting test case $c=11+197i$ with initial guess $k_0=6$. In python define

def FracLud(z): return (((29*z+261)*z+255)*z+22)/(((7*z+165)*z+324)*z+71)

and iterate

c=11+197j

k=6
for _ in range(4): z=c/k**3; k *= FracLud(z); print k, abs(c-k**3)
(5.07502720254+2.80106967751j) 2.5577949438
(5.09495916653+2.81659441015j) 1.39979625158e-11
(5.09495916653+2.81659441015j) 1.71121351099e-13
(5.09495916653+2.81659441015j) 8.54314994314e-14

k=6
for _ in range(4): z=c/k**3; k *= FracPad1(z); print k, abs(c-k**3)
(4.67251486867+3.25849790265j) 60.612721727
(5.09919052684+2.81457943777j) 0.476738798014
(5.09495916512+2.8165944087j) 2.05734637949e-07
(5.09495916653+2.81659441015j) 8.54314994314e-14

k=6
for _ in range(4): z=c/k**3; k *= FracPad2(z); print k, abs(c-k**3)
(5.16659218119+2.75059104398j) 9.95653080095
(5.09495916898+2.81659441176j) 2.97864153108e-07
(5.09495916653+2.81659441015j) 2.89169943033e-14
(5.09495916653+2.81659441015j) 8.54314994314e-14

k=6
for _ in range(4): z=c/k**3; k *= FracPad3(z); print k, abs(c-k**3)
(5.08204659353+2.82547419372j) 1.59145673695
(5.09495916653+2.81659441015j) 8.54314994314e-14
(5.09495916653+2.81659441015j) 2.89169943033e-14
(5.09495916653+2.81659441015j) 8.54314994314e-14

so the order 5 Ludenir/Luderian method converges slightly faster than the computationally faster order 5 Padé2 method, but noticeably slower than the order 7 Padé3 method that has the same computational complexity.

The computationally simplest Halley/Padé1 method needs 4 steps where Padé3 needs 2. In terms of operations, there are 4*(2 div + 2 mul) (discounting addition, multiplication with small constants) versus 2*(2 div + (2+6) mul). Counting 1 complex division as about equal to 2 complex multiplications, this is an equal effort.

• +1 for the neat analysis. FWIW the author credited on the linked page is Ludenir Santos who once posted it on MSE in an answer but appears to have never followed up on it. – dxiv Nov 12 '16 at 20:14
• This is the same original research as in WP, neither really a source for the other. No further derivation given, very weak presentation. The only additional property I found is $F(8)=2$ (with order 1 there), which with order 5 in $z=1$ still leaves one degree of freedom in a fraction of cubic polynomials. – LutzL Nov 12 '16 at 21:49
• @Lutzl. I remember that you have answered me in another question from the same autor here: math.stackexchange.com/questions/1825062/… – Fractional Inquirer Nov 12 '16 at 22:00
• Yes, there he was re-inventing 3rd order methods like Halley's method and related ones. – LutzL Nov 12 '16 at 23:10

Cube Root

We can start with the Taylor Series $$\sqrt[3]{1+x}=1+\frac x3-\frac{x^2}9+\frac{5x^3}{81}-\frac{10x^4}{243}+\frac{22x^5}{729}-\frac{154x^6}{6561}+O\!\left(x^7\right)\tag{1}$$ then use the Extended Euclidean Algorithm to get the Padé Approximation $$\sqrt[3]{1+x}=\frac{81+135x+63x^2+7x^3}{81+108x+36x^2+2x^3}+O\!\left(x^7\right)\tag{2}$$ Now we can simply substitute $x\mapsto x-1$ and get $$f(x)=\frac{2+30x+42x^2+7x^3}{7+42x+30x^2+2x^3}=\sqrt[3]{x}+O\!\left((x-1)^7\right)\tag{3}$$ This approximation is nice in that it gives $\sqrt[\large3]{\frac1x}=\dfrac1{\sqrt[3]x}$. However, the maximum it can return is $\frac72$.

The sequence defined by $u_0=1$ and $u_{n+1}=u_nf\!\left(x/u_n^3\right)$, converges to $\sqrt[3]{x}$, and the fact that the error is $O\!\left((x-1)^7\right)$ means that if $u_n$ is accurate to $d$ places, $u_{n+1}$ should be accurate to $7d$ places. For example, for $x=8$, \begin{align} u_1&=1.981746273197444478247642227&\text{2 places}\\ u_2&=1.999999999999997579278630790&\text{15 places}\\ u_3&=2.000000000000000000000000000 \end{align}

The Approximation in the Question $$\frac{29z^3+261z^2+255z+22}{7z^3+165z^2+324z+71}=\sqrt[3]{z}+O\left((z-1)^5\right)$$ This will not converge as fast as the approximation given in $(3)$ (and uses larger coefficients). I am not sure how this approximation was derived.

Square Root

The same thing can be done to accelerate the square root iteration using $$f(x)=\frac{1+21x+35x^2+7x^3}{7+35x+21x^2+x^3}=\sqrt{x}+O\!\left((x-1)^7\right)$$ The sequence $u_0=1$ and $u_{n+1}=u_nf\!\left(x/u_n^2\right)$ converges in the same manner to $\sqrt{x}$.

• @robjhon Thanks for this answer. Indeed, the fraction that you offer is more simple and fast than the Luderian. But I see that this is something well know, and is not propertly a discover, as the autor of the formula says. Do you agree? – Fractional Inquirer Nov 12 '16 at 21:26
• Yes. Padé approximations are well known. The author of the article may well have discovered the approximation they cite, but the Padé approximation is better. – robjohn Nov 12 '16 at 21:34

COMMENT.- Just as a complement to the exhaustive answer given by @LutzL, we give here the following which I think it is pertinent (is it not, am I wrong?) :

Let $z_0\in\Bbb C$; how many functions $f:\Bbb C\to \Bbb C$ are there such that $f(z_0)=z_0$ and $f$ being a contraction mapping?

Likely there is a lot; if $f$ is one of them then for Banach fixed-point theorem, takink and arbitrary $z\in\Bbb C$, the limit of $$\underbrace{f\circ f\circ\cdots f}_{n\text{ times }}(z)$$ when $n\to\infty$ is $z_0$.

Furthermore if $\underbrace{f\circ f\circ\cdots f}_{n\text{ times }}(z)=z_n$ then the sequence $\{z_n\}$ converges fastly to $z_0$ because of the contractivity of $f$ one has $$d(z_n,z_0)\le \frac{K^n}{1-K}d(z,z_0)$$ where $K\in[0,1)$ is the constant of contraction of $f$ (so the smaller the $K$, the faster the convergence).

What I wanted to mean is there are many approximations to $\sqrt[3]{c}$ analogues to proposed by the O.P. and that the same holds for any complex constant.