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 19581 days, 17 hours, 1 minutes ago.