SICP-RS: Exercise 1.7 - Square Roots, Newton Style

Part 2 of the SICP-RS Series

The next coding exercise in Chapter 1 of Structure and Interpretation of Computer Programs is Exercise 1.7:

The good-enough? test used in computing square roots will not be very effective for finding the square roots of very small numbers. Also, in real computers, arithmetic operations are almost always performed with limited precision. This makes our test inadequate for very large numbers. Explain these statements, with examples showing how the test fails for small and large numbers. An alternative strategy for implementing good-enough? is to watch how guess changes from one iteration to the next and to stop when the change is a very small fraction of the guess. Design a square-root procedure that uses this kind of end test. Does this work better for small and large numbers? ^[1]

In the chapter, there is an implementation of Newton's method for calculating square roots. But the implementation provided isn't adequate for very large and small numbers. The task is to create a new version that tracks the difference of each guess to know when to stop.

Step 1: Comparing Floats

A key piece of the implementation will be knowing when to stop. If the function can't improve the guess, that will be the closest approximation of the square root we can calculate.

Rust has two float types, f32 (32 bit) and f64 (64 bit), each with its own level of precision of what numbers it can represent. So, to compare if the guess improved, we need a way to check if the delta between two guesses is smaller than that precision can represent.

/// Compare if two float values are equal (or as close as we can get)
fn f64_eq(x: f64, y: f64) -> bool {
    // Check if they are equal, or as close as we can get with floating point precision
    x == y || (x - y).abs() <= f64::EPSILON
}

We first check if we can get a straight equality check for free. If not, then check if the absolute difference between the two values is less than the machine epsilon.

Step 2: Improving our Guess

According to Newton's method, the way to improve our guess is to average the guess with the original number divided by the guess.

So let's write an improve function, along with a few tests to make sure it is doing what we expect:

/// Improve square root guess by averaging the guess
/// with the number (x) divided by the current guess
fn improve(guess: f64, x: f64) -> f64 {
    (guess + x / guess) / 2.0
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn test_improve() {
        assert!(f64_eq(improve(1.0, 2.0), 1.5));
        assert!(f64_eq(improve(1.5, 2.0), 1.416_666_666_666_666_5));
        assert!(f64_eq(
            improve(1.416_666_666_666_666_5, 2.0),
            1.414_215_686_274_509_7
        ));
    }
}

We use our f64_eq function here in the tests so that we can check if we are getting the float values we expect. Including some small numbers, since that is part of the goal of this exercise.

Step 3: Square Root

With these helper functions in place, we can write our square root calculation.

First, we should have some tests to assert the behavior we expect. In this case, we have a great, built-in way to see if we are getting close. The Rust Standard Library already provides a sqrt method on f64, so let's see if our method can provide the same value! I'm sure the included method is more performant, but we can at least see if we are generating the same values.

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_sqrt() {
        let numbers: [f64; 5] = [
            9.0,
            0.0001,
            10_000_000_000_000.000_1,
            100_000_000_000_000_000_000.0,
            0.000_000_000_000_1,
        ];

        for n in numbers {
            assert!(f64_eq(sqrt(n), n.sqrt()));
        }
    }
}

There is an array of values to test, some very small and some very large. The test then loops over them and makes the same assertion for each.

Now to implement sqrt ourselves!

/// Find the square root of a given number using Newton's method
pub fn sqrt<T>(x: T) -> f64
where
    f64: From<T>,
{
    // Convert x into a f64
    let x = f64::from(x);

    // Initial guess
    let mut guess = 1.0;

    // 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
}

This allows the user to pass in any value, so long as the value can convert into an f64 value. We then cache x as an f64 value so that we don't have to convert it on every iteration.

The guess begins at 1.0, and then the function enters into a loop. On each loop iteration we improve the guess, and then compare the new guess to the current one. If they are equal, we are at the limit of what our precision can calculate. We are as close as we are going to get, and break out of the loop.

Otherwise, we assign the new guess to the guess variable and continue onto the next iteration.

Once outside the loop, the function returns the final guess.

Conclusion

This sqrt method can return the same values as the one already provided by Rust. There isn't much use in reusing this function in the future, but it was still helpful to implement it ourselves.

Since SICP is dealing with Scheme, let's also compare this Rust implementation to the Scheme version:

(define (average x y) (/ (+ x y) 2))

(define (improve guess x)
    (average guess (/ x guess)))

(define (good-enough? guess x)
    (= (improve guess x) guess))
    
(define (sqrt-iter guess x)
    (if (good-enough? guess x)
        guess
        (sqrt-iter (improve guess x) x)))

(define (sqrt x)
    (sqrt-iter 1.0 x))

It isn't too different. The biggest difference though is that instead of a loop, this implementation uses recursion to call the same method, sqrt-iter, over and over. And this is a core difference between the two languages.

Scheme uses recursive functions for iterative behavior. Rust can call functions recursively, but because it is lacking tail call optimizations, this can run into performance limitations (see my very bad fibonacci solution for an example).

Thanks to Rust's borrow checker and ownership model though, it is quite easy, and safe, to mutate our existing guess without fear and use a loop. Both methods have a very similar control flow anyway. We continue to call the same code, and break out once we reach our desired state.

The full Rust implementation is below if you would like to see it all together. You can also view it on my github repository, along with the Scheme code as well.

//! Exercise 1.7: The good-enough? test used in computing square roots will not
//! be very effective for finding the square roots of very small numbers. Also,
//! in real computers, arithmetic operations are almost always performed with
//! limited precision. This makes our test inadequate for very large numbers.
//! Explain these statements, with examples showing how the test fails for
//! small and large numbers. An alternative strategy for implementing
//! good-enough? is to watch how guess changes from one iteration to the next
//! and to stop when the change is a very small fraction of the guess. Design a
//! square-root procedure that uses this kind of end test. Does this work
//! better for small and large numbers?

/// Compare if two float values are equal (or as close as we can get)
fn f64_eq(x: f64, y: f64) -> bool {
    // Check if they are equal, or as close as we can get with floating point precision
    x == y || (x - y).abs() <= f64::EPSILON
}

/// Improve square root guess by averaging the guess
/// with the number (x) divided by the current guess
fn improve(guess: f64, x: f64) -> f64 {
    (guess + x / guess) / 2.0
}

/// Find the square root of a given number using Newton's method
pub fn sqrt<T>(x: T) -> f64
where
    f64: From<T>,
{
    // Convert x into a f64
    let x = f64::from(x);

    // Initial guess
    let mut guess = 1.0;

    // 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_improve() {
        assert!(f64_eq(improve(1.0, 2.0), 1.5));
        assert!(f64_eq(improve(1.5, 2.0), 1.416_666_666_666_666_5));
        assert!(f64_eq(
            improve(1.416_666_666_666_666_5, 2.0),
            1.414_215_686_274_509_7
        ));
    }

    #[test]
    fn test_sqrt() {
        let numbers: [f64; 5] = [
            9.0,
            0.0001,
            10_000_000_000_000.000_1,
            100_000_000_000_000_000_000.0,
            0.000_000_000_000_1,
        ];

        for n in numbers {
            assert!(f64_eq(sqrt(n), n.sqrt()));
        }
    }
}

Thanks to the scheme community wiki for some great test case examples. Ensuring my rust implementation worked correctly was much easier with a standard test set.