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

Question

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?

Answer 1

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) 
def FracPad1(z): return (2*z+1)/(z+2)
def FracPad2(z): return ((14*z+35)*z+5)/((5*z+35)*z+14)
def FracPad3(z): return (((7*z+42)*z+30)*z+2)/(((2*z+30)*z+42)*z+7)

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.

Answer 2

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.

Answer 3

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