Computing symbolic continued fractions for rational functions with respect to a variable

There are quite a few questions here about continued fractions, so this might be a duplicate, but I honestly could not find what I want.

What I want is, having two polynomials f[t], g[t], to build from them a continued fraction following the Euclidean algorithm in the standard way. That is, for f[t]==q[t] g[t]+r[t] with degree of r[t] strictly less than that of g[t], write f[t]/g[t]==q[t]+1/(g[t]/r[t]) and then do the same with g[t] in place of f[t] and r[t] in place of g[t] repeatedly until reaching zero.

Is there some efficient way to do it? I came up with my own code that seems to work but I am not sure whether it always gives correct answers and whether it can be improved:

confrac[f_,g_,t_]:=If[Exponent[f,t]<Exponent[g,t] \[Or] Exponent[f,t]<1,Factor[f/g],
With[{qr=Map[Factor,PolynomialQuotientRemainder[f,g,t]]},
First[qr]+If[Last[qr]===0,0,1/confrac[g,Last[qr],t]]]]


As requested in a comment, here is some sample output from this:

confrac[t^5 + 2 t + 1, t^2 - 1, t]
(*
t (1+t^2)+1/(1/9 (-1+3 t)-8/(9 (1+3 t)))
*)


confrac[4x^6 t+3x^6 t^2+2x^6 t^3+x^2 t^4+2t+1,x^4 t^2-2x^2 t^3-1,x]
(*
((4+3 t+2 t^2) (2 t+x^2))/t
+1/(-((t^3 (17 t+14 t^2+8 t^3+32 t^5+26 t^6+16 t^7-4 x^2-3 t x^2-2 t^2 x^2
-16 t^4 x^2-13 t^5 x^2-8 t^6 x^2))
/(4+3 t+2 t^2+16 t^4+13 t^5+8 t^6)^2)
-(t (16+24 t+25 t^2+12 t^3-21 t^4-62 t^5-46 t^6-20 t^7-32 t^8-74 t^9
-55 t^10-24 t^11))
/((4+3 t+2 t^2+16 t^4+13 t^5+8 t^6)^2
(9 t+8 t^2+4 t^3+4 x^2+3 t x^2+2 t^2 x^2+16 t^4 x^2+13 t^5 x^2+8 t^6 x^2)))
*)


• Please add a few concrete example pairs for f[t] and g[t] along with desired outputs for each case.
– Syed
Commented Aug 28, 2023 at 3:04
• @Syed You are right, added Commented Aug 28, 2023 at 7:28

To solve this, we first apply the Euclidean algorithm. This gives a list of terms of the continued fraction. To display the fractions in a continuous fraction form we add the inverse from the end using "Fold".

As an example we divide "divid" by "divis";

divid = 1 + 2 t + 3 t^2 + 4 t^3 + 5 t^4 + 6 t^5 + 7 t^6;
divis = t^3 - 2;

terms = Reap[
Do[
{q, r} = PolynomialQuotientRemainder[divid, divis, t];
Sow[q];
divid = divis; divis = r;
, {4}]
][[2, 1]]


To display "terms" in continued fraction form, we add the inverse beginning at the end:

sol= Fold[1/#1 + #2 &, Reverse@terms]


To check if this is correct, we simplify it:

sol // Simplify


The same with "While":

r = 1;
terms = Reap[
While[r =!= 0, {q, r} =
PolynomialQuotientRemainder[divid, divis, t];
Sow[q];
divid = divis;
divis = r]
][[2, 1]]

• Would it be possible to automate this somehow? Presumably the {4} in Do can be replaced by some While statement or? Commented Aug 28, 2023 at 10:49
• Simple, why do you not try it yourself? Commented Aug 28, 2023 at 11:29
• I do not think it is simple. Commented Aug 28, 2023 at 13:50
• Sorry, I copied the old code in the Addendum, it is fixed now. You see how simple it is. Commented Aug 28, 2023 at 16:24
• Very nice, thanks! I will wait a little for other possible answers but I find this one definitely worth accepting Commented Aug 28, 2023 at 20:33