Pathological floating point problems

Q&A's, tips, howto's
Locked
cameyo
Posts: 183
Joined: Sun Mar 27, 2011 3:07 pm
Location: Italy
Contact:

Pathological floating point problems

Post by cameyo »

From https://rosettacode.org/wiki/Pathologic ... t_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 all

(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 all

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

Code: Select all

(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 all

(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
Posts: 607
Joined: Fri Apr 08, 2005 7:13 pm
Location: Front Royal, Virginia

Re: Pathological floating point problems

Post by 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 all

(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))
(λx. x x) (λx. x x)

cameyo
Posts: 183
Joined: Sun Mar 27, 2011 3:07 pm
Location: Italy
Contact:

Re: Pathological floating point problems

Post by cameyo »

Thanks rickyboy.
I'll use your rational functions ;-)

Locked