continuedFraction.frink

View or download 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
}



View or download 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 17592 days, 15 hours, 43 minutes ago.