Pathological floating point problems

[cameyo]

From https://rosettacode.org/wiki/Pathological_floating_point_problems">//https://rosettacode.org/wiki/Pathological_floating_point_problems

Problem

The Chaotic Bank Society is offering a new investment account to their customers.

You first deposit $ (e - 1) where "e" is 2.7182818 (the base of natural logarithms).

After each year, your account balance will be multiplied by the number of years that have passed, and $1 in service charges will be removed.

Example:

after 1 year, your balance will be multiplied by 1 and $1 will be removed for service charges.

after 2 years your balance will be doubled and $1 removed.

after 3 years your balance will be tripled and $1 removed.

...

after 10 years, multiplied by 10 and $1 removed, and so on.

What will your balance be after 25 years?

The correct result is:

   Starting balance: $ (e - 1)

   Balance = (Balance * year) - 1 for 25 years

   Balance after 25 years: $ 0.0399387296732302

Solution (using floating point math):

Code Select Expand


(define (bank)
  (local (e balance)
    ;definiamo il numero e
    (setq e (exp 1))
    (setq balance (sub e 1))
    (for (i 1 25)
      (setq balance (sub (mul balance i) 1))
      (println i { } balance)
    )
    balance
  )
)

(bank)
;-> 1 0.7182818284590451
;-> 2 0.4365636569180902
;-> 3 0.3096909707542705
;-> ...
;-> 16 0.05924783418595325
;-> 17 0.007213181161205284
;-> 18 -0.8701627390983049
;-> 19 -17.53309204286779
;-> 20 -351.6618408573559
;-> 21 -7385.898658004473
;-> 22 -162490.7704760984
;-> 23 -3737288.720950264
;-> 24 -89694930.30280632
;-> 25 -2242373258.570158
;-> -2242373258.570158

The result is wrong, as the rounding of the floating point operations diverge the calculations.

To solve the problem we can use fractions, ie we perform all the calculations with fractions (integers) and we use the division only to get the value of the result as a floating point. To do this we must also represent the number "e" with a fraction:

Code Select Expand

e = 106246577894593683/39085931702241241
The functions to use the four operations of the fractions are the following:

Code Select Expand


(define (semplifica frac)
  (local (num den n d temp, nums dens)
    (setq num (first frac))
    (setq den (last frac))
    (setq n (first frac))
    (setq d (last frac))
    ; calcola il numero massimo che divide esattamente numeratore e denominatore
    (while (!= d 0)
      (setq temp d)
      (setq d (% n temp))
      (setq n temp)
    )
    (setq nums (/ num n))
    (setq dens (/ den n))
    ; controllo del segno
    (cond ((or (and (< dens 0) (< nums 0)) (and (< dens 0) (> nums 0)))
           (setq nums (* nums -1))
           (setq dens (* dens -1))
          )
    )
    (list nums dens)
  )
)

(define (+f frac1 frac2 redux)
  (local (num den n1 d1 n2 d2)
    (setq n1 (first frac1))
    (setq d1 (last frac1))
    (setq n2 (first frac2))
    (setq d2 (last frac2))
    (setq num (+ (* n1 d2) (* n2 d1)))
    (setq den (* d1 d2))
    (if redux (list num den)
          (semplifica (list num den))
    )
  )
)

(define (-f frac1 frac2 redux)
  (local (num den n1 d1 n2 d2)
    (setq n1 (first frac1))
    (setq d1 (last frac1))
    (setq n2 (first frac2))
    (setq d2 (last frac2))
    (setq num (- (* n1 d2) (* n2 d1)))
    (setq den (* d1 d2))
    (if redux (list num den)
          (semplifica (list num den))
    )
  )
)

(define (*f frac1 frac2 redux)
  (local (num den n1 d1 n2 d2)
    (setq n1 (first frac1))
    (setq d1 (last frac1))
    (setq n2 (first frac2))
    (setq d2 (last frac2))
    (setq num (* n1 n2))
    (setq den (* d1 d2))
    (if redux (list num den)
          (semplifica (list num den))
    )
  )
)

(define (/f frac1 frac2 redux)
  (local (num den n1 d1 n2 d2)
    (setq n1 (first frac1))
    (setq d1 (last frac1))
    (setq n2 (first frac2))
    (setq d2 (last frac2))
    (setq num (* n1 d2))
    (setq den (* d1 n2))
    (if redux (list num den)
          (semplifica (list num den))
    )
  )
)


Solution (using integer math):

Code Select Expand


(define (bank)
  (local (e balance)
    (setq e '(106246577894593683L 39085931702241241L))
    (setq balance (-f e '(1 1)))
    (for (i 1 25)
      (setq balance (-f (*f balance (list i 1)) '(1 1)))
      (println i { } balance { } (div (first balance) (last balance)))
    )
    balance
  )
)

(bank)
;-> 1 (28074714490111201L 39085931702241241L) 0.7182818284590452
;-> 2 (17063497277981161L 39085931702241241L) 0.4365636569180905
;-> 3 (12104560131702242L 39085931702241241L) 0.3096909707542714
;-> ...
;-> 16 (2433979885881703L 39085931702241241L) 0.06227253080274239
;-> 17 (2291726357747710L 39085931702241241L) 0.05863302364662064
;-> 18 (2165142737217539L 39085931702241241L) 0.05539442563917152
;-> 19 (2051780304892000L 39085931702241241L) 0.05249408714425882
;-> 20 (1949674395598759L 39085931702241241L) 0.04988174288517631
;-> 21 (1857230605332698L 39085931702241241L) 0.04751660058870241
;-> 22 (1773141615078115L 39085931702241241L) 0.04536521295145283
;-> 23 (1696325444555404L 39085931702241241L) 0.04339989788341503
;-> 24 (1625878967088455L 39085931702241241L) 0.04159754920196069
;-> 25 (1561042474970134L 39085931702241241L) 0.03993873004901714
;-> (1561042474970134L 39085931702241241L)

Now the result is correct.

[rickyboy]

That's very cool.  And a good thing to keep in mind when doing computer arithmetic.



But, your code could be much shorter -- it would then be more readable to others (and to you, months or years later :).

Code Select Expand

(define (rat n d)
  (let (g (gcd n d))
    (map (curry * 1L)
         (list (/ n g) (/ d g)))))

(define (+rat r1 r2)
  (rat (+ (* (r1 0) (r2 1))
          (* (r2 0) (r1 1)))
       (* (r1 1) (r2 1))))

(define (-rat r1 r2)
  (rat (- (* (r1 0) (r2 1))
          (* (r2 0) (r1 1)))
       (* (r1 1) (r2 1))))

(define (*rat r1 r2)
  (rat (* (r1 0) (r2 0))
       (* (r1 1) (r2 1))))

(define (/rat r1 r2)
  (rat (* (r1 0) (r2 1))
       (* (r1 1) (r2 0))))

(define (bank)
  (letn (e (rat 106246577894593683L
                39085931702241241L)
         balance (-rat e (rat 1 1)))
    (for (i 1 25)
      (setq balance (-rat (*rat balance (rat i 1))
                          (rat 1 1)))
      (println i { } balance { }
               (div (first balance)
                    (last balance))))
    balance))


[cameyo]

Thanks rickyboy.

I'll use your rational functions ;-)