@@ -661,22 +661,13 @@ from its default.
 .. _option-flag-invdiscrete:
 
 :at:`invdiscrete`: :at-val:`[disable, enable], "disable"`
-   This dual-purpose flag enables discrete-time inverse dynamics and disables :ref:`midpoint integration<geMidpoint>`.
-
-   Enable discrete-time inverse dynamics
-     This flag **enables** discrete-time inverse dynamics with :ref:`mj_inverse` for all
-     :ref:`integrators<option-integrator>` other than ``RK4``. Recall from the :ref:`numerical
-     integration<geIntegration>` section that the one-step integrators (``Euler``, ``implicit`` and ``implicitfast``),
-     modify the mass matrix :math:`M \rightarrow M-hD`. This implies that finite-differenced accelerations
-     :math:`(v_{t+h} - v_t)/h` will not correspond to the continuous-time acceleration ``mjData.qacc``. When this flag
-     is enabled, :ref:`mj_inverse` will interpret ``qacc`` as having been computed from the difference of two sequential
-     velocities, and undo the above modification.
-
-   Disable midpoint integration
-     Additionally and relatedly, this flag **disables** :ref:`midpoint integration<geMidpoint>` for free bodies, which
-     would otherwise break the linear relationship between finite-differenced velocities and forces assumed by discrete
-     inverse dynamics. Note that disabling midpoint integration might be useful for debugging or for other reasons,
-     regardless or whether inverse dynamics are used.
+   This flag enables discrete-time inverse dynamics with :ref:`mj_inverse` for all
+   :ref:`integrators<option-integrator>` other than ``RK4``. Recall from the
+   :ref:`numerical integration<geIntegration>` section that the one-step integrators (``Euler``, ``implicit`` and
+   ``implicitfast``), modify the mass matrix :math:`M \rightarrow M-hD`. This implies that finite-differenced
+   accelerations :math:`(v_{t+h} - v_t)/h` will not correspond to the continuous-time acceleration ``mjData.qacc``.
+   When this flag is enabled, :ref:`mj_inverse` will interpret ``qacc`` as having been computed from the difference of
+   two sequential velocities, and undo the above modification.
 
 
 .. _option-flag-multiccd:
 
@@ -7,6 +7,14 @@ Upcoming version (not yet released)
 
 General
 ^^^^^^^
+- Replaced midpoint integration of free bodies with :ref:`gyroscopic derivatives<geFreeBody>` in the ``implicitfast``
+  :ref:`integrator<geIntegrators>`: the bias-force derivative of every standalone free body is applied via a local
+  unsymmetric solve of its decoupled block, making ``implicitfast`` identical to ``implicit`` for such bodies.
+  Unlike midpoint integration, which required vacuum and no constraints, this applies in all environments (contacts,
+  fluid, constraints), and is compatible with discrete-time inverse dynamics. Spinning free bodies no
+  longer gain energy, but tumbling motion is now mildly damped; models requiring long-horizon energy conservation of
+  tumbling bodies in vacuum should use ``RK4``. The :ref:`invdiscrete<option-flag-invdiscrete>` flag no longer has any
+  effect on forward dynamics.
 - Added Nesterov momentum extrapolation with adaptive gradient restart (O'Donoghue-Candès) to the PGS solver,
   significantly improving convergence. Overall PGS now requires ~2x fewer iterations.
 - Added the Newton decrement -- the quadratic model's predicted cost improvement of the next iteration -- as a third
@@ -233,7 +241,7 @@ General
    improving performance by ~20%.
 3. :commit:`b9c1877e` Added support for :ref:`elastic2d<flex-elasticity-elastic2d>` for trilinear and quadratic flex
    :ref:`dofs<body-flexcomp-dof>`.
-4. :commit:`910b3336` :ref:`Midpoint integration<geMidpoint>` is now restricted to the ``implicitfast``
+4. :commit:`910b3336` Midpoint integration is now restricted to the ``implicitfast``
    :ref:`integrator<geIntegrators>` and is disabled when fluid forces are active
    (nonzero :ref:`density<option-density>` or :ref:`viscosity<option-viscosity>`).
    Midpoint integration treats external forces as zero-order-hold constants, which causes
@@ -335,7 +343,7 @@ General
    scalar arrays (``jnt_stiffness``, ``dof_damping``, etc.) continue to hold the linear coefficient and are unchanged.
    The polynomial order is defined by the new constant :ref:`mjNPOLY<glNumericSizes>`. A future breaking C-API change
    may unify the linear and higher-order coefficients into a single array.
-4. :commit:`0c337799` Added :ref:`midpoint integration<geMidpoint>` for standalone free bodies in ``implicit`` and
+4. :commit:`0c337799` Added midpoint integration for standalone free bodies in ``implicit`` and
    ``implicitfast`` :ref:`integrators<geIntegrators>`. This applies the implicit midpoint rule to the rotational
    dynamics of free bodies with no children, conserving kinetic energy to machine precision in the absence of external
    torques. The :ref:`invdiscrete<option-flag-invdiscrete>` flag now also disables midpoint integration, providing an
 
@@ -580,50 +580,24 @@ Solving for :math:`v_{t+h}`, we obtain the implicit-in-velocity update
        \widehat{M} &\equiv M-h D
    \end{aligned}
 
-.. _geMidpoint:
+.. _geFreeBody:
 
-Midpoint integration for free bodies in vacuum
-   The implicit-in-velocity update :eq:`eq_implicit_update` treats the acceleration as a function of velocity and
-   linearizes. While effective for damping-like forces, it is sub-optimal for rotational dynamics, where
-   Coriolis and gyroscopic forces are *quadratic* in angular velocity. For this case, a better approach is to directly
-   discretize the rotational equations of motion using the *midpoint method*.
+Gyroscopic derivatives for free bodies
+   The ``implicitfast`` integrator described :ref:`below<geIntegrators>` excludes the derivatives of centripetal,
+   Coriolis and gyroscopic forces from :math:`D`, so that :math:`\widehat M` remains symmetric and can be factorized
+   with the faster Cholesky decomposition. However integrating gyroscopic forces explicitly can lead to
+   energy gain and divergence of fast-spinning free bodies with asymmetric inertia.
 
-   Consider a rigid body rotating in its principal-axis frame with angular velocity
-   :math:`\omega \in \mathbb{R}^3` and diagonal inertia tensor :math:`I = \text{diag}(I_1, I_2, I_3)`. The rotational
-   dynamics are given by `Euler's rotation equation
-   <https://en.wikipedia.org/wiki/Euler%27s_equations_(rigid_body_dynamics)>`__:
+   Therefore for *standalone free bodies* (free joints whose body has no children), these derivatives are reinstated.
+   The rows of :math:`\widehat M` corresponding to such a body form a :math:`6\times 6` block which is decoupled from
+   the rest of the system. After the global Cholesky solve, this block is re-assembled with the exact derivative of the
+   body's bias force and re-solved with an optimized :math:`6\times 6` LU routine. For standalone free bodies,
+   ``implicitfast`` and ``implicit`` therefore compute identical updates.
 
-   .. math::
-      I \dot{\omega} + \omega \times I\omega = \tau
-
-   where :math:`\tau` is the external torque in the principal-axis frame.
-   Evaluating the velocities at the midpoint, :math:`\omega_\text{mid} = (\omega_t + \omega_{t+h})/2`, gives:
-
-   .. math::
-      \frac{2}{h} I (\omega_\text{mid} - \omega_t) + \omega_\text{mid} \times I \omega_\text{mid} = \tau
-
-   This is a system of 3 nonlinear equations in 3 unknowns :math:`\omega_\text{mid}`, solved at each timestep using
-   Newton's method with a backtracking line search. After solving, the new velocity is recovered as
-   :math:`\omega_{t+h} = 2\omega_\text{mid} - \omega_t`.
-
-   **Properties.** The midpoint method preserves all `quadratic first integrals
-   <https://doi.org/10.1007/3-540-30666-8>`__ of the ODE. For Euler's equations, these are the
-   kinetic energy :math:`H = \frac{1}{2}\omega^T I\omega` and the squared angular momentum
-   :math:`C = \frac{1}{2}|I\omega|^2`, both conserved exactly in the absence of external torque. Since :math:`C` is the
-   Casimir function of the `Lie-Poisson <https://en.wikipedia.org/wiki/Poisson_bracket>`__ structure, the midpoint
-   method is a symmetric (time-reversible) and second-order accurate *Poisson integrator*.
-
-   **Eligibility.** Midpoint integration is only applied when using the ``implicitfast`` integrator, to
-   free bodies with no child bodies, and only when the medium has zero :ref:`density<option-density>` and
-   :ref:`viscosity<option-viscosity>`.
-
-   **Performance.** While the midpoint method carries computational overhead, we've found it to be
-   negligible compared to the rest of the pipeline, on the order of 1% in the worst case.
-
-   **Disabling.** Because midpoint integration solves a nonlinear equation for the next velocity, it breaks the linear
-   relationship between finite-differenced velocities and forces assumed by discrete inverse dynamics. Therefore,
-   setting the :ref:`invdiscrete<option-flag-invdiscrete>` flag disables midpoint integration, and also provides a
-   general opt-out mechanism for this integrator.
+   **Properties.** The kinetic energy of a spinning free body is non-increasing in the absence of applied force.
+   Steady spins about principal axes are conserved almost exactly; tumbling motion is mildly damped, at a rate
+   scaling like :math:`(h|\omega|)^2` per step. Systems requiring long-horizon energy conservation of tumbling
+   bodies should use the ``RK4`` integrator.
 
 .. _geIntegrators:
 
@@ -661,8 +635,8 @@ Fast implicit-in-velocity (``implicitfast``)
    scenarios which are not common and already well-handled by the Runge-Kutta integrator (see below). Because the RNE
    derivatives are also the main source of asymmetry of :math:`D`, by dropping them and symmetrizing, we can use the
    faster :math:`L^TL` rather than :math:`LU` decomposition.
-   The ``implicitfast`` integrator applies :ref:`midpoint integration<geMidpoint>` to eligible free bodies in vacuum,
-   providing exact energy conservation for spinning objects at negligible additional cost.
+   For standalone free bodies, the dropped :ref:`gyroscopic derivatives<geFreeBody>` are reinstated with a local
+   unsymmetric solve, preventing energy gain of spinning bodies at negligible additional cost.
 
 4th-order Runge-Kutta (``RK4``)
    One advantage of our continuous-time formulation is that we can use higher order integrators such as Runge-Kutta or
@@ -696,11 +670,10 @@ Fast implicit-in-velocity (``implicitfast``)
      increased stability, and is therefore a strict improvement. It is the recommended integrator for most models.
     **implicit**:
      The benefit over ``implicitfast`` is the implicit integration of Coriolis and centripetal forces for *coupled*
-     rotational systems such as multi-link pendula. Note that ``implicit`` does not apply :ref:`midpoint
-     integration<geMidpoint>` (only ``implicitfast`` does), but its RNE derivatives provide comparable stability
-     for free-body rotation. For example, `gyroscopic.xml <../_static/gyroscopic.xml>`__ shows an ellipsoid rolling
-     on an inclined plane; both ``implicitfast`` and ``implicit`` handle this case well, while ``Euler`` quickly
-     diverges.
+     rotational systems such as multi-link pendula. For standalone free bodies the two integrators coincide, since
+     ``implicitfast`` applies the :ref:`gyroscopic derivatives<geFreeBody>` to such bodies. For example,
+     `gyroscopic.xml <../_static/gyroscopic.xml>`__ shows an ellipsoid rolling on an inclined plane; both
+     ``implicitfast`` and ``implicit`` handle this case well, while ``Euler`` quickly diverges.
     **RK4**:
      This integrator is best for systems which are energy conserving, or almost energy-conserving. `pendulum.xml
      <../_static/pendulum.xml>`__ shows a complicated pendulum mechanism which diverges quickly using ``Euler`` or
 
@@ -19,6 +19,7 @@
 #include <mujoco/mjsan.h>  // IWYU pragma: keep
 #include "engine/engine_core_util.h"
 #include "engine/engine_crossplatform.h"
+#include "engine/engine_inline.h"
 #include "engine/engine_memory.h"
 #include "engine/engine_passive.h"
 #include "engine/engine_sleep.h"
@@ -705,6 +706,186 @@ static void mjd_rne_vel(const mjModel* m, mjData* d) {
 }
 
 
+// 3x3 sub-blocks of (d qfrc_bias / d qvel) for a standalone free body
+//   outputs the two 3x3 blocks lin and rot such that the rotational columns
+//   of the full 6x6 bias Jacobian B are  [-mass*lin; rot]  (linear columns are zero)
+//
+// derivation: let R = xmat, s = xipos - xpos, w = R*qvel[rot] (world angular velocity),
+// Iw = ximat * diag(body_inertia) * ximat' (world inertia about the CoM). with qacc = 0,
+// the CoM acceleration is w x (w x s) and the world bias force/torque at the CoM are
+//   f = mass * w x (w x s),   tau = w x Iw*w
+// projected onto the joint coordinates: bias = [f;  R'*(s x f + tau)]. differentiating
+// w.r.t. the rotational dofs (through w = R*qvel[rot]), with K = [w x s]_x + [w]_x [s]_x:
+//   d f / d w   = -mass * K            =>  lin = K * R
+//   d tau / d w = [w]_x Iw - [Iw*w]_x  =>  rot = R' * (-mass*[s]_x K + d tau/d w) * R
+static void freeBias_vel_blocks(mjtNum mass, const mjtNum R[9], const mjtNum Xi[9],
+                                const mjtNum inertia[3], const mjtNum s[3],
+                                const mjtNum qvel_rot[3], mjtNum lin[9], mjtNum rot[9]) {
+  // world-frame angular velocity
+  mjtNum w[3];
+  mji_mulMatVec3(w, R, qvel_rot);
+
+  // world-frame inertia about CoM: Iw = Xi * diag(inertia) * Xi^T
+  mjtNum Xi_I[9];
+  for (int i=0; i < 3; i++) {
+    Xi_I[3*i+0] = Xi[3*i+0] * inertia[0];
+    Xi_I[3*i+1] = Xi[3*i+1] * inertia[1];
+    Xi_I[3*i+2] = Xi[3*i+2] * inertia[2];
+  }
+  mjtNum Iw[9];
+  Iw[0] = Xi_I[0]*Xi[0] + Xi_I[1]*Xi[1] + Xi_I[2]*Xi[2];
+  Iw[4] = Xi_I[3]*Xi[3] + Xi_I[4]*Xi[4] + Xi_I[5]*Xi[5];
+  Iw[8] = Xi_I[6]*Xi[6] + Xi_I[7]*Xi[7] + Xi_I[8]*Xi[8];
+  Iw[1] = Iw[3] = Xi_I[0]*Xi[3] + Xi_I[1]*Xi[4] + Xi_I[2]*Xi[5];
+  Iw[2] = Iw[6] = Xi_I[0]*Xi[6] + Xi_I[1]*Xi[7] + Xi_I[2]*Xi[8];
+  Iw[5] = Iw[7] = Xi_I[3]*Xi[6] + Xi_I[4]*Xi[7] + Xi_I[5]*Xi[8];
+
+  // intermediate vectors: ws = w x s  (CoM offset velocity),  Iww = Iw * w  (angular momentum)
+  mjtNum ws[3], Iww[3];
+  mji_cross(ws, w, s);
+  mji_mulMatVec3(Iww, Iw, w);
+
+  // K = [w x s]_x + [w]_x [s]_x = s w^T - (w . s) I + [ws]_x
+  mjtNum w_dot_s = w[0]*s[0] + w[1]*s[1] + w[2]*s[2];
+  mjtNum K[9];
+  K[0] = s[0]*w[0] - w_dot_s;
+  K[1] = s[0]*w[1] - ws[2];
+  K[2] = s[0]*w[2] + ws[1];
+
+  K[3] = s[1]*w[0] + ws[2];
+  K[4] = s[1]*w[1] - w_dot_s;
+  K[5] = s[1]*w[2] - ws[0];
+
+  K[6] = s[2]*w[0] - ws[1];
+  K[7] = s[2]*w[1] + ws[0];
+  K[8] = s[2]*w[2] - w_dot_s;
+
+  // lin = K * R
+  mji_mulMatMat3(lin, K, R);
+
+  // C = -mass * [s]_x K + [w]_x Iw - [Iww]_x,  column by column
+  // the last term (-[Iww]_x) is the negated cross-product matrix, added via ternaries
+  mjtNum C[9];
+  for (int c=0; c < 3; c++) {
+    mjtNum s_x_K_row0 = s[1]*K[6+c] - s[2]*K[3+c];
+    mjtNum s_x_K_row1 = s[2]*K[c] - s[0]*K[6+c];
+    mjtNum s_x_K_row2 = s[0]*K[3+c] - s[1]*K[c];
+
+    mjtNum w_x_Iw_row0 = w[1]*Iw[6+c] - w[2]*Iw[3+c];
+    mjtNum w_x_Iw_row1 = w[2]*Iw[c] - w[0]*Iw[6+c];
+    mjtNum w_x_Iw_row2 = w[0]*Iw[3+c] - w[1]*Iw[c];
+
+    C[c]     = -mass * s_x_K_row0 + w_x_Iw_row0 + (c == 1 ? Iww[2] : (c == 2 ? -Iww[1] : 0));
+    C[3 + c] = -mass * s_x_K_row1 + w_x_Iw_row1 + (c == 0 ? -Iww[2] : (c == 2 ? Iww[0] : 0));
+    C[6 + c] = -mass * s_x_K_row2 + w_x_Iw_row2 + (c == 0 ? Iww[1] : (c == 1 ? -Iww[0] : 0));
+  }
+
+  // rot = R^T * C * R
+  mjtNum tmp[9];
+  mji_mulMatTMat3(tmp, R, C);
+  mji_mulMatMat3(rot, tmp, R);
+}
+
+
+// 6x6 block B = d qfrc_bias / d qvel for a standalone free body
+//   assembles the full 6x6 from the 3x3 sub-blocks computed by freeBias_vel_blocks
+//   rows/cols ordered like the free joint dofs: [linear(3); rotational(3)]
+//   linear columns are zero: the bias force does not depend on linear velocity
+void mjd_freeBias_vel(const mjModel* m, const mjData* d, int jnt, mjtNum B[36]) {
+  int body = m->jnt_bodyid[jnt];
+  int adr = m->jnt_dofadr[jnt];
+  mjtNum mass = m->body_mass[body];
+  const mjtNum* R = d->xmat + 9*body;    // body  -> world
+  const mjtNum* Xi = d->ximat + 9*body;  // inertia -> world
+  const mjtNum* inertia = m->body_inertia + 3*body;
+
+  // CoM offset from joint origin, world frame
+  mjtNum s[3];
+  mji_sub3(s, d->xipos + 3*body, d->xpos + 3*body);
+
+  mjtNum lin[9], rot[9];
+  freeBias_vel_blocks(mass, R, Xi, inertia, s, d->qvel + adr + 3, lin, rot);
+
+  mju_zero(B, 36);
+  for (int r=0; r < 3; r++) {
+    for (int c=0; c < 3; c++) {
+      B[6*r + 3+c] = -mass * lin[3*r+c];
+      B[6*(3+r) + 3+c] = rot[3*r+c];
+    }
+  }
+}
+
+
+// 6x6 block A = M - h * (d qfrc_smooth / d qvel) for the free joint of a standalone body
+//   returns 1 and writes A if jnt is the free joint of a standalone awake body, 0 otherwise
+//   requires valid d->qDeriv rows for the block, computed with flg_bias = 0; the bias
+//   derivative excluded from qDeriv is added here via freeBias_vel_blocks
+int mjd_freeMhat(const mjModel* m, const mjData* d, int jnt, mjtNum h, mjtNum A[36]) {
+  // must be a free joint
+  if (m->jnt_type[jnt] != mjJNT_FREE) {
+    return 0;
+  }
+
+  int body = m->jnt_bodyid[jnt];
+  int adr = m->jnt_dofadr[jnt];
+  int tree = m->dof_treeid[adr];
+  mjtNum mass = m->body_mass[body];
+
+  // must be a standalone 6-DOF tree with no children, awake
+  if (m->tree_dofnum[tree] != 6 ||
+      m->body_subtreemass[body] != mass ||
+      !d->tree_awake[tree]) {
+    return 0;
+  }
+
+  // D rows of a standalone free body are exactly the 6x6 block (D sparsity is tree-local);
+  // guard the gathers below against any violation of this invariant
+  if (m->D_rownnz[adr] != 6) {
+    return 0;
+  }
+
+  // A = M block (gather from sparse lower triangle)
+  mju_zero(A, 36);
+  for (int r=0; r < 6; r++) {
+    int rowadr = m->M_rowadr[adr+r];
+    int rownnz = m->M_rownnz[adr+r];
+    for (int k=0; k < rownnz; k++) {
+      int c = m->M_colind[rowadr+k] - adr;
+      A[6*r+c] = A[6*c+r] = d->M[rowadr+k];
+    }
+  }
+
+  // A -= h * qDeriv block (actuator and passive derivatives)
+  for (int r=0; r < 6; r++) {
+    int rowadr = m->D_rowadr[adr+r];
+    int rownnz = m->D_rownnz[adr+r];
+    for (int k=0; k < rownnz; k++) {
+      int c = m->D_colind[rowadr+k] - adr;
+      A[6*r+c] -= h * d->qDeriv[rowadr+k];
+    }
+  }
+
+  // A -= h * d(qfrc_smooth)/d(qvel) for the bias term missing from qDeriv;
+  // qfrc_smooth includes -qfrc_bias, so subtracting its derivative adds +h*B
+  mjtNum s[3];
+  mji_sub3(s, d->xipos + 3*body, d->xpos + 3*body);
+
+  mjtNum lin[9], rot[9];
+  freeBias_vel_blocks(mass, d->xmat + 9*body, d->ximat + 9*body,
+                      m->body_inertia + 3*body, s, d->qvel + adr + 3, lin, rot);
+
+  mjtNum h_mass = -h * mass;
+  for (int r=0; r < 3; r++) {
+    for (int c=0; c < 3; c++) {
+      A[6*r + 3+c] += h_mass * lin[3*r+c];
+      A[6*(3+r) + 3+c] += h * rot[3*r+c];
+    }
+  }
+
+  return 1;
+}
+
+
 //--------------------- utility functions for (d force / d vel) Jacobians --------------------------
 
 // add J'*B*J to qDeriv
 
@@ -43,6 +43,15 @@ MJAPI void mjd_passive_vel(const mjModel* m, mjData* d);
 // subtract (d qfrc_bias / d qvel) from qDeriv (dense version)
 MJAPI void mjd_rne_vel_dense(const mjModel* m, mjData* d);
 
+// 6x6 block B = d qfrc_bias / d qvel for the free joint of a standalone body
+MJAPI void mjd_freeBias_vel(const mjModel* m, const mjData* d, int jnt,
+                            mjtNum B[36]);
+
+// 6x6 block A = M - h * (d qfrc_smooth / d qvel) for the free joint of a standalone body
+//   returns 1 and writes A if jnt is the free joint of a standalone awake body, 0 otherwise
+//   requires valid d->qDeriv rows for the block, computed with flg_bias = 0
+MJAPI int mjd_freeMhat(const mjModel* m, const mjData* d, int jnt, mjtNum h, mjtNum A[36]);
+
 // compute res += (s1 + s2*damping) * J'*K*J * vec, for all interpolated flexes
 //   K_rot_cache: if non-NULL, use pre-cached K_rot (same layout as m->flex_stiffness)
 MJAPI void mjd_flexInterp_mul(const mjModel* m, mjData* d, mjtNum* res, const mjtNum* vec,