  The inspiral of two compact objects in gravitational wave astronomy is
described by a post-Newtonian expansion in powers of $(v/c)$. In most cases, it
is believed that the post-Newtonian expansion is asymptotically divergent. A
standard technique for accelerating the convergence of a power series is to
re-sum the series by means of a rational polynomial called a Pad\'e
approximation. If we liken this approximation to a matrix, the best convergence
is achieved by staying close to a diagonal Pad\'e approximation. This broadly
presents two subsets of the approximation : a super-diagonal approximation
$P^M_N$ and a sub-diagonal approximation $P_M^N$, where $M = N+\epsilon$, and
$\epsilon$ takes the values of 0 or 1. Left as rational polynomials, the
coefficients in both the numerator and denominator need to be re-calculated as
the order of the initial power series approximation is increased. However, the
sub-diagonal Pad\'e approximant is computationally advantageous as it can be
expressed in terms of a Gauss-like continued fraction. Once in this form, each
coefficient in the continued fraction is uniquely determined at each order.
This means that as we increase the order of approximation of the original power
series, we now have only one new additional coefficient to calculate in the
continued fraction. While it is possible to provide explicit expressions for
the continued fraction coefficients, they rapidly become unwieldy at high
orders of approximation. It is also possible to numerically calculate the
coefficients by means of ratios of Hankel determinants. However, these
determinants can be ill-conditioned and lead to numerical instabilities. In
this article, we present a method for calculating the continued fraction
coefficients at arbitrary orders of approximation.
