Continued Fractions
Continued FractionsMoved from RationalApproximationsOfPi. See also ValueOfPi
Continued fractions provide the "best" approximation to decimals -- but I can't seem to recall the criteria for "best" that it provides; remind me, someone? Anyway, they're very easy to implement. Here's a quick little C program I hacked out to approximate pi. -- DougMerritt
Output (note that the first four entries are indeed on the previous list on RationalApproximationsOfPi of best rational approximations of pi):
1 steps: 3/ 1 = 3.00000000000000000000 (0.14159265358979311600) 2 steps: 22/ 7 = 3.14285714285714279370 (0.00126448926734967770) 3 steps: 333/ 106 = 3.14150943396226400850 (0.00008321962752910750) 4 steps: 355/ 113 = 3.14159292035398252096 (0.00000026676418940497) 5 steps: 103993/ 33102 = 3.14159265301190249176 (0.00000000057789062424) 6 steps: 104348/ 33215 = 3.14159265392142117435 (0.00000000033162805835) 7 steps: 208341/ 66317 = 3.14159265346743676872 (0.00000000012235634728) 8 steps: 312689/ 99532 = 3.14159265361893647039 (0.00000000002914335440) 9 steps: 833719/ 265381 = 3.14159265358107786525 (0.00000000000871525074) 10 steps: 1146408/ 364913 = 3.14159265359140382756 (0.00000000000161071156) 11 steps: 4272943/ 1360120 = 3.14159265358938899482 (0.00000000000040412118) 12 steps: 5419351/ 1725033 = 3.14159265358981532046 (0.00000000000002220446) 13 steps: 80143857/25510582 = 3.14159265358979267191 (0.00000000000000044409) 14 steps:245850922/78256779 = 3.14159265358979311600 (0.00000000000000000000)The actual value of pi is 3.1415926535897932384626...Basically we run out of bits in our floating point. The output cannot be better than the internal representation. Bignums could of course be used to give us any accuracy we desire.
Source:
typedef struct {
int top;
int bottom;
} Frac;
void invert(Frac *f) {
/* Inverting a fraction just means swapping numerator/denominator */
int tmp = f->bottom;
f->bottom=f->top;
f->top=tmp;
}
void add(Frac *f, int i) {
/* To add int to fraction, first multiply by denominator: */
f->top += (i * f->bottom);
}
#define PI 3.1415926535897932384626
continued_fraction(int limit) {
double approx, pi = PI;
int i, series[50];
Frac f;
/* Find the continued fraction up to "limit" terms: */
for (i=0; i<limit; i++) {
series[i] = (int) pi;
pi = 1.0 / (pi - (int) pi);
}
/* Convert the continued fraction to a plain rational: */
f.top = 1;
f.bottom = 0;
for (i=limit-1; i>=0; i--) {
invert(&f);
add(&f, series[i]);
}
/* Show results: */
approx = ((double) f.top) / ((double) f.bottom);
printf("%2d steps: %8d/%8d = %1.20f (%1.20f)\n", limit, f.top, f.bottom,
approx, fabs(PI-approx));
}
main() {
int i;
for (i=1; i<15; i++)
continued_fraction(i);
}
Some discussion of continued fractions as representations of pi:
http://www.isi.edu/~johnh/ABOUT/FEATURES/RATIONAL_PI/fast_convergence.txt
http://www.phy6.org/outreach/edu/contfrac.htm
The same program working on e instead of pi:
1 steps: 2/ 1 = 2.00000000000000000000 (0.71828182845904509080) 2 steps: 3/ 1 = 3.00000000000000000000 (0.28171817154095490920) 3 steps: 8/ 3 = 2.66666666666666651864 (0.05161516179237857216) 4 steps: 11/ 4 = 2.75000000000000000000 (0.03171817154095490920) 5 steps: 19/ 7 = 2.71428571428571441260 (0.00399611417333067820) 6 steps: 87/ 32 = 2.71875000000000000000 (0.00046817154095490920) 7 steps: 106/ 39 = 2.71794871794871806259 (0.00033311051032702821) 8 steps: 193/ 71 = 2.71830985915492950866 (0.00002803069588441787) 9 steps: 1264/ 465 = 2.71827956989247310204 (0.00000225856657198875) 10 steps: 1457/ 536 = 2.71828358208955211950 (0.00000175363050702870) 11 steps: 2721/ 1001 = 2.71828171828171827329 (0.00000011017732681751) 12 steps: 23225/ 8544 = 2.71828183520599253598 (0.00000000674694744518) 13 steps: 25946/ 9545 = 2.71828182294394959939 (0.00000000551509549140) 14 steps: 49171/ 18089 = 2.71828182873569579314 (0.00000000027665070235) 15 steps: 517656/ 190435 = 2.71828182844540133800 (0.00000000001364375279) 16 steps: 566827/ 208524 = 2.71828182847058386074 (0.00000000001153876994) 17 steps: 1084483/ 398959 = 2.71828182845856325400 (0.00000000000048183679) 18 steps: 13580623/ 4996032 = 2.71828182845906507481 (0.00000000000001998401) 19 steps: 14665106/ 5394991 = 2.71828182845902821541 (0.00000000000001687539) 20 steps: 28245729/10391023 = 2.71828182845904597897 (0.00000000000000088818) 21 steps:325368125/119696244= 2.71828182845904509080 (0.00000000000000000000)The actual value of e is 2.718281828459045235360287471352...once again we ran out of floating point bits.
Incidentally, some things, like square roots, have a lot of pattern in continued fraction representation that are not apparent in the seemingly random decimal representation. After 7 steps this program finds 239/169 as a representation for the square root of two...but if we look at the continued fraction form it found before reducing, it turns out to be 1 + 1/(2 + 1/(2 + 1/(2 + 1/(2 + 1/(2 + ...))))) or 1 2 2 2 2 2 2 2 2 2...
Square root of 3 is 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2...sqrt(5): 2 4 4 4 4 4 4 4 4 4 4 4 ...sqrt(6): 2 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 2
ContinuedFractions are our friends. :-)
---
Working in a language with such constrained numeric types can obscure this, I think. When you have language support for rationals, it is clearer. Consider this simple code to work out the cf expansion of any target:
(defun cf (target num-steps)
(loop with x = target
for n below num-steps
collect
(multiple-value-bind (quotient remainder) (round x)
(setf x (/ remainder))
quotient)))
which returns a list of the terms in the expansion, e.g:
CL-USER> (cf pi 3) (3 7 16)which is to say 3 + 1/( 7 + 1/16) = 335/113 so we can write another little function to sum up a c.f. in this form
(defun eval-cf (list)
(loop with leading-term = (car list)
with fractional-terms = (reverse (cdr list))
for x in fractional-terms
for sum = (/ x) then (/ (+ x sum))
finally (return (+ leading-term sum))))
CL-USER> (eval-cf (cf pi 3))
355/113
then make a 'table' like above very simply.
CL-USER> (loop for n from 2 upto 10 collecting (eval-cf (cf pi n))) (22/7 355/113 104348/33215 312689/99532 1146408/364913 5419351/1725033 80143857/25510582 245850922/78256779 817696623/260280919) CL-USER> (mapcar #'float *) (3.142857 3.141593 3.1415927 3.1415927 3.1415927 3.1415927 3.1415927 3.1415927 3.1415927) CL-USER> (mapcar (lambda (x) (float x 1.0d0)) **) (3.142857142857143d0 3.1415929203539825d0 3.141592653921421d0 3.1415926536189365d0 3.141592653591404d0 3.1415926535898153d0 3.1415926535897927d0 3.141592653589793d0 3.141592653589793d0)which shows how quickly we converge to within single and double precision floats....
Agreed. Very nice...except...is using "for" and "loop" truly the nicest approach? :-( -- dm
No, good point and I should have said something about it. Pedagogically giving an iterative solution seemed sensible here, given that most people here seem to like these idioms. The problem is *naturally* recursive though so what I would actually use would probably be more like the following. Again assuming common lisp (and untested)...
(defun cf (x depth)
(unless (<= depth 0)
(multiple-value-bind (q r) (round x)
(cons q (cf2 (/ r) (1- depth))))))
(defun eval-cf (cf)
(labels ((r-eval (x cf)
(if (null cf)
x
(r-eval (/ (+ x (car cf))) (cdr cf)))))
(+ (car cf) (r-eval 0 (reverse (cdr cf))))))
There's a simple algorithm for evaluating continued fractions that works forwards rather than backwards, so that you can (e.g.) use it to calculate convergents one by one, and stop when you're done. Fits nicely with lazy evaluation.
(defun eval-cf (cf)
(let ((a 0) (b 1) (c 1) (d 0))
(loop for n in cf do
(shiftf a c (+ a (* n c)))
(shiftf b d (+ b (* n d))))
(/ c d)))
On the other hand, if you want to work backwards and don't mind consuming an amount of space proportional to the length of the continued fraction, you can simplify the recursion above. (It won't be a tail recursion any more, but you get to lose the call to REVERSE, which conses O(n) space too.)
(defun eval-cf (cf)
(let ((int-part (first cf))
(tail (rest cf)))
(+ int-part (if tail (/ (eval-cf tail)) 0)) ))
What we really want to write for that is
(defun eval-cf (cf)
(if (null cf)
infinity
(+ (first cf) (/ (eval-cf (rest cf))))))
but there's no such thing as INFINITY in CL :-).
It isn't only roots of quadratic equations that have simple patterns in their continued fraction expansions, though a periodic expansion implies that the number is the root of a quadratic equation. For instance, the continued fraction for e goes [2,1,2,1,1,4,1,1,6,1,1,...], which is better written as [1,0,1, 1,2,1, 1,4,1, 1,6,1, ...].
The sizes of coefficients in the continued fraction expansion are closely related to how closely a number can be approximated by rationals. (Because if a large number appears, then truncating the c.f. just before that number gives a good approximation. The 292 that occurs early in the expansion for pi is the reason why 355/113 is so close.) So sqrt(2) = [1,2,2,2,...] is especially badly approximated by rationals. That's probably why the musical interval of an augmented fourth sounds so nasty :-). The smallest coefficients of all come from [1,1,1,1,...] = (1+sqrt(5)/2) ~= 1.618033989, the so-called "golden ratio". I can confirm that playing two musical tones whose frequencies are in that ratio does indeed sound unpleasant.
Doug asked about the criterion for "best". If p/q is a continued fraction convergent to x, then |p'-q'x| > |p-qx| for any p',q' with 0<q'<q. Note that this is stronger than saying that |p'/q'-x| > |p/q-x|. The converse holds, too: any p,q such that the stronger property holds must be a continued fraction convergent to x.
EditText of this page
(last edited April 28, 2008)
or FindPage with title or text search