(* :Copyright: Copyright 1988-1999, Wolfram Research, Inc. *)

(* :Mathematica Version: 4.0 *)

(* :Package Version: 2.1 *)

(* :Name: NumberTheory`ContinuedFractions` *)

(* :Author: Mark Sofroniou *)

(* :Summary:
This package provides functions for computing continued fraction convergents
of numbers and some typesetting functions for use in conjunction with
ContinuedFraction and RealDigits.
*)

(* :History:
V2.1 May 1998 by John M. Novak.
    Corrected some bugs in formatting; completely revised PeriodicForm;
    added formatting of periodic continued fractions; fixed bugs in
    Convergents.
Modified April 1998 by Mark Sofroniou.
    Moved functionality for periodic base expansions into RealDigits. 
Considerably extended January, September 1997.
    Moved ContinuedFractions into the kernel and added the
    Convergents function (suggested by Stan Wagon).
Totally rewritten March 1996 by Mark Sofroniou.
	Added PeriodicForm and ToPeriodicForm.
	Rewrote ContinuedFraction to unravel recursion and added
	unrestricted continued fractions for rationals.
Modified May 1992 by John M. Novak:
	Allow nonrational nonfloating-point numbers (estimates precision)
	Force machine numbers to bignums (machine precision) to reduce
	or eliminate roundoff error (with suggestions from J. Keiper).
Modified March 1992 by John M. Novak:
	Allow rational numbers. Fix behavior when partial fraction
	is finite and less than the order requested.
*)

(* :Keywords:
number theory, continued fractions, convergents, rational numbers,
decimal expansions.
*)

(* :Requirements: None. *)

(* :Discussion:
    Continued fraction convergents are often used to approximate
    irrational numbers by rational ones. They have many interesting
    properties such as exponential rate of convergence in the number
    of terms and alternating approximation from above and below.
    Furthermore a convergent p/q  of a simple continued fraction
    is better than any other rational approximation with denominator
    less than q.

	The base expansion of any rational number can be written in terms
    of pre-periodic and periodic parts, either of which may be zero in
    length. For example, the decimal expansion of 1/8 terminates, but
    the expansion of 1/7 does not. 1/6 is a simple example of a rational
    number with non-zero pre-periodic and periodic parts
    (1 and 6 respectively).
*)

(* :Source:
    R. P. Brent, A. J. van der Porten and H. J. J. te Riele,
    A Comparative Study of Algorithms for Computing Continued
    Fractions of Algebraic Numbers. in Proceedings of the Second
    International Symposium on Algorithmic Number Theory (ANTS II),
    edited by H. Cohen, Springer Verlag, 1996, pp. 35--47.
	A. Ya. Khinchin, Continued Fractions,
	University of Chicago Press, Chicago, 1964.
	H. E. Rose, A Course in Number Theory,
	Second Edition, Oxford University Press, Oxford, 1994.
	K. H. Rosen, Elementary number theory and its applications,
	Third edition, Addison Wesley, Reading Mass., 1993.
	P. Shiu, Computation of continued fractions without input values,
	Math. Comp. 64 211, 1307--1317, 1995.
	S. Wagon, Mathematica in Action, W. H. Freeman, New York, 1991.
*)

(* :Limitations:
*)

BeginPackage["NumberTheory`ContinuedFractions`"]

Unprotect[ ContinuedFractionForm, Convergents, PeriodicForm,
QuadraticIrrationalQ ]

ContinuedFractionForm::usage =
"ContinuedFractionForm[ cf ] represents a continued fraction as given
by ContinuedFraction.  The representation may be converted to a
rational number using Normal."

Convergents::usage =
"Convergents[x] gives a list of the continued fraction
convergents for the number or continued fraction expansion x.
Convergents[x, n] gives a list of the first n continued
fraction convergents for the number or continued fraction expansion x.
Note that the order returned may be less than n if the continued fraction
terminates before n steps."

PeriodicForm::usage = 
"PeriodicForm[{{d1, d2, ...,dn, {p1, p2, ..., pn}}, exp}] represents the
decimal expansion of a rational in terms of the preperiodic and periodic
parts as given by RealDigits. PeriodicForm[digits, b] represents the
expansion in base b. The representation may be converted to a rational
number using Normal."

QuadraticIrrationalQ::usage =
"QuadraticIrrationalQ[x] gives True if x is a quadratic irrational
and False otherwise."

Begin["`Private`"]

Convergents::terms =
"`` does not have a terminating or periodic continued fraction expansion;
specify an explicit number of terms to generate."

PMIntegerQ[x_]:= Developer`MachineIntegerQ[x] && Positive[x];

(* Predicate for testing quadratic irrationals. *)

qi[_Integer] = True;
qi[_Rational] = True;
qi[Power[x_, (1/2 | -1/2)]]:= qi[x];
qi[x:(_Plus|_Times)]:= Apply[ And, Map[qi, Apply[List, x] ] ];
qi[_] = False;

QuadraticIrrationalQ[n_]:=
  Block[{$RecursionLimit = Infinity},
    If[ !NumericQ[n] || NumberQ[n], False, qi[RootReduce[n]] ]
  ];

(**** Convergents of continued fractions ****)

ScalarQ[_List]:= False;
ScalarQ[_]:= True;

CFQ[_?VectorQ]:= True;
CFQ[{___?ScalarQ, _?VectorQ}]:= True;
CFQ[_] = False;

ValidCFNumberQ[_Integer]:= True;
ValidCFNumberQ[_Rational]:= True;
ValidCFNumberQ[_Real]:= True;
ValidCFNumberQ[_]:= False;

Convergents[x_?ValidCFNumberQ]:=
  Module[{res},
    res = cnvgts[ ContinuedFraction[ x ] ];
    res /; (res =!= $Failed)
  ];

Convergents[ x_?QuadraticIrrationalQ ]:=
  Module[{xr, res},
    xr = RootReduce[ x ];
    If[ NumberQ[xr], (* RootReduce found a reduction to a number *)
      If[ ValidCFNumberQ[xr],
        res = Convergents[ xr ],
        res = $Failed
      ],
      If[ NumericQ[xr],
        res = cnvgtsqiinf[ xr, ContinuedFraction[ xr ] ],
        res = $Failed (* RootReduce failed. *)
      ];
    ];
    res /; (res =!= $Failed)
  ];

Convergents[x_?(NumericQ[#]&&!NumberQ[#]&)]:=
  x /; (Message[Convergents::terms,x]; False);

Convergents[{}] = {0};

Convergents[{{}}] = {0};

Convergents[cf_?CFQ]:=
  Module[{res, pre}, 
    If[ VectorQ[ cf ],
      res = cnvgts[ cf ],
      If[ Last[ cf ]==={},
        res = cnvgts[ Drop[ cf, -1 ] ],
        res = FromContinuedFraction[ cf ]; (* Attempt quadratic irrational case *)
        If[ NumericQ[ res ],
          res = cnvgtsqiinf[ res, cf ],
          res = $Failed
        ]
      ]
    ];
    res /; res =!= $Failed
  ];

(* Specified number of terms *)

Convergents[ x_?ValidCFNumberQ, n_?PMIntegerQ]:=
  Module[{res},
    res = cnvgts[ ContinuedFraction[ x, n ] ];
    res /; (res =!= $Failed)
  ];

Convergents[ x_?QuadraticIrrationalQ, n_?PMIntegerQ]:=
  Module[{xr, res},
    xr = RootReduce[ x ];
    If[ NumberQ[xr], (* RootReduce found a reduction to a number *)
      If[ ValidCFNumberQ[xr],
        res = Convergents[ xr, n ],
        res = $Failed
      ],
      If[ NumericQ[xr],
        res = cnvgtsqi[ ContinuedFraction[ xr ], n ],
        res = $Failed (* RootReduce failed. *)
      ];
    ];
    res /; (res =!= $Failed)
  ];

(* This definition must come after the attempted reduction of
 numeric expressions by RootReduce. *)

Convergents[ x_?(NumericQ[#]&&!NumberQ[#]&), n_?PMIntegerQ]:=
  Module[{res},
    res = cnvgts[ ContinuedFraction[ x, n ] ];
    res /; (res =!= $Failed)
  ];

(* A specified continued fraction. *)

Convergents[{}, _?PMIntegerQ]:= {0};

Convergents[{{}}, _?PMIntegerQ]:= {0};

Convergents[cf_?CFQ, n_?PMIntegerQ]:=
  If[ VectorQ[ cf ],
    cnvgts[ cf, n ],
    If[ Last[ cf ] === {},
      cnvgts[ Drop[ cf, -1 ], n ],
      cnvgtsqi[ cf, n ] (* Quadratic irrational case *)
    ]
  ];

(* allow convergents of a ContinuedFractionForm object *)
Convergents[ContinuedFractionForm[arg_], rest___] :=
    Convergents[arg, rest]

(* Continued fraction convergents of numbers. *)

(* Compute the convergents from the list of partial quotients. *)

cterms[cfl_] :=
  Module[{mat, next},
    mat = {{First[cfl], 1}, {1, 0}};
    FoldList[
      (next = {#2, 1}.mat; mat = {next, #1}; next)&,
      First[mat],
      Rest[cfl]
    ]
  ];


(* Use the list of partial quotients if it was specified. *)

divide[n_Integer, m_Integer] = System`Private`RationalNoReduce[n, m];
divide[n_, m_]:= Divide[n, m];

Default[cnvgts,2] = Infinity;

cnvgts[{}, ___] := {}

cnvgts[x_?ListQ, n_.]:=
  Module[{cf},
    If[n===Infinity,
      cf = x,
      cf = Take[ x, Min[ Length[x], n ] ]
    ];
    Apply[ divide, cterms[cf], 1 ]
  ];

cnvgts[_, _.]:= $Failed; (* ContinuedFraction failed *)

(* Convergents of quadratic irrationals. *)

cnvgtsqiinf[ qi_, pcf_ ]:=
  Module[{tmpcf},
    tmpcf = Join[ Drop[pcf,-1], Drop[Last[pcf],-1] ];
    Join[ cnvgts[ tmpcf ], {qi} ]
  ];

(* Finite number of terms. *)

cnvgtsqi[ pcf_, n_ ]:=
  Module[{len, len1, len2, quo, rem, tpcf1, tpcf2},
    len = n;
    tpcf1 = Drop[pcf,-1]; (* Preperiodic part *)
    len1 = Length[tpcf1];
    If[ len <= len1,
      Return[ cnvgts[Take[tpcf1,len] ] ],
      len -= len1
    ];
    tpcf2 = Last[pcf]; (* Periodic part *)
    len2 = Length[tpcf2];
    quo = Quotient[len,len2];
    rem = len - len2*quo;
    If[ quo > 0, tpcf1 = Join[ tpcf1, Flatten[ Table[tpcf2,{quo}] ] ] ];
    If[ rem > 0, tpcf1 = Join[ tpcf1, Take[ tpcf2, rem ] ] ];
    cnvgts[ tpcf1 ]
  ];

(***** ContinuedFractionForm formatting wrapper *****)

(* Function to reconstruct a rational from a continued fraction expansion. *)
ContinuedFractionForm /:
  Normal[ ContinuedFractionForm[args__] ]:=
    Module[{res},
      res = FromContinuedFraction[args];
      res /; Head[res] =!= FromContinuedFraction
    ];

ContinuedFractionForm /:
   FromContinuedFraction[ContinuedFractionForm[args__]] :=
       FromContinuedFraction[args]

(* One problem to avoid is that the normal typesetting rules
 reduce the font size of deeply nested fractions. *)

ContinuedFractionForm /:
  MakeBoxes[ cf:ContinuedFractionForm[a_?VectorQ], fmt_]:=
    InterpretationBox[#, cf]&[
      StyleBox[
        buildCFboxes[a],
        ScriptSizeMultipliers -> 1
      ]
    ]

(* formatting for cf with periodic part; notation due to Lou D'Andria *)

ContinuedFractionForm /:
  MakeBoxes[ cf:ContinuedFractionForm[{a___, per_?VectorQ}], fmt_]:=
    InterpretationBox[#, cf]&[
      StyleBox[
        If[Length[{a}] === 0,
            RowBox[{"(", buildCFboxes[Append[per, "\[Ellipsis]"]], ")"}],
            buildCFboxes[{a,
                RowBox[{"(", buildCFboxes[Append[per, "\[Ellipsis]"]], ")"}]}
            ]
        ],
        ScriptSizeMultipliers -> 1, SpanSymmetric -> False
      ]
    ]

buildCFboxes[{a___, b_}] :=
    Fold[RowBox[{#2, "+", FractionBox[1, #1]}]&,
          b,
          Reverse[{a}]
    ]/. (n_Integer :> ToString[n])

(***** PeriodicForm formatting wrapper *****)

(* Function to reconstruct a rational from a periodic base expansion. *)
PeriodicForm /:
  Normal[ PeriodicForm[
        digs:{{___Integer, {___Integer}} | {___Integer}, _Integer},
        rest___] ]:=
    Module[{res},
      res = FromDigits[digs, rest];
      res /; Head[res] =!= FromDigits
    ];

(* allow alternate form without exponent *)

PeriodicForm/:
  Normal[ PeriodicForm[digs:({___Integer, {___Integer}} | {___Integer}),
          rest___]] :=
      Normal[PeriodicForm[{digs, 0}, rest]]

(* Formatted version of periodic numbers (based on RealDigits output) *)
(* allow form without exponent specified; note that the Interpretation box
   needs to interpret as the original pf, not the one with the 0 exponent *)

PeriodicForm/:
  MakeBoxes[pf:PeriodicForm[digs:{___Integer, {___Integer}}, base_Integer:10],
            rest___] :=
       MakeBoxes[PeriodicForm[{digs, 0}, base], rest]/.
           HoldPattern[InterpretationBox][f_, p_, r___] :>
               InterpretationBox[f, pf, r]

(* allow form without periodic part *)
PeriodicForm/:
  MakeBoxes[pf:PeriodicForm[{{digs___Integer},exp_}, base_Integer:10],
            rest___] :=
      MakeBoxes[PeriodicForm[{{digs, {}}, exp}, base]]/.
           HoldPattern[InterpretationBox][f_, p_, r___] :>
               InterpretationBox[f, pf, r]

(* general form, with periodic part and exponent -
   note that 'ZeroWidthTimes' spacing is used in the row boxes if the
   base is <= 36 (which means a single character per digit); note that
   the product with the exponent is not added unless the exponent is
   between -6 and 6 (unless the nonperiodic part doesn't support moving
   the decimal appropriately). *)

PeriodicForm/:
  MakeBoxes[pf:PeriodicForm[
          digs:{{pre___Integer, per:{___Integer}}, exp_Integer},
          base_Integer:10],
      rest___] :=
      InterpretationBox[#, pf]&[
         If[base <= 36, StyleBox[#, ZeroWidthTimes -> True], #]&[
             If[(-6 <= exp <=Length[{pre}]) && (exp <= 6),
                 RowBox[pform[{pre}, per, exp, base]],
                 RowBox[{RowBox[pform[{pre}, per, 0, base]], "\[Times]",
                         SuperscriptBox[ToString[base], ToString[exp]]}]
             ]
          ]
      ]

(* formatting subutility, assembles list of digits with decimal point
    and base subscript (via fixbases) *)
(* empty number *)
pform[{}, {}, _, _] := {"0"}

(* following rule is redundant, but should be most called so slight
   performance improvement by including it, I think *)
pform[pre_, per_, 0, base_] :=
    fixbases[Flatten[{0, ".", pre, Overscribe[per]}], base]

(* not called in the following form unless n >= Length[pre]; shift
   decimal point right *)
pform[pre_, per_, n_?Positive, base_] :=
    fixbases[Flatten[{Take[pre, n], ".", Drop[pre, n], Overscribe[per]}],
             base]

(* shift decimal point left *)
pform[pre_, per_, n_, base_] :=
    fixbases[Flatten[{0, ".", Table[0, {Abs[n]}], pre, Overscribe[per]}],
             base]

(* fixbases also handles stringification of numbers, as well as
   formatting digits to the base and adding a subscript noting the base *)
fixbases[digs_, 10] := digs/. (i_Integer :> ToString[i])

(* for bases between 1 and 36 (other than 10), use letters for digits *)
$digitreplacements = Dispatch[
{0 -> "0", 1 -> "1", 2 -> "2", 3 -> "3", 4 -> "4", 5 -> "5", 6 -> "6",
 7 -> "7", 8 -> "8", 9 -> "9", 10 -> "A", 11 -> "B", 12 -> "C", 13 -> "D",
 14 -> "E", 15 -> "F", 16 -> "G", 17 -> "H", 18 -> "I", 19 -> "J", 20 -> "K",
 21 -> "L", 22 -> "M", 23 -> "N", 24 -> "O", 25 -> "P", 26 -> "Q", 27 -> "R",
 28 -> "S", 29 -> "T", 30 -> "U", 31 -> "V", 32 -> "W", 33 -> "X", 34 -> "Y",
 35 -> "Z"}
];

fixbases[digs_, b_?(# <= 36 &)] :=
    Append[Drop[#,-1], SubscriptBox[Last[#], ToString[b]]]&[
        digs/.$digitreplacements
    ]

fixbases[digs_, b_] :=
    Append[Drop[#,-1], SubscriptBox[Last[#], ToString[b]]]&[
        digs/.(d_Integer :> UnderscriptBox[ToString[d], "\[UnderBracket]"])
    ]

Overscribe[{}] = {};

Overscribe[per_]:=
  Map[OverscriptBox[#, "_"]&, per]

End[]  (* NumberTheory`ContinuedFractions`Private` *)

SetAttributes[{ContinuedFractionForm, PeriodicForm}, NHoldAll];

SetAttributes[QuadraticIrrationalQ, Listable];

SetAttributes[ {ContinuedFractionForm, Convergents, PeriodicForm,
QuadraticIrrationalQ}, ReadProtected ]

Protect[ ContinuedFractionForm, Convergents, PeriodicForm,
QuadraticIrrationalQ ]

EndPackage[]  (* NumberTheory`ContinuedFractions` *)

(*:Examples:

QuadraticIrrationalQ[ 1/4 (5 + Sqrt[31]) ]

Convergents[ 101/9801 ]
Convergents[ 10201/96059601, 10 ]
Convergents[ N[ Pi ] ]
Convergents[ {2, {5, 1, 1, 5, 3}} ]
Convergents[ 2 + Sqrt[3] - Sqrt[2] - Sqrt[5 - 2*Sqrt[6]], 20 ]
Convergents[ Pi, 10 ]

Normal[ ContinuedFractionForm[{1, 1, 2, 3}] ]
Normal[ ContinuedFractionForm[{1, 1, {2, 3}}] ]

PeriodicForm[ {{4, 2, {0, 3, 1, 7, 4}}, 0} ]
PeriodicForm[ {{4, 2, {0, 3, 1, 7, 4}}, 0}, 8 ]

Normal[ PeriodicForm[ {1, 2, {3, 4}} ] ]
Normal[ PeriodicForm[ {1, 2, {3, 4}}, 5 ] ]
*)
