Search code examples
schemesicpmit-scheme

`1.0+1e-100=1.` in MIT-scheme

I am reading SICP and read this recitation about chapter 1. I use MIT/GNU Scheme as MIT course 6.5151 (6.905) does with version 12.1.

I am stuck at Problem 3.

Write a procedure that computes e.

I tried to use . So I have code

(define (e x)
  (newline)
  (display (/ 1.0 x))
  (newline)
  (display (+ 1.0 (/ 1.0 x)))
  (newline)
  (display x)
  (expt (+ 1.0 (/ 1.0 x)) x))
(e (expt 10 100))

It has answer 1.. Here I use display to check where the resolution is not as expected. It turns out (display (+ 1.0 (/ 1.0 x))) outputs 1. which is unexpected.

The instructor solution uses Taylor series which is different from the above idea.

Solution

  • The problem is that 1 + 10-100 differs from 1 in the 100th decimal place. That means that, as a float, in order for it to be distinct from 1 you need at least 100 decimal places, or at least about 332 bits of precision. IEEE 754 double-precision floats have a precision of 53 bits, so there is no way that this can work using double floats.

    In general if you're adding numbers of hugely differing magnitudes like this you are going to run into trouble with floating-point, or indeed with any kind of machine integer: if you tried to represent these things as 64-bit integers you still can't do anything useful here.

    Using Racket's 'bigfloat' package you can control the precision of floats and see the answer appear.

    (require math/bigfloat)

(define (bfe x (precision (bf-precision)))
  (define bfx (bf x))
  (parameterize ([bf-precision precision])
    (bfexpt (bf+ 1.bf (bf/ 1.bf bfx)) bfx)))

    Now:

    > (bf-precision)
128
> (bfe (expt 10 100))
(bf 1)
> (bfe (expt 10 100) 332)
(bf 1)
> (bfe (expt 10 100) 333)
(bf #e3.13612321207875129576038816881464108475667088571290143008146478665151883781810987619137532041768112489)
> (bfe (expt 10 100) 400)
(bf #e2.7182818284590452353645832758305382066800249084717283867982178383777754898191243109391113320627898418259308185781451392861)

    So you can see that to get anything other than 1 you need 333 bits of precision and to get any reasonable answer you need many more bits than that.