SICP-RS: Exercise 1.8 - Cube Roots

Part 3 of the SICP-RS 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 cube-root procedure analogous to the square-root 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/sicp-rs$ 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://rust-lang.github.io/rust-clippy/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 built-in method on floats I was unaware of, mul_add. According to the docs:

Fused multiply-add. Computes (self * a) + b with only one rounding error, yielding a more accurate result than an unfused multiply-add.
Using mul_add may be more performant than an unfused multiply-add if the target architecture has a dedicated fma 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://rust-lang.github.io/rust-clippy/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.rust-lang.org/std/primitive.f64.html#method.cbrt
            assert!((n.cbrt() - cbrt(n)).abs() < 1e-10);
        }
    }
}

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 built-in 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 (cbrt-iter guess x)
    (if (good-enough? guess x)
        guess
        (cbrt-iter (improve guess x) x)))

(define (improve guess x)
    (/ (+ (/ x (* guess guess))
          (* 2 guess))
        3))

(define (good-enough? guess x)
    (= (improve guess x) guess))

(define (cbrt x)
    (cbrt-iter 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 low-level 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 cube-root procedure analogous to the
//! square-root 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://rust-lang.github.io/rust-clippy/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.rust-lang.org/std/primitive.f64.html#method.cbrt
            assert!((n.cbrt() - cbrt(n)).abs() < 1e-10);
        }
    }
}