(* :Name: Algebra`PolynomialContinuedFractions` *)

(* :Title: Polynomial continued fractions *)

(* :Author: Mark Sofroniou.  *)

(* :Summary:
This package provides functions for computing continued fraction
expansions of rational functions of polynomials.
*)

(* :Context: Algebra`PolynomialContinuedFractions` *)

(* :Package Version: 1.0 *)

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

(* :History:
First version January 1997 by Mark Sofroniou.
*)

(* :Keywords:
Pade approximants, continued fractions, rational functions of polynomials.
*)

(* :Discussion:
	The continued fraction expansion of a rational function of polynomials
    is the most stable and efficient general method for numerical evaluation,
    just as Horner form is for polynomials.
*)

(* :Source:
	A. Ya. Khinchin, Continued Fractions,
	University of Chicago Press, Chicago, 1964.

    N. J. Higham, Accuracy and Stability of Numerical Algorithms,
    SIAM, Philadelphia, 1996.
*)

(* :Mathematica Version: 4.0 *)

(* :Limitations:
	The speed of polynomial division and Expand.
*)

BeginPackage["Algebra`PolynomialContinuedFractions`"]

Unprotect[ PolynomialContinuedFraction, FromPolynomialContinuedFraction ]

PolynomialContinuedFraction::usage =
"PolynomialContinuedFraction[p, x] gives the continued fraction
representation for the rational function of polynomials p in the variable x.
PolynomialContinuedFraction[p, x, n] gives the continued fraction
representation through order n."

FromPolynomialContinuedFraction::usage =
"FromPolynomialContinuedFraction[{a0,{b1,a1},{b2,a2},...}] reconstructs a
nested polynomial from the list of its continued fraction terms."

Begin["`Private`"]

(* Locally defined predicates *)

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

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

(* Don't restrict these to integers; we need the capability of transforming
 symbolic forms for the quadratic irrational case and for converting
 continued fractions of polynomials. *)

CFPolyQ[{}] = True;

CFPolyQ[x_]:=
  MatchQ[x,{_?ScalarQ,{_?ScalarQ,_?ScalarQ}...}];

(* Continued fractions of rational functions of polynomials. *)

PolynomialContinuedFraction[p_, v_?ScalarQ]:=
  Module[{res},
    res = (!NumberQ[v] && RatPolyQ[p, v]);
    res /; (res && (res = cfrp[p, v])=!=$Failed)
  ];

PolynomialContinuedFraction[p_, v_?ScalarQ, n_?PMIntegerQ]:=
  Module[{res},
    res = (!NumberQ[v] && RatPolyQ[p, v]);
    res /; (res && (res = cfrp[p, v, n])=!=$Failed)
  ];

RatPolyQ[p_, v_]:=
  (PolynomialQ[Numerator[p], v] && PolynomialQ[Denominator[p], v]);

polydiv = Developer`PolynomialDivision;

Default[cfrp, 3] = Infinity;

cfrp[p_, v_, n_.] :=
  Module[{i, ll, leadquo, pd, pn, pqr},
    ll = {};
    i = 0;
    pn = normalform[ Numerator[p], v];
	pd = normalform[ Denominator[p], v];

    pqr = polydiv[ pn, pd, v ];
    If[ !ListQ[ pqr ], Return[ $Failed ] ];
	pqr = normalform[pqr, v];
    leadquo = First[pqr];
    pn = pd;
	pd = Last[ pqr ];

    While[ UnsameQ[pd, 0] && i++ < n,
      pqr = polydiv[ pn, pd, v ];
      If[ !ListQ[pqr], Return[ $Failed ] ];
      pqr = normalpqr[pqr, v];
      ll = {ll, First[pqr]};
	  pn = pd;
	  pd = Last[ pqr ];
    ];

    Join[ {leadquo}, Partition[ Flatten[ll], 2 ] ]
  ];

SetAttributes[normalform, Listable];
normalform[p_, v_]:= Collect[p, v, Together];

(* Normalize according to the leading coefficient of the quotient. *)

normalpqr[{pq_, pr_}, v_]:=
  With[{lc = Coefficient[pq, v, Exponent[pq, v]]},
    If[ UnsameQ[lc, 0],
	  {{Divide[1, lc], normalform[pq/lc, v]}, normalform[pr/lc, v]},
	  {{1, pq}, pr}
	]
  ];

(* Conversion to nested form *)

FromPolynomialContinuedFraction[{}] = 0;

FromPolynomialContinuedFraction[{x_?ScalarQ}] := x;

FromPolynomialContinuedFraction[pre_?CFPolyQ] :=
  First[pre] + Apply[Divide,Last[pre]] /; Length[pre]==2;

FromPolynomialContinuedFraction[pre_?CFPolyQ] :=
  With[{leadquo = First[pre], terms = Rest[pre]},
    leadquo +
      Fold[
        (First[#2]/(Last[#2] + #1))&,
        Apply[ Divide, Last[terms] ],
        Reverse[ Drop[terms, -1] ]
      ]
  ];

End[] (* End `Private` context. *)

SetAttributes[
  {PolynomialContinuedFraction, FromPolynomialContinuedFraction},
  ReadProtected
]

Protect[ PolynomialContinuedFraction, FromPolynomialContinuedFraction ]

EndPackage[] (* End package context. *)

(*:Examples:

PolynomialContinuedFraction[ (1 + 1/2 x + 3/4 x^2)/(2 - 5/6 x + 7/8 x^2), x ]
PolynomialContinuedFraction[ (a0 + a1 x + a2 x^2)/(b0 + b1 x + b2 x^2), x ]

FromPolynomialContinuedFraction[{x}]
FromPolynomialContinuedFraction[{x, {y, z}}]
*)
