Why does Math.cbrt(1728) produce a more accurate result than Math.pow(1728, 1/3)?

Question

In JavaScript, Math.cbrt(1728) evaluate to the exact result of 12.

However, the seemingly-equivalent expression Math.pow(1728, 1/3) evaluates to 11.999999999999998.

Why do these results vary in precision?

Answer 1

A couple of general remarks up front:

1. As explained in this seminal paper, due to finite precision and range limits, floating-point arithmetic is sufficiently different from real mathematics (for example, lack of associativity) that mathematically equivalent expressions are not necessarily equivalent when evaluated in floating-point arithmetic.

2. Standards for computer languages do not typically guarantee any particular accuracy for math functions, or identical error bounds between different math functions such as cbrt() or pow(). But math libraries that deliver correctly rounded results for a given precision do exist, such as CRlibm.

In this case however, cbrt(x) will deliver more accurate results than pow(x,1.0/3.0) even when both functions are correctly rounded for all inputs.

The issue is that 1.0/3.0 cannot be represented exactly as a floating-point number, whether in binary or decimal. The IEEE-754 double precision number closest to one third is 3.3333333333333331e-1 (or 0x1.5555555555555p-2 when expressed in the C/C++ hexadecimal floating-point format). The relative representational error is -5.5511151231257827e-17 (-0x1.0000000000000p-54), meaning the best double-precision representation of 1/3 is somewhat smaller than the desired mathematical value.

This initial error in one of the inputs of pow() is not only passed through to the output, it is magnified due to the error magnification property of exponentiation. As a result, pow(x,1.0/3.0) will generally deliver results that are too small compared to the desired cube root, even if pow() delivers correctly rounded results. For the example in the question, the correctly rounded results are

cbrt(1728.0)        = 1.2000000000000000e+1  (0x1.8000000000000p+3)
pow(1728.0,1.0/3.0) = 1.1999999999999998e+1  (0x1.7ffffffffffffp+3)

that is, the result from pow() is one ulp smaller than the result from cbrt(). For arguments large in magnitude, the difference will be much larger. For example, if x is 2¹⁰²², the respective results differ by 94 ulps:

x              = 4.4942328371557898e+307  (0x1.0000000000000p+1022)
cbrt(x)        = 3.5553731598732904e+102  (0x1.965fea53d6e3dp+340)
pow(x,1.0/3.0) = 3.5553731598732436e+102  (0x1.965fea53d6ddfp+340)

The relative error in the result of pow() in this example is 1.3108e-14, demonstrating the magnification of the relative error mentioned above.

For reasons of both accuracy and performance, math libraries that implement cbrt() therefore typically do not map cbrt(x) to pow(x,1.0/3.0) but use alternative computational schemes. While implementations will differ, a commonly used approach is to start with an initial low-precision approximation followed by one or several steps of Halley's method which has cubic convergence.

As a rule of thumb, when a computer language offers both a dedicated cube root functionality and general exponentiation functionality, the former should be preferred to the latter for the computation of cube roots.