0

This code is taken from Sussman and Wisdom's Structure and Interpretation of Classical Mechanics, its purpose is to derive (close to) the smallest positive floating point the host machine supports. https://github.com/hnarayanan/sicm/blob/e37f011db68f8efc51ae309cd61bf497b90970da/scmutils/src/kernel/numeric.scm

Running it in DrRacket results in 2.220446049250313e-016 on my machine.

My question, what causes this to even return a value? This code is tail recursive, and it makes sense at some point the computer can no longer divide by 2. Why does it not throw?

(define *machine-epsilon*
  (let loop ((e 1.0))
     (if (= 1.0 (+ e 1.0))
         (* 2 e)
         (loop (/ e 2)))))
*machine-epsilon*
Jack Fox
  • 953
  • 6
  • 20

4 Answers4

4

This code is tail recursive, and it makes sense at some point the computer can no longer divide by 2. Why does it not throw?

No, the idea is different: at some point the computer still can divide by 2, but the result (e) becomes indistinguishable from 0 [upd: in the context of floating-point addition only - very good point mentioned in the comment] (e + 1.0 = 1.0, this is exactly what if clause is checking). We know for sure that the previous e was still greater than zero "from the machine point of view" (otherwise we wouldn't get to the current execution point), so we simply return e*2.

Konstantin Strukov
  • 2,899
  • 1
  • 10
  • 14
  • 1
    `e` does not become indistinguishable from `0`: `(+ 1.0 e)` becomes indistinguishable from `1.0`, which is a very different thing. –  Jan 22 '20 at 17:58
  • @tfb the question was not "why if clause works", right? :) But seriously speaking you are absolutely correct - `(= 1.0 (+ e 1.0))` doesn't imply `(= e 0)` – Konstantin Strukov Jan 22 '20 at 20:19
1

It becomes clearer when you get rid of the confusing named let notation.

(define (calculate-epsilon (epsilon 1.0))
  (if (= 1.0 (+ 1.0 epsilon))
      (* epsilon 2)
      (calculate-epsilon (/ epsilon 2))))

(define *machine-epsilon* (calculate-epsilon))

Is what the code does actually. So now we see for what the named let expression is good. It defines locally the function and runs it. Just that the name of the function as loop was very imprecise and confusing and the naming of epsilon to e is a very unhappy choice. Naming is the most important thing for readable code.

So this example of SICP should be an example for bad naming choices. (Okay, maybe they did it by intention to train the students).

The named let defines and calls/runs a function/procedure. Avoiding it would lead to better code - since clearer.

In common lisp such a construct would be much clearer expressed:

(defparameter *machine-epsilon*
  (labels ((calculate-epsilon (&optional (epsilon 1.0))
              (if (= 1.0 (+ 1.0 epsilon))
                  (* epsilon 2)
                  (calculate-epsilon (/ epsilon 2)))))
    (calculate-epsilon)))

In CLISP implementation, this gives: 1.1920929E-7

Gwang-Jin Kim
  • 9,303
  • 17
  • 30
  • `*machine-epsilon*` as defined in the question *is* a number, and it's also a global variable. –  Jan 22 '20 at 17:26
  • @tfb yes, `(define *machine-epsilon* ... )` is a variable definition. What I meant was that the named `(let ...)` construct is obscuring two things: a function definition and its call. But true, at the beginning when I wrote, I was confused. - Okay I deleted a big part. – Gwang-Jin Kim Jan 22 '20 at 22:10
1

This form of let-binding is syntactic sugar for recursion.

You may avoid using too much syntax until you master the language and write as much as possible using the kernel language, to focus on essential problem. For example, in full SICP text, never is specified this syntactic sugar for iteration.

The r6rs definition for iteration is here.

alinsoar
  • 15,386
  • 4
  • 57
  • 74
1

The purpose of this code is not to find the smallest float that the machine can support: it is to find the smallest float, epsilon such that (= (+ 1.0 epsilon) 1.0) is false. This number is useful because it's the upper bound on the error you get from adding numbers In particular what you know is that, say, (+ x y) is in the range [(x+y)*(1 - epsilon), (x+y)*(1 + epsilon)], where in the second expression + &c mean the ideal operations on numbers.

In particular (/ *machine-epsilon* 2) is a perfectly fine number, as is (/ *machine-epsilon* 10000) for instance, and (* (/ *machine-epsilon* x) x) will be very close to *machine-epsilon* for many reasonable values of x. It's just the case that (= (+ (/ *machine-epsilon* 2) 1.0) 1.0) is true.

I'm not familiar enough with floating-point standards, but the number you are probably thinking of is what Common Lisp calls least-positive-double-float (or its variants). In Racket you can derive some approximation to this by

(define *least-positive-mumble-float*
  ;; I don't know what float types Racket has if it even has more than one.
  (let loop ([t 1.0])
    (if (= (/ t 2) 0.0)
        t
        (loop (/ t 2)))))

I am not sure if this is allowed to raise an exception: it does not in practice and it gets a reasonable-looking answer.