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 30, 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- or (preferably) four-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 4-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.2.0 (April 23, 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.1664647 + (1.025745 + 0.032798 * x) * sqrt(x) - 0.2250075 * 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); }

_cbrtl - cube root function for C

#include <math.h> #include <float.h> long double _cbrtl(long double d); /* The _cbrtl(d) function compute the cube root of d.*/ /* _cbrtl Rev.2.0 (April 23, 2025) (c) 2003 Takayuki HOSODA */ long double _cbrtl(long double d) { long double r, x, a, c; int e; if (d == 0) return 0; /* Range reduction */ r = frexpl(fabsl(d), &e); x = ldexpl(r, e % 3); /* Initial approximation for 0.125 <= x <= 8, |err| < 0.0049 */ a = 0.1664647 + (1.025745 + 0.032798 * x) * sqrtl(x) - 0.2250075 * 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 = fabsl(x); if (a * a * a * (1.0 - LDBL_EPSILON) > r) a *= 1.0 - 0.5 * LDBL_EPSILON ; if (a * a * a * (1.0 + 2.0 * LDBL_EPSILON) < r) a *= 1.0 + LDBL_EPSILON ; /* Range expansion */ return ldexpl(d > 0.0 ? a : -a, e / 3); }

Note:
In an empirical evaluation conducted using GCC 7.5, the custom _cbrt and _cbrtl functions were tested on 1,000,000 randomly generated double-precision and long double-precision values, respectively. The relative error, computed as _cbrt(d)³ / d - 1 and _cbrtl(d)³ / d - 1, was found to be within ±4.44089209850063e-16 and ±2.16840434497101e-19, approximately twice the machine epsilon for IEEE 754 double and long double precision, respectively.

Both _cbrt(d) and _cbrtl(d) also support negative arguments: for d < 0, they return -_cbrt(-d) and -_cbrtl(-d), respectively, consistent with the fact that the cube root, being an odd-degree root, is continuously defined over the entire set of real numbers.

Initial Approximation for Real Cube Roots

In the implementation example for real-valued cube roots, the initial estimate uses the approximation formula (10), which incorporates a square root term. This approach is based on the following idea:

First, a Padé approximation of order (2, 0) is applied to the cube root of x² around x ≈ 1. By subsequently replacing x with √x in the resulting expression, an approximation for the cube root of x itself is obtained.
This substitution effectively expands the valid approximation range, as the target domain is squared, leading to reduced approximation errors across a broader region.

The resulting approximation was then adjusted so that it satisfies the condition cbrt(1) = 1, while minimizing the maximum absolute approximation error over the interval x ∈ [0.125, 8]. This yields the final form of approximation formula (10) (see Fig.4).

Fig.4: Approximation of cube root over the squared domain (x ∈ [0.125, 8])

Other Examples of Initial Approximations for Cube Roots

As a reference, several higher-order approximation formulas designed with the same idea are presented below, along with a comparison graph.

These formulas approximate the cube root using 4- to 5-term series expansions that include square root terms.
Approximation formulas:

Fig.5: Approximation errors of initial cube root approximations using 4- to 5-term expansions with square root terms (x ∈ [0.125, 8])

Conclusion

This article has explored a method for approximating the cube root using only square roots, along with implementation examples. For the initial approximation, a Puiseux-like expansion was used in the real-valued version, and an exponential form x^m/2ⁿ was used in the complex-valued version, followed by a Padé-type iterative refinement to improve accuracy.

The proposed approximations and iterative methods, while relatively simple in form, demonstrated good convergence and accuracy for both real and complex inputs. Since the computational cost of evaluating square roots is comparable to that of division, the approaches presented here, though not widely established, may serve as practical options for cube root computation.

Further research into higher-order initial approximations and comparisons with alternative iterative methods may lead to even higher precision or faster algorithms in the future.

Note: The initial approximation method for cube roots described in this article — namely, constructing a Pade approximation of x ^2/3 near x = 1, then substituting x with √x to transform it into a function of x ^1/3, thereby expanding the effective approximation range — is a technique independently developed by the author based on the idea of Puiseux series. As of the time of writing, no similar approach has been found in the literature.

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.
"cbrt.c", GNU C Library — (accessed April 11, 2025, version 2.41) glibc-2.41
Jean-Michel Muller, "Elementary Functions: Algorithms and Implementation", Birkhauser, 3rd Edition, 2016.
Okumura, Haruhiko; Shudo, Kazuyuki; Sugiura, Masaki; Tsuchimura, Nobuyuki; Tsuru, Kazuo; Hosoda, Takayuki; Matsui, Yoshimitsu; & Mitsunari, Shigeo. (2003). Java ni yoru Arugorizumu Jiten [Encyclopedia of Algorithms with Java]. Gijutsu Hyoron Sha. ISBN4-7741-1729-3 [in Japanese]

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

External Site

"Cube root" — Wikipedia
"Padé approximant" — Wikipedia
"Puiseux series" — Wikipedia