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

Question

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)))
*)

Answer 1

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

Addendum

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]]