Cube Root Approximations Using Square Roots with Implementation Examples

for HP 42S and C99
Cube Root Approximations Using Square Roots
with Implementation Examples

April 4, 2003
Reviced April 16, 2025
Takayuki HOSODA

This article presents cube root approximations for both real and complex numbers, tailored for environments where only square root operations are available. The real-valued implementation uses a three-term Puiseux-like approximation optimized for x ∈ [0.125,8], followed by a Padé-type recurrence with cubic convergence. For the complex-valued version, a broader dynamic range is addressed by constructing the initial estimate using nested square roots, effectively approximating exponentiation in the form x^{m / 2ⁿ}, to ensure stable convergence across magnitudes.

Approximating the Cube Root Using Square Roots

The square root of x is given by x^1/2 and The cube root of x is given by x^1/3. When approximating the exponent 1/3 using the form m / 2ⁿ, where m = round(2ⁿ / 3), the approximation error is given by:
m / 2ⁿ × 3 - 1

which results in an error of ±1/2n (with the negative sign when n is even). The specific values are:
  1 /   2 × 3 - 1 = +0.5
  1 /   4 × 3 - 1 = -0.25
  3 /   8 × 3 - 1 = +0.125
  5 /  16 × 3 - 1 = -0.0625
 11 /  32 × 3 - 1 = +0.03125
 21 /  64 × 3 - 1 = -0.015625
 43 / 128 × 3 - 1 = +0.0078125
 85 / 256 × 3 - 1 = -0.00390625
Fig.1 Approximation error when cube root x ^1/3 is approximated by x ^m/2ⁿ
Among these, it was found that x^11/32 is particularly suitable for achieving an initial approximation error of about half a digit, corresponding to a relative error of approximately ±1/√10 ≈ ±0.316, in the range [1e-10, 1e+10].

$x_0 &=& \sqrt{ \sqrt{a \sqrt{ \sqrt{a \sqrt{a} } } } } = a^{11/32} = a^{0.34375}$

The maximum relative error of the initial approximation formula (1) over the range x ∈ [1e-12, 1e+12] is approximately +0.366 at x = 1e+12, which corresponds to a precision of about 0.437 decimal digits. Since the recurrence relation (6) described later exhibits cubic convergence, four iterations improve the precision by a factor of 3^4 = 81. Thus, the expected precision becomes:
0.437 × 81 ≈ 35.3 digits > 16 digits (IEEE 754 double precision)
Therefore, double-precision accuracy can typically be achieved within four iterations.

Computing the Cube Root Using a Recurrence Formula

Using the initial approximation obtained from equation (1), we refine the cube root using a recurrence relation similar to Newton's method.
$x = \sqrt[3]{a}$
Defining the relative error h as:
$h = 1 - \frac{a}{x^3}$
The cube root function x^1/3 can be expanded into an infinite series for x ≠ 0. Thus, we can express (1 + h)^1/3 as the following series expansion:
$(1 + h)^{1/3} \approx 1 + \frac{1}{3}h - \frac{1}{9} h^2 + \frac{5}{81} h^3 - \frac{10}{243} h^4 + \frac{22}{729} h^5 + O(h^6)$
To obtain a more efficient iterative formula, we replace this series with a (1, 1) Padé approximation, which ensures that the first-order and second-order coefficients match those of the series expansion:
$(1 + h)^{1/3} \approx \frac{2 h + 3}{h + 3}$
From equations (3) and (5), we derive the following recurrence relation:
$x_{i+1} &\leftarrow& x_i\left(\displaystyle\frac{x^3_i + 2a}{2x^3_i + a}\right)$

In the vicinity of the solution, Newton-Raphson iteration exhibits quadratic convergence. However, since the (1,1) Padé approximation matches the series expansion up to the second-order term, this recurrence relation achieves cubic convergence. Moreover, compared to the Newton-Raphson method, which can become unstable when the initial estimate is far from the true value, this Padé-based recurrence formula exhibits better numerical stability in such cases. Although the derivation process is different, the asymptotic formula is equivalent to that by Halley's method.

This method can be extended to the nth root of a (where a ≠ 0), leading to the following asymptotic formula with cubic convergence:

Similarly to Equation (6), when using Padé approximants of order (2, 2) or (3, 3), the recurrence formulas (8) and (9) are obtained. These formulas exhibit fifth- and seventh-order convergence, respectively.
$x_{i+1} &\leftarrow& x_{i}\left(\frac{5x_i^6 + 35ax_i^3 + 14a^2}{14x_i^6 + 35 a x_i^3 + 5a^2}\right)$

$x_{i+1} &\leftarrow& x_{i}\left(\frac{2x^9 + 30ax^6 + 42a^2x^3 + 7a^3}{ 7x^9 + 42ax^6 + 30a^2x^3 + 2a^3}\right)$

Initial approximation error and convergence of asymptotic calculations

Fig.2 cbrt Initial approximatin error and its convergence behavior
(Computed under IEEE 754 Double precision environment)

Implementation Example of Complex Cube Root (HP 42S and C99 Source Code)

Formulas used

Initial approximation :
$x_0 &=& \sqrt{ \sqrt{a \sqrt{ \sqrt{a \sqrt{a} } } } } = a^{11/32} = a^{0.34375}$

Recurrence formula of third order convergence :
$x_\mathrm{n+1} &\leftarrow& x_\mathrm{n} \cdot \frac{2a + x_\mathrm{n}^3}{a + 2x_\mathrm{n}^3}$

cbrt - for HP 42S
Rev.1.01 (March 22, 2025) ccbrt - complex cube root function for C99

00 { 66-Byte Prgm }
01▸LBL "cbrt"
02 ENTER
03 X=0?
04 RTN
05 ENTER
06 ENTER
07 SQRT
08 ×
09 SQRT
10 SQRT
11 ×
12 SQRT
13 SQRT
14▸LBL 00
15 ENTER
16 ENTER
17 X↑2
18 ×
19 ENTER
20 +
21 LASTX
22 RCL+ ST T
23 RCL+ ST T
24 X<>Y
25 RCL+ ST T
26 ÷
27 X<>Y
28 ×
29 LASTX
30 RCL÷ ST Y
31 1
32 -
33 ABS
34 1E-11
35 X>Y?
36 GTO 01
37 R↓
38 R↓
39 GTO 00
40▸LBL 01
41 R↓
42 R↓
43 RTN
44 .END.

/* REFERNCE ONLY:
   FOR EXPLANATION OF THE ALGORITHM AND NOT INTENDED FOR ACTUAL USE.*/

#include <math.h>
#include <float.h>
#include <complex.h>

double complex ccbrt(double complex x);
/* The ccbrt(x) function compute the cube root of x. */

/* ccbrt Rev.1.43 (April 10, 2025) (c) 2003 Takayuki HOSODA */
double complex ccbrt(double complex a) {
    double complex x, xp, c;
    double         r, rp, s;
    if (a == 0) return 0;
    s = 1.0;
    //Simple range reduction
    while (cabs(a) < 0x1.0p-195) { a *= 0x1.0p+195; s *= 0x1.0p-65; }
    while (cabs(a) > 0x1.0p+195) { a *= 0x1.0p-195; s *= 0x1.0p+65; }
    while (cabs(a) < 0x1.0p-39)  { a *= 0x1.0p+39;  s *= 0x1.0p-13; }
    while (cabs(a) > 0x1.0p+39)  { a *= 0x1.0p-39;  s *= 0x1.0p+13; }
    //Initial approximation by a^(11/32)
    x = csqrt(csqrt(a*csqrt(csqrt(a*csqrt(a)))));
    r = 0x1.0p+39;
    do { //Asymptotic formula
        xp = x; rp = r;
        c = x * x * x;
        #if defined( USE_SEVENTH )
        x += x * (a - c) * (5 * (a * a  + c * c) + 17 * a * c )
            / (a * a * (a + a + 30 * c) + c * c * (42 * a + 7 * c));
        #elif defined( USE_FIFTH )
        x += 9 * x * (a + c) * (a - c) 
            / (14 * c * c + 35 * a * c + 5 * a * a);
        #else
        x += x * (a - c) / (c + c + a); //Cubic convergence
        #endif
        r = cabs(x - xp); //Correction radious
    } while (r != 0 && rp > r); //Convergence check by Urabe's method
    return s * x;
}

Note: Practical use may require exception handling and range reduction.

Download : cbrt.raw — (raw program file for Free42)

Implementation Example of Real Cube Root (ANSI C Source Code)

In the case of real numbers, if your environment supports the sqrt function along with ldexp and frexp for directly manipulating the exponent part of floating-point numbers, range reduction becomes straightforward. Using a Puiseux-series-like approximation based on powers of 1/2 (i.e., square roots), we present the following 3-term initial approximation formula (10), optimized for the range reduced to [0.125, 8]:
$x_0 &=& 0.20493 + 0.88852 \sqrt{a} - 0.09345 a, \quad \ x \in [0.125, 8]$

The maximum approximation error of the initial formula (9) is approximately 0.0148 at x ≈ 3.065, which corresponds to about 1.83 decimal digits of precision. Therefore, applying the previously discussed recurrence formula (6) of third-order convergence twice gives: 1.83 × 3² = 16.46 digits of precision, which is sufficient to achieve IEEE 754 double-precision accuracy.

Fig.3: Approximation error in cube root estimation using 3- to 5-term expressions involving square roots (x ∈ [0.125, 8])

If formula (10) is used in 80-bit extended precision arithmetic (such as the x87 floating-point format), applying the cubic-converging recurrence formula (6) twice is not sufficient to reach the full precision of approximately 19 decimal digits. In such cases, the following four-term initial approximation formula (11), which is optimized over the interval [0.125, 8], provides improved accuracy. The maximum relative error of this approximation is approximately 0.0049, corresponding to a precision of about 2.31 decimal digits. Thus, applying the same cubic-converging recurrence formula (6) twice yields:
2.31 × 3² ≈ 20.79 digits > 19 digits (extended precision)
and hence allows the desired accuracy to be achieved.

$x_0 &=& c_1 + c_2 a^{1/2} + c_3 a + c_4 a^{3/2}, \quad x \in [0.125, 8] &(11)\\ \text{where:}&& c_1= 0.166466,\ c_2 = 1.025738, c_3 = - 0.22500,\ c_4 = 0.032796$
cbrt - cube root function for C
#include <math.h>
#include <float.h>

double _cbrt(double d); /* The _cbrt(d) function compute the cube root of d.*/

/* _cbrt Rev.1.9 (April 16, 2025) (c) 2003 Takayuki HOSODA */
double _cbrt(double d) {
    double r, x, a, c;
    int    e;
    if (d == 0) return 0;
    /* Range reduction */
    r = frexp(fabs(d), &e);
    x = ldexp(r, e % 3);
    /* Initial approximation for 0.125 <= x <= 8, |err| < 0.0049 */
    a = 0.166466 + (1.025738 + 0.032796 * x) * sqrt(x) - 0.22500 * x;
    /* Apply the asymptotic formula of cubic convergence twice */
    c = a * a * a; a += a * (x - c) / (c + c + x);
    c = a * a * a; a += a * (x - c) / (c + c + x);
    /* A 1-ulp correction to ensure result consistency. */
    r = fabs(x);
    if (a * a * a * (1.0 - DBL_EPSILON) > r) a *= 1.0 - 0.5 * DBL_EPSILON ;
    if (a * a * a * (1.0 + 2.0 * DBL_EPSILON) < r) a *= 1.0 + DBL_EPSILON ;
    /* Range expansion */
    return ldexp(d > 0.0 ? a : -a, e / 3);
}
Note:
In an empirical evaluation conducted using GCC 7.5, the custom _cbrt function was tested on 1,000,000 randomly generated double-precision floating-point values. The relative error computed as _cbrt(d)³ / d - 1, was found to be within ±4.44089209850063e-16 in all cases, which is approximately two times the machine epsilon for IEEE 754 double precision.

The _cbrt(d) function also supports negative arguments. For d < 0, it returns -_cbrt(-d), consistent with the fact that the cube root is an odd-degree root and thus is defined continuously over the entire set of real numbers.

Appendix

Usage of cbrt (HP 42S)

Operation

[XEQ] "cbrt"   — cube root function

Input

REG X      : argument value a

Output

REG y      : a
REG X      : cbrt(a)

Calculation example of cbrt (HP 42S)

 2 "cbrt" → 1.25992104989
-3 "cbrt" → 0.72112478515 + 1.24902476648 i
0 ENTER 1 COMPLEX "cbrt" → 0.866025403784 + 0.5 i
0 ENTER 10 COMPLEX "cbrt" → 1.86579517236 + 1.07721734502 i
-16 ENTER 16 COMPLEX "cbrt" → 2 + 2 i
1e9 ENTER "cbrt" → 1000

Reference

"Convergence of Numerical Iteration in Solution of Equations.", Urabe Minoru, J. Sci. Hiroshima Univ. Ser. A 19 (1956), no. 3, 479--489. doi:10.32917/hmj/1556071264.
Puiseux series — Wikipedia
Padé approximant — Wikipedia
cbrt.c — glibc