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.

"_... returns a more precise answer...._" -- Racket uses double-precision IEEE floating-point numbers for inexact reals by default; you are bumping up against the limits of precision, which is about 15 significant figures for doubles. Welcome to the world of numeric computation ;) — ad absurdum, Nov 15 '20 at 19:11
There are no errors in your code, it's just the way floating-point arithmetic works. Mandatory reading: [What Every Computer Scientist Should Know About Floating-Point Arithmetic](https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) — Óscar López, Nov 15 '20 at 19:27
0.249999875000001 is not less precise than 0.24999987500000073. Both are just printed representations of some IEEE 64 bit double. That printed represenation happens to chop off trailing zeros, so that 0.2499998750000010 becomes 0.249999875000001. What you have is simply two *different* results, due to different calculation paths showing different error accumulation behavior. — Kaz, Nov 16 '20 at 02:39
@ÓscarLópez Is that still all mandatory these days? Most of it is about memory details and the IEEE standard. Should I be paying attention to a particular section? — J. Mini, Jan 14 '21 at 22:22

coredump · Accepted Answer · 2020-11-16T12:48:30.080

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.

Is there any obvious reason why the recursive version adds the terms in a different order than the iterative version? — J. Mini, Nov 15 '20 at 22:07
The sum is done on the accumulator in the iterative version. I mean, what I want to demonstrate by showing the resulting source code is that the way it is written makes the operation happen in a different order. — coredump, Nov 15 '20 at 22:11
Oh, that was even more obvious than I was expecting. I was foolish not to see why the order of operations is different. — J. Mini, Nov 15 '20 at 22:13
Is all of _What Every Computer Scientist Should Know About Floating-Point Arithmetic_ still considered important reading? Most of it is about memory details and the IEEE standard. Should I be paying attention to a particular section? — J. Mini, Jan 14 '21 at 22:23
I'd recommend having at least a quick look at each section, not necessarily in details, to see the different kind of issues associated with floats. In particular The Compiler and Languages section, and the conclusion (Summary). You don't need to know exactly how it works but know that there's a standard float format used by everyone which has some gotchas wrt precision, rounding errors etc. The article is quite exhaustive but you can bookmark it for when you really need to dive more deeply — coredump, Jan 14 '21 at 23:49

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

1 Answers1