  We analyze Implicit Repeated Squaring ($\textbf{IRS}$) -- an algorithm that
implicitly computes $(A^{-1}B)^{2^p}$ for $A,B \in {\mathbb C}^{n \times n}$
and $p \in {\mathbb N}_+$. Originating with inverse-free divide-and-conquer
eigensolvers, this routine uses only matrix multiplication and QR to produce
$A_p$ and $B_p$ satisfying $A_p^{-1}B_p = (A^{-1}B)^{2^p}$. By offering the
choice to work with $A_p$ and $B_p$ individually or form $A_p^{-1}B_p$,
$\textbf{IRS}$ is applicable to both implicit and explicit computations
involving $(A^{-1}B)^{2^p}$. Moreover, it exhibits stability. In this paper, we
present the first rigorous and general analysis of the algorithm, showing that
a finite-arithmetic implementation of $\textbf{IRS}$ -- using any fast matrix
multiplication routine -- satisfies a mixed forward/backward stability bound
provided $O(\log(1/\epsilon) + \log(\mu(n)) + (p-1)
\log(\kappa_\text{IRS}(A,B,p)))$ bits of precision are used, where $\epsilon$
is the desired (sufficiently small) forward error, $\mu(n)$ is a polynomial in
$n$, and $\kappa_{\text{IRS}}$ is an appropriately chosen condition number.
Hence, fast $\textbf{IRS}$ is stable in relatively modest precision when $p$ is
logarithmic in $n$ or smaller. We also demonstrate on a handful of examples
that $\textbf{IRS}$ offers improved accuracy over a naive algorithm that
explicitly forms and squares $A^{-1}B$. Along the way, we develop a spectral
norm perturbation bound for full QR that may be of independent interest.
