The method of Viskovatov
While smooth real functions can have at most one convergent power series expansion in the neighborhood of a point (Taylor series), expanding the same function as a continued fraction presents many options. A basic technique (Viskovatov, 1806) operates on the ratio of a pair of power series (or polynomials) to produce a Thiele-type continued fraction or with slight variation, a corresponding-type continued fraction:
$$ p(x)/q(x) = c_0 + \cfrac{x^{\alpha_1}}{c_1 + \cfrac{x^{\alpha_2}}{c_2 +
\cfrac{x^{\alpha_3}}{c_3 + \cdots}}} $$
where it is assumed $q(0)$ is nonzero.
Specifically, after subtracting off the constant term $c_0 = p(0)/q(0)$:
$$ p(x)/q(x) = c_0 + \frac{p(x) - c_0 q(x)}{q(x)} $$
If the latter numerator is not identically zero, then $p(x) - c_0 q(x) = x^{\alpha_1} r(x)$ since it has a root at $x=0$ with multiplicity $\alpha_1 \ge 1$. Thus in that case we expand:
$$ p(x)/q(x) = c_0 + \frac{x^{\alpha_1}}{q(x)/r(x)} $$
By construction $r(0)$ is nonzero, so the method applies to $q(x)/r(x)$ and may be continued until potentially an identically zero numerator is reached at some stage. Terminating fractions are then essentially rational expressions in $x$.
Truncating such an expansion, omitting the ellipsis at some stage, gives rational expressions analogous to nested multiplication forms of polynomials in that divisions replace the multiplies, especially in the simple case where exponents $\alpha_i$ are always one:
c0 + x/(c1 + x/(c2 + x/(c3 + ...)))
Some regularity can be seen in the Thiele fraction expansion around $x=0$ of the function $f(x)/x^3 = \sum_{k=0}^{\infty} \frac{B_k x^k}{k!(k+3)}$:
$$ \frac{f(x)}{x^3} = \frac{1}{3} + \cfrac{x}{-8 + \cfrac{x}{-\frac{15}{16} +
\cfrac{x}{8 + \cfrac{x}{\frac{7}{5} + \cfrac{x}{-8 + \cfrac{x}{-\frac{27}{16} +
\cfrac{x}{8 + \cfrac{x}{\frac{275}{161} + \cdots}}}}}}}} $$
The curious reader will find a few more of the constants:
$$ c_9 = -8, c_{10} = -\frac{144417}{97264} $$
$$ c_{11} = 8, c_{12} = \frac{36954241}{25202043} $$
$$ c_{13} = -8, c_{14} = -\frac{1763915461119}{1034996631248} $$
$$ c_{15} = 8, c_{16} = \frac{53785750138088275}{33572832373135209} $$
$$ c_{17} = -8, c_{18} = -\frac{559934841935054771682651}{394380619248992328976208} $$
$$ c_{19} = 8,
c_{20} = \frac{123408586919942234000562495331463}{74004329695604288642520453854805} $$
seem to dispose of easy conjectures about the even index constants $c_{2k}$ while also confirming the alternating appearances of $\pm 8$ in odd indices. For the profoundly curious reader I have an Appendix below that describes my implementation in Maxima CAS of the method of Viskovatov.
For the sake of comparison with the Maclaurin approximations given in the Background section, I've chosen Thiele-type expansions of orders 4, 12, and 20 (as multiplied by removed factor $x^3$) to match up with polynomial degrees 7, 15, and 23 previously shown. In one sense this is a conservative choice, as the higher odd index Bernoulli coefficients being zero would allow those same Maclaurin polynomials to be construed as approximations of degree 8, 16, and 24.
Fig. 2 Debye function of 3rd order and its Thiele-type fractions
The plot demonstrates that increasing the order of expansion does not "hit a wall" in the way the power series did, that the domain of approximation appears to extend smoothly with increasing orders.
Appendix: Maxima CAS implementation
I always anticipate barking my shins at every turn when trying to perform symbolic calculations with a CAS, and this experience with Maxima was in line with that. Maxima has no datatype for continued fractions, but it does provide Taylor series constructions and a function to compute Pade approximations. The latter was useful as a check on the results of my Viskovatov method.
This is a user-defined Maxima function that returns a rational expression:
/*
Maxima function that performs Thiele fraction expansion by
applying Viskovatov's method to the ratio of two power series.
viskovatovRat(Variable,PowerSeries0,PowerSeries1,Depth)
returning a rational expression
Expands about Variable=0; assumes Depth is nonnegative integer.
Assumes value of PowerSeries1 at Variable=0 is nonzero.
*/
viskovatovRat(x,p0,p1,n) :=
block ( [a,c,p2],
c : ev(p0,[x=0])/ev(p1,[x=0]),
if n <= 0
then return(c),
p2 : ratsimp(p0 - c*p1),
a : 0,
if p2 = 0
then return(c),
while ev(p2,[x=0]) = 0
do (
a : a+1,
p2 : ratsimp(p2/x)
),
c + (x^a/viskovatovRat(x,p1,p2,n-1))
)$
I also created a variant that returns a list of constants plus powers of $x$ which was convenient for compactness of display.
Best Answer
Gosper found efficient ways to do arithmetic with continued fractions (without converting them to ordinary fractions or decimals). Here is a page with links to Gosper's work, but also with an exposition of Gosper's methods.
See also this older m.se question, Faster arithmetic with finite continued fractions