Download or view continuedFraction.frink in plain text format
// Generate continued fraction representation for a number, or turn a
// continued fraction back into a number. This also finds the closest fraction
// to a number in the most efficient way.
// See the book _Continued Fractions_ by A. Ya. Khinchin
// Most of the notes and references to theorems and algorithms are taken
// from this book.
//
// See http://mathworld.wolfram.com/ContinuedFraction.html
// Generate a continued fraction (as an array) for the number x.
// No more than "limit" terms will be produced.
//
// As part of this process, we build up the represented value (the
// "convergent") as we go to determine if we've exceeded the available
// working precision.
continuedFraction[x, limit=10] :=
{
r = x
result = new array
count = 0
a = floor[r]
var p
var q
p2 = 1
q2 = 0
p1 = a
q1 = 1
while (r != 0) and count < limit
{
result.push[a]
denom = r - a
if denom == 0
break
r = 1 / denom
if count > 0
{
p = a p1 + p2
q = a q1 + q2
// println["$p / $q"]
if q != 0 and p/q - x == 0 // Error is zero; we're done
break
p2 = p1
q2 = q1
p1 = p
q1 = q
}
a = floor[r]
count = count + 1
}
return result
}
// Turns an array indicating a continued fraction into a single
// fraction (or array of fractions) indicating its value.
// This performs the calculation top-down.
//
// Since the result for a canonical continued fraction will always be an
// exact rational number, it is safe to do the calculation top-down or
// bottom-up.
//
// if asArray is true, this will return the results as an array of successive
// approximations to the value (also called convergents of the continued
// fraction.)
continuedFractionToFraction[a, asArray=false] :=
{
// n is the order of the continued fraction (the index of the last element,
// starting with 0)
n = length[a] - 1
if n == 0
return a@0
if asArray == true
results = new array
// The final value will be p/q. For canonical continued fractions, p and q
// will both always be integers.
//
// An interesting note is that p/q is *already* always in simplest form,
// that is, the fraction will never need reducing. Or, in other terms,
// gcd[p,q] = 1, or in even other terms, p and q share no common factors.
var p
var q
// p1 and q1 are p and q from the k-1 th (i.e. the previous) step.
p1 = a@0
q1 = 1
// p2 and q2 are p(k-2) and q(k-2) (i.e. two steps ago.)
// These values are necessary to make the convergent work for k=1
// See Khinchin p.5 (remark)
p2 = 1
q2 = 0
if asArray == true
results.push[p1]
for k = 1 to n
{
// See Khinchin Theorem 1
p = a@k p1 + p2
q = a@k q1 + q2
// Move current values to previous values
p2 = p1
q2 = q1
p1 = p
q1 = q
if asArray == true
results.push[p/q]
}
if asArray
return results
else
return p/q
}
/** This is a "bottom-up" calculation of a continued fraction, turning it into
a fraction. It's more obvious than the above "top-down" calculation, but
this approach doesn't work with infinite continued fractions, and doesn't
show the convergence of continued fractions. (All canonical integer-valued
continued fractions will converge, see Khinchin, theorem 10)
Since for canonical continued fractions the result is an exact rational
number, calculations can be done "top-down" or "bottom-up" without loss of
accuracy. This version is preserved for obviousness and clarity.
**/
/*
continuedFractionToFraction[a] :=
{
n = length[a] - 1
frac = 0
for i = n to 1 step -1
frac = 1/(a@i+frac)
return a@0 + frac
}
*/
// Returns the closest fraction to x as a fraction by creating a continued
// fraction with no more than "limit" terms.
closestFraction[x, limit=10] :=
{
cf = continuedFraction[x,limit]
return continuedFractionToFraction[cf]
}
// Turns an array representing a continued fraction back into an array
// of fractions. Each fraction in the array shows the results of using more
// and more terms from the continued fraction representation.
continuedFractionToArray[cf] :=
{
continuedFractionToFraction[cf, true]
}
// Return an array which is a series of fractions representing
// successively-better approximations to the number x. No more than "limit"
// terms will be produced.
approximations[x, limit=10] := continuedFractionToArray[continuedFraction[x,limit]]
// Return an array which is a series of fractions representing
// successively-better approximations to the number x. No more than "limit"
// terms will be produced.
approximationsWithErrors[x, limit=10] :=
{
result = new array
approximations = continuedFractionToArray[continuedFraction[x,limit]]
for a = approximations
result.push[ [a, (a-x)/x] ]
return result
}
Download or view continuedFraction.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19973 days, 20 hours, 10 minutes ago.