## 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.

- ^ Continued fraction on Wikipedia.

blog comments powered by Disqus