SICPRS: Exercise 1.8  Cube Roots
Part 3 of the SICPRS Series
The next coding exercise in Chapter 1 of Structure and Interpretation of Computer Programs is Exercise 1.8, which is a slight variation on the 1.7 exercise:
Exercise 1.8: Newton’s method for cube roots is based on the fact that if y is an approximation to the cube root of x, then a better approximation is given by the value (x / y^2 + 2y) / 3. Use this formula to implement a cuberoot procedure analogous to the squareroot procedure. ^{[1]}
Step 1: New improve
Function
The basic calculation is almost identical to the last one, with a minor change to how we improve our guess.
For a cube root, we have the following equation:
(x / y^2 + 2y) / 3
So let's implement it:
use crate::ex1_3::square;
/// Improve cube root guess using Newton's method for approximating a better guess
fn improve(guess: f64, x: f64) > f64 {
(x / square(guess) + 2.0 * guess) / 3.0
}
It is mostly a translation of the equation into Rust syntax. Notice, that there is also a use
statement to bring in our square
function from exercise 1.3.
After writing this though, I got a helpful warning from Clippy, Rust's linter:
~/dev/sicprs$ cargo clippy
warning: multiply and add expressions can be calculated more efficiently and accurately
> ch1/src/ex1_8.rs:14:5

14  (x / square(guess) + 2.0 * guess) / 3.0
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ help: consider using: `2.0f64.mul_add(guess, x / square(guess))`
= help: for further information visit https://rustlang.github.io/rustclippy/master/index.html#suboptimal_flops
warning: `ch1` (lib) generated 1 warning
This is one of the reasons I love writing code in Rust, and why I especially like using Clippy! Both Clippy and the compiler often look for ways you can improve your code, which helps you discover better ways to use the Rust language.
In this case, it looks like there is a builtin method on floats I was unaware of, mul_add
. According to the docs:
Fused multiplyadd. Computes
(self * a) + b
with only one rounding error, yielding a more accurate result than an unfused multiplyadd.
Usingmul_add
may be more performant than an unfused multiplyadd if the target architecture has a dedicatedfma
CPU instruction. However, this is not always true, and will be heavily dependant on designing algorithms with specific target hardware in mind.
So Clippy is letting us know there is a more accurate and possibly more performant way to achieve the same result. Here's the improve
function of this new approach. Not that we are annotating 2.0
with an f64
annotation so the compiler knows which type we are using for this method.
/// Improve cube root guess using Newton's method for approximating a better guess
fn improve(guess: f64, x: f64) > f64 {
// Should be the same as: (x / square(guess) + 2.0 * guess) / 3.0
// https://rustlang.github.io/rustclippy/master/index.html#suboptimal_flops
2.0f64.mul_add(guess, x / square(guess)) / 3.0
}
Step 2: cbrt
To write the cbrt
function, we can use sqrt
as a starting point and make a few changes.
use crate::ex1_7::f64_eq;
/// Find the cube root of a given number using Newton's method
pub fn cbrt<T>(x: T) > f64
where
f64: From<T>,
{
// Convert x into a f64
let x = f64::from(x);
// Initial guess
let mut guess = 1.1;
// improve the guess until it is good enough with our precision
loop {
let next_guess = improve(guess, x);
// If these two are equal, break out of the loop, we are done.
if f64_eq(next_guess, guess) {
break;
}
// Otherwise update guess and try again.
guess = next_guess;
}
guess
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_sqrt() {
let numbers: [f64; 5] = [5.0, 2.0, 27.0, 0.0, 100_000_000_000_000.000_1];
for n in numbers {
// Level of precision used in docs for cbrt()
// https://doc.rustlang.org/std/primitive.f64.html#method.cbrt
assert!((n.cbrt()  cbrt(n)).abs() < 1e10);
}
}
}
It is almost identical. The only difference is the name and the initial guess. When I wrote the tests, it seemed to have an issue calculating the cube root of 2.0
if the initial guess was 1.0
, so I increased it to 1.1
.
Also, note the difference in the tests. When comparing our implementation with the builtin version, the precision wasn't quite the same, so f64_eq
was too strict an assertion for this implementation. In the cbrt()
docs, this seems to be expected as they only check that it is accurate to within 1e10
, so I used the same assertion here.
Conclusion
It was nice to be able to leverage existing code for this exercise. This exercise also gave us a chance to see how Clippy can be helpful in our implementations and give us hints on how to improve our code.
As always, here is the comparison with the Scheme version. As with 1.7, this approach uses recursion instead of a loop:
(define (cbrtiter guess x)
(if (goodenough? guess x)
guess
(cbrtiter (improve guess x) x)))
(define (improve guess x)
(/ (+ (/ x (* guess guess))
(* 2 guess))
3))
(define (goodenough? guess x)
(= (improve guess x) guess))
(define (cbrt x)
(cbrtiter 1.1 x))
The minor changes were almost identical in both Scheme and Rust, so we can derive mostly the same conclusion as the last exercise.
However, it is also apparent that Rust is concerned much more with performance, and helping programmers find more performant solutions. In this example, Rust provides us with a way to use lowlevel arithmetic operations to make sure our implementation is as performant and accurate as the hardware allows. Which is a huge win!
The full Rust implementation is below if you would like to see it all together. You can view it on my GitHub repository, along with the Scheme code as well.
//! Exercise 1.8: Newton’s method for cube roots is based on the fact that if y
//! is an approximation to the cube root of x, then a better approximation is
//! given by the value
//!
//! (x / y^2 + 2y) / 3
//!
//! Use this formula to implement a cuberoot procedure analogous to the
//! squareroot procedure.
use crate::{ex1_3::square, ex1_7::f64_eq};
/// Improve cube root guess using Newton's method for approximating a better guess
fn improve(guess: f64, x: f64) > f64 {
// Should be the same as: (x / square(guess) + 2.0 * guess) / 3.0
// https://rustlang.github.io/rustclippy/master/index.html#suboptimal_flops
2.0f64.mul_add(guess, x / square(guess)) / 3.0
}
/// Find the cube root of a given number using Newton's method
pub fn cbrt<T>(x: T) > f64
where
f64: From<T>,
{
// Convert x into a f64
let x = f64::from(x);
// Initial guess
let mut guess = 1.1;
// improve the guess until it is good enough with our precision
loop {
let next_guess = improve(guess, x);
// If these two are equal, break out of the loop, we are done.
if f64_eq(next_guess, guess) {
break;
}
// Otherwise update guess and try again.
guess = next_guess;
}
guess
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_sqrt() {
let numbers: [f64; 5] = [5.0, 2.0, 27.0, 0.0, 100_000_000_000_000.000_1];
for n in numbers {
// Level of precision used in docs for cbrt()
// https://doc.rustlang.org/std/primitive.f64.html#method.cbrt
assert!((n.cbrt()  cbrt(n)).abs() < 1e10);
}
}
}
Hal Abelson's, Jerry Sussman's, and Julie Sussman's Structure and Interpretation of Computer Programs, Second Edition (MIT Press, 1996; ISBN 0262510871) ↩︎