(* :Title: Continued Fractions *)

(* :Authors: Stephen Wolfram, Ilan Vardi *)

(* :Summary:
  Generates continued fraction approximations to real numbers,
  represented as ContinuedFractionForm objects.  Extends
  the Normal function to convert such objects to rational numbers.
*)

(* :Context: NumberTheory`ContinuedFractions` *)

(* :Package Version: 1.3 *)

(* :Copyright: Copyright 1990  Addison-Wesley Publishing Company. *)

(* :History:
  Extended to rational numbers by Silvio Levy, February 1990.
  Revised by Ilan Vardi (Wolfram Research), September 1989.
  Version 1.1 by Stephen Wolfram (Wolfram Research), May 1987.
*)

(* :Keywords: 
  number theory, continued fraction, rational approximation, 
  partial quotient 
*)

(* :Source:
  Hardy and Wright: "An Introduction to the Theory of Numbers",
  Cambridge University Press.

  (This field may be repeated or absent.)
*)

(* :Warning: Redefines the function Normal. 

  (This field specifies any aspects of this package that might interfere
  with ordinary Mathematica operations. It may be repeated or absent.)
*)

(* :Mathematica Version: 1.2 *)

(* One or more :Requirement: fields may be present here,
stating packages on which this one depends (unless they are
part of the standard Mathematica distribution),
and machine or operating system dependencies. *)

(* :Limitation: Affected by floating-point roundoff in input numbers. *)
(* :Limitation:
   Does not define mathematical operations on 
   ContinuedFractionForm objects.
*)

(* :Discussion: 
The continued fraction of a real number  x  is 
ContinuedFractionForm[{a0,a1,a2,...}] if  x  equals
a0+1/(a1+1/(a2+...)).

(This field is optional.  It may elaborate on the summary, contain
detailed usage instructions, or explanations about the algorithms.)
*)

BeginPackage["NumberTheory`ContinuedFractions`"]

(* Include usage messages for each externally-visible function,
as well as error or warning messages as appropriate. *)

ContinuedFraction::usage = "\n
  ContinuedFraction[x,n] generates the continued fraction\n
  representation for the number x to order n."

ContinuedFractionForm::usage = "\n
  ContinuedFractionForm[{a0,a1,...}] is a continued fraction;\n
  it is converted to an ordinary rational using Normal."

(* Use a private context for the actual definitions.*)

Begin["`Private`"]

(* :Note:
   In ContinuedFraction, we use the symbol Infinity to flag for a
   terminating continued fraction (which corresponds to a rational
   number); in this case the length of the resulting list may be less
   than the second argument to ContinuedFraction.

(Use this field to explain any non-obvious implementation details.)
*)

next[Infinity] := Infinity
next[x_] := If[x==Floor[x], Infinity, 1/(x-Floor[x])]

ContinuedFraction[x_?NumberQ,n_Integer?Positive] :=
  ContinuedFractionForm[Floor[Select[NestList[next,x,n-1],NumberQ]]] /;
    Im[x]==0

ContinuedFractionForm/: Normal[ContinuedFractionForm[a_List]] :=
  Last[Accumulate[ 1/#1+#2 &, Reverse[a]]]

End[]
EndPackage[]

***********************************

(* Test files are numbered sequentially. *)
(* :Tests: 1 *)

(* Repeat essential information here. *)
(* :Title: Continued Fractions *)
(* :Authors: Stephen Wolfram, Ilan Vardi *)
(* :Context: NumberTheory`ContinuedFractions` *)
(* :Package Version: 1.2 *)

(* A Mathematica/ session containing the tests.  You should start by
using ResetMedium and loading your package. *)

In[1]:= ResetMedium[FormatType->InputForm, PageWidth->78]

In[2]:= Needs["NumberTheory`ContinuedFractions`"]

In[3]:= ContinuedFraction[N[1234/56789],7]

Out[3]= ContinuedFractionForm[{0, 46, 49, 2, 1, 3, 2}]

In[4]:= Normal[%]

Out[4]= 1234/56789

In[5]:= ContinuedFraction[N[Pi,100],30]

Out[5]= ContinuedFractionForm[{3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 
    1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13, 1, 4}]

In[6]:= Normal[%]

Out[6]= 30246273033735921/9627687726852338

In[7]:= N[%-Pi,80]

Out[7]= 4.5616784501006042255845466532971749289855821495*10^-33

