Your solution has two problems, both related to the use of floating-point arithmetic. The first issue is that Python float
s only carry roughly 16 significant decimal digits of precision, so as soon as your answer requires more than 16 significant digits or so (so more than 6 digits before the point, and 10 digits after), you've very little hope of getting the correct trailing digits. The second issue is more subtle, and affects even small values of n
. That's that your approach of rounding to 11 decimal digits and then dropping the last digit suffers from potential errors due to double rounding. For an example, take n = 33
. The cube root of n
, to 20 decimal places or so, is:
3.20753432999582648755...
When that's rounded to 11 places after the point, you end up with
3.20753433000
and now dropping the last digit gives 3.2075343300
, which isn't what you wanted. The problem is that that round to 11 decimal places can end up affecting digits to the left of the 11th place digit.
So what can you do to fix this? Well, you can avoid floating-point altogether and reduce this to a pure integer problem. We need the cube root of some integer n
to 10 decimal places (rounding the last place down). That's equivalent to computing the cube root of 10**30 * n
to the nearest integer, again rounding down, then dividing the result by 10**10
. So the essential task here is to compute the floor of the cube root of any given integer n
. I was unable to find any existing Stack Overflow answers about computing integer cube roots (still less in Python), so I thought it worth showing how to do so in detail.
Computing cube roots of integers turns out to be quite easy (with the help of a tiny bit of mathematics). There are various possible approaches, but one approach that's both efficient and easy to implement is to use a pure-integer version of the Newton-Raphson method. Over the real numbers, Newton's method for solving the equation x**3 = n
takes an approximation x
to the cube root of n
, and iterates to return an improved approximation. The required iteration is:
x_next = (2*x + n/x**2)/3
In the real case, you'd repeat the iteration until you reached some desired tolerance. It turns out that over the integers, essentially the same iteration works, and with the right exit condition it will give us exactly the correct answer (no tolerance required). The iteration in the integer case is:
a_next = (2*a + n//a**2)//3
(Note the uses of the floor division operator //
in place of the usual true division operator /
above.) Mathematically, a_next
is exactly the floor of (2*a + n/a**2)/3
.
Here's some code based on this iteration:
def icbrt_v1(n, initial_guess=None):
"""
Given a positive integer n, find the floor of the cube root of n.
Args:
n : positive integer
initial_guess : positive integer, optional. If given, this is an
initial guess for the floor of the cube root. It must be greater
than or equal to floor(cube_root(n)).
Returns:
The floor of the cube root of n, as an integer.
"""
a = initial_guess if initial_guess is not None else n
while True:
d = n//a**2
if a <= d:
return a
a = (2*a + d)//3
And some example uses:
>>> icbrt_v1(100)
4
>>> icbrt_v1(1000000000)
1000
>>> large_int = 31415926535897932384626433
>>> icbrt_v1(large_int**3)
31415926535897932384626433
>>> icbrt_v1(large_int**3-1)
31415926535897932384626432
There are a couple of annoyances and inefficiencies in icbrt_v1
that we'll fix shortly. But first, a brief explanation of why the above code works. Note that we start with an initial guess that's assumed to be greater than or equal to the floor of the cube root. We'll show that this property is a loop invariant: every time we reach the top of the while loop, a
is at least floor(cbrt(n))
. Furthermore, each iteration produces a value of a
strictly smaller than the old one, so our iteration is guaranteed to eventually converge to floor(cbrt(n))
. To prove these facts, note that as we enter the while
loop, there are two possibilities:
Case 1. a
is strictly greater than the cube root of n
. Then a > n//a**2
, and the code proceeds to the next iteration. Write a_next = (2*a + n//a**2)//3
, then we have:
a_next >= floor(cbrt(n))
. This follows from the fact that (2*a + n/a**2)/3
is at least the cube root of n
, which in turn follows from the AM-GM inequality applied to a
, a
and n/a**2
: the geometric mean of these three quantities is exactly the cube root of n
, so the arithmetic mean must be at least the cube root of n
. So our loop invariant is preserved for the next iteration.
a_next < a
: since we're assuming that a
is larger than the cube root, n/a**2 < a
, and it follows that (2a + n/a**2) / 3
is smaller than a
, and hence that floor((2a + n/a**2) / 3) < a
. This guarantees that we make progress towards the solution at each iteration.
Case 2. a
is less than or equal to the cube root of n
. Then a <= floor(cbrt(n))
, but from the loop invariant established above we also know that a >= floor(cbrt(n))
. So we're done: a
is the value we're after. And the while loop exits at this point, since a <= n // a**2
.
There are a couple of issues with the code above. First, starting with an initial guess of n
is inefficient: the code will spend its first few iterations (roughly) dividing the current value of a
by 3
each time until it gets into the neighborhood of the solution. A better choice for the initial guess (and one that's easily computable in Python) is to use the first power of two that exceeds the cube root of n
.
initial_guess = 1 << -(-n.bit_length() // 3)
Even better, if n
is small enough to avoid overflow, is to use floating-point arithmetic to provide the initial guess, with something like:
initial_guess = int(round(n ** (1/3.)))
But this brings us to our second issue: the correctness of our algorithm requires that the initial guess is no smaller than the actual integer cube root, and as n
gets large we can't guarantee that for the float-based initial_guess
above (though for small enough n
, we can). Luckily, there's a very simple fix: for any positive integer a
, if we perform a single iteration we always end up with a value that's at least floor(cbrt(a))
(using the same AM-GM argument that we used above). So all we have to do is perform at least one iteration before we start testing for convergence.
With that in mind, here's a more efficient version of the above code:
def icbrt(n):
"""
Given a positive integer n, find the floor of the cube root of n.
Args:
n : positive integer
Returns:
The floor of the cube root of n, as an integer.
"""
if n.bit_length() < 1024: # float(n) safe from overflow
a = int(round(n**(1/3.)))
a = (2*a + n//a**2)//3 # Ensure a >= floor(cbrt(n)).
else:
a = 1 << -(-n.bit_length()//3)
while True:
d = n//a**2
if a <= d:
return a
a = (2*a + d)//3
And with icbrt
in hand, it's easy to put everything together to compute cube roots to ten decimal places. Here, for simplicity, I output the result as a string, but you could just as easily construct a Decimal
instance.
def cbrt_to_ten_places(n):
"""
Compute the cube root of `n`, truncated to ten decimal places.
Returns the answer as a string.
"""
a = icbrt(n * 10**30)
q, r = divmod(a, 10**10)
return "{}.{:010d}".format(q, r)
Example outputs:
>>> cbrt_to_ten_places(2)
'1.2599210498'
>>> cbrt_to_ten_places(8)
'2.0000000000'
>>> cbrt_to_ten_places(31415926535897932384626433)
'315536756.9301821867'
>>> cbrt_to_ten_places(31415926535897932384626433**3)
'31415926535897932384626433.0000000000'