Why does my iterative higher-order procedure give more precise results than my equivalent recursive procedure?

Question

Exercise 1.30 of SICP invites us to rewrite the following code as an iterative (read: tail recursive) procedure:

(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

To assist us, we are given the following template:

(define (sum term a next b)
  (define (iter a result)
    (if <??>
        <??>
        (iter <??> <??>)))
  (iter <??> <??>))

After one false start, I produced the following answer, which appears correct:

(define (sumIT term a next b)
  (define (iter a result)
    (if (> a b)
        result
        (iter (next a) (+ (term a) result))))
  (iter a 0))

To test this, I copied the following code from just above where this exercise is given:

(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))
(integral cube 0 1 0.01)
(integral cube 0 1 0.001) 

And quickly made this version using sumIT:

(define (integralIT f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sumIT f (+ a (/ dx 2.0)) add-dx b)
     dx))
(integralIT cube 0 1 0.01)
(integralIT cube 0 1 0.001) 

Then, as I am running in DrRacket's #lang sicp mode, which does not have a cube procedure by default, I also defined:

(define (cube x) (* x x x))

Finally, I ran the code.

Unsurprisingly, (integral cube 0 1 0.01) and (integralIT cube 0 1 0.01) both produce the exact same result: 0.24998750000000042. But, to my shock and horror, (integral cube 0 1 0.001) returns 0.249999875000001 whereas (integralIT cube 0 1 0.001) returns a more precise answer of 0.24999987500000073.

Why is this? Why would rewriting a recursive higher-order procedure to be tail recursive increase the precision of my result? I can't think of any part or error in my code that would cause this.

Answer 1

Floating-point addition and multiplication are commutative, but not associative. This means (+ (+ a b) c) might not be equal to (+ a (+ b c), due to error propagation. As alredy mentioned by Óscar López, see What Every Computer Scientist Should Know About Floating-Point Arithmetic .

I rewrote your example in Common Lisp, and I can observe the same behaviour.

Then I changed a little bit the code by inserting backquotes and commas, so that the functions returns a tree instead of a number. See Macros: defining your own if you are not familiar with the synax.

The tree is made of numbers and symbols (*, +, ...) and corresponds to the operations that are made in each version, as written. This should help understand how and when intermediate results are computed.

Recursive sum:

(defun sum (term a next b)
  (labels ((it (a)
             (if (> a b)
                 0
                 `(+ ,(it (funcall next a))
                     ,(funcall term a)))))
    (it a)))

Tail-recursive sum:

(defun sum-it (term a next b)
  (labels ((it (a result)
             (if (> a b)
                 result
                 (it (funcall next a)
                     `(+ ,(funcall term a) ,result)))))
    (it a 0)))

Recursive integral:

(defun integral (f a b dx)
  (flet ((add-dx (x) (+ x dx)))
    `(* ,(sum f (+ a (/ dx 2.0d0)) #'add-dx b) ,dx)))

Tail-recursive integral:

(defun integral-it (f a b dx)
  (flet ((add-dx (x) (+ x dx)))
    `(* ,(sum-it f (+ a (/ dx 2.0d0)) #'add-dx b) ,dx)))

And cube:

(defun cube (x) (* x x x))

Take care when testing: with the target precision argument being too precise, you'll have a large result.

You can see that both approaches eventually computes the result in a different order:

CL-USER> (integral #'cube 0 1 0.4d0)
(* (+ (+ (+ 0 1.0d0) 0.21600000000000008d0) 0.008000000000000002d0) 0.4d0)

Compared to:

CL-USER> (integral-it #'cube 0 1 0.4d0)
(* (+ 1.0d0 (+ 0.21600000000000008d0 (+ 0.008000000000000002d0 0))) 0.4d0)

In the above example, the result does not differ, but for some inputs like the one you found, there will be a difference.