(* ::Package:: *)
(* :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-2007, 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.
*)
Message[General::obspkg, "Algebra`PolynomialContinuedFractions`"]
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`"]
issueObsoleteFunMessage[fun_, context_] :=
(Message[fun::obspkgfn, fun, context];
)
(* 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]:=
(issueObsoleteFunMessage[PolynomialContinuedFraction,"Algebra`PolynomialContinuedFractions`"];
Module[{res},
res = (!NumberQ[v] && RatPolyQ[p, v]);
res /; (res && (res = cfrp[p, v])=!=$Failed)
]);
PolynomialContinuedFraction[p_, v_?ScalarQ, n_?PMIntegerQ]:=
(issueObsoleteFunMessage[PolynomialContinuedFraction,"Algebra`PolynomialContinuedFractions`"];
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 = PolynomialQuotientRemainder;
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}] :=
(issueObsoleteFunMessage[FromPolynomialContinuedFraction,"Algebra`PolynomialContinuedFractions`"];
x);
FromPolynomialContinuedFraction[pre_?CFPolyQ] :=
(issueObsoleteFunMessage[FromPolynomialContinuedFraction,"Algebra`PolynomialContinuedFractions`"];
First[pre] + Apply[Divide,Last[pre]] /; Length[pre]==2);
FromPolynomialContinuedFraction[pre_?CFPolyQ] :=
(issueObsoleteFunMessage[FromPolynomialContinuedFraction,"Algebra`PolynomialContinuedFractions`"];
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}}]
*)