// 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
}