In Featherstone's book "Rigid Body Dynamics Algorithms", there is a section of Chapter 9 (specifically, 9.6) dedicated to explaining how to incorporate gears into a given dynamic model.
At a high level, you need to incorporate the gear ratios into the dynamics via a similar process to how a kinematic loop would be handled, by introducing additional joint variables for the additional moving parts and then mapping them onto the standard (typical) joint variables in the algorithm above.
The example shared in the book describes a simple geared system where there are geared motors driving the angular positions of two links, which can be described as the following set of four moving components (assuming small gears drive large gears at each joint and the stator of the first motor is fixed):
- The first link, the first large gear, and the stator of the second motor;
- The first small gear and the rotor of the first motor;
- The second link and the second large gear; and
- The second small gear and the rotor of the second motor.
This means that joints 1 and 3 are the typical joints you would expect, while joints 2 and 4 are rotations within the motors driving the gears, and leads to a relationship between the joint positions in $q$ and the gear ratios $n_1$ and $n_2$. Specifically, if we define independent variables as $y_1 = q_1$ and $y_2 = q_3$, then we can define a function $\gamma$ to map from $y$ to $q$ as $\gamma = [\begin{matrix} y_1 & n_1 y_1 & y_2 & n_2 y_2 \end{matrix}]^T$.
Once we've done that, we can define the matrix $G = {\partial \gamma}/{\partial y}$ and the vector $g = \dot G \dot y$ to get to the following dynamic equation of motion (where $H$ and $C$ are the coefficients of the equation of motion without gearing):
$$G^T H G \ddot y + G^T (C+Hg) = G^T \tau$$.
Note: In this formulation, $\tau_1 = \tau_3 = 0$ while the motor torques will appear in $\tau_2$ and $\tau_4$. Using this equation, you can then apply the inverse dynamics equations as needed, such as the ABA you've outlined above.