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 xm / 2n, to ensure stable convergence across magnitudes.
The square root of x is given by x1/2 and The cube root of x is given by x1/3. When approximating the exponent 1/3 using the form m / 2n, where m = round(2n / 3), the approximation error is given by:
m / 2n × 3 - 1which results in an error of ±1/2n (with the negative sign when n is even). The specific values are:
Among these, it was found that x11/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].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/2n![]()
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.
Using the initial approximation obtained from equation (1), we refine the cube root using a recurrence relation similar to Newton's method.
Defining the relative error h as:
The cube root function x1/3 can be expanded into an infinite series for x ≠ 0. Thus, we can express (1 + h)1/3 as the following series expansion:
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:
From equations (3) and (5), we derive the following recurrence relation:![]()
![]()
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.
![]()
![]()
Fig.2 cbrt Initial approximatin error and its convergence behavior
(Computed under IEEE 754 Double precision environment)
Formulas usedInitial approximation :Recurrence formula of third order convergence :
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)
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]:
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 × 32 = 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 × 32 ≈ 20.79 digits > 19 digits (extended precision)and hence allows the desired accuracy to be achieved.
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)3 / 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.
Usage of cbrt (HP 42S)OperationCalculation example of cbrt (HP 42S)Input
[XEQ] "cbrt" — cube root function
Output
REG X : argument value a
REG y : a REG X : cbrt(a)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