About me

About me

Feeds

RSS feed

Rational approximations

16th November 2019

The following program finds a series of progressively better rational approximations to a real number, using continued fractions [1]. I wrote it a while ago, but use it so frequently that I thought I'd share it in case it's useful.

First, cf converts a rational number to a continued fraction:

(defun cf (n)
  (let ((n (rational n))
        result)
    (loop
     (multiple-value-bind (x y) (truncate n)
       (push x result)
       (cond
        ((zerop y) (return (reverse result)))
        (t (setq n (/ 1 y))))))))

For example:

> (cf pi)
(3 7 15 1 292 1 1 1 2 1 3 1 14 3 3 2 1 3 3 7 2 1 1 3 2 42 2)

Next, rat converts a continued fraction back into a rational number:

(defun rat (f)
  (destructuring-bind (car . cdr) f
    (if (null cdr) car (+ car (/ 1 (rat cdr))))))

For example:

> (float (rat '(3 7 15 1 292 1 1 1 2 1 3 1 14 3 3 2 1 3 3 7 2 1 1 3 2 42 2)))
3.1415928

Finally, approximations uses cf and rat to find a series of progressively better approximations to its real argument:

(defun approximations (n)
  (reverse
   (maplist #'(lambda (x) 
                (let ((r (rat (reverse x))))
                  (cons r (float (- n r)))))
            (reverse (cf n)))))

For each approximation it shows the error. For example:

> (pprint (approximations pi))

((3 . 0.14159265358979312D0)
 (22/7 . -0.0012644892673496777D0)
 (333/106 . 8.32196275291075D-5)
 (355/113 . -2.667641894049666D-7)
 (103993/33102 . 5.778906242426274D-10)
 (104348/33215 . -3.3162805834763276D-10)
 (208341/66317 . 1.2235634727630895D-10)
 (312689/99532 . -2.914335439641036D-11)
 (833719/265381 . 8.715250743307479D-12)
 (1146408/364913 . -1.6107115641261772D-12)
 (4272943/1360120 . 4.04121180963557D-13)
 (5419351/1725033 . -2.220446049250313D-14)
... 

There's the schoolchild approximation to pi, 22/7, and the far better one, 355/113.

I also use it to find good approximations for units conversions; for example, km to miles:

> (pprint (approximations 1.609344))

((1 . 0.609344)
 (2 . -0.390656)
 (3/2 . 0.109344006)
 (5/3 . -0.05732262)
 (8/5 . 0.009343982)
 (29/18 . -0.0017671585)
 (37/23 . 6.483793E-4)
 (103/64 . -3.0994416E-5)
 (758/471 . 2.1457672E-6)
 (861/535 . -1.7881394E-6)
 (1619/1006 . 1.1920929E-7)
 (21908/13613 . 0.0)
...

37/23 is a pretty good approximation. 

  1. ^ Continued fraction on Wikipedia.

Next: Recursive functions without a base case

Previous: Maximum number of duplicates in a list

blog comments powered by Disqus