Add detailed comments for quaternion multiplication

Author: flferretti

src/jaxsim/math/quaternion.py

@@ -87,7 +87,17 @@ def derivative(
         # Construct pure quaternion: (scalar damping term, angular velocity components)
         ω_quat = jnp.hstack([K * safe_norm(ω) * (1 - safe_norm(quaternion)), ω])
 
-        # Apply quaternion multiplication based on frame representation
+        # Quaternion multiplication using index tables.
+        # This approach avoids using the explicit quaternion multiplication formula
+        # by encoding the necessary element-wise products and signs via indexed operations.
+        # Given two quaternions q and w, their Hamilton product q ⊗ w can be written
+        # as a combination of q[i] * w[j] terms with appropriate signs.
+        # i_idx and j_idx define which elements of the outer product q ⊗ w to select.
+        # For example, i_idx[1][2] = 2 and j_idx[1][2] = 3 means: take q[2] * w[3] for this term.
+        # sign_matrix[i][j] gives the sign (+1 or -1) to apply to each q[i] * w[j] term,
+        # depending on quaternion multiplication rules.
+        # This indexed summation reproduces the Hamilton product of quaternions in a
+        # vectorized way, and is suitable for use with JAX.
         i_idx = jnp.array([[0, 1, 2, 3], [0, 1, 2, 3], [0, 2, 3, 1], [0, 3, 1, 2]])
         j_idx = jnp.array([[0, 1, 2, 3], [1, 0, 3, 2], [2, 0, 1, 3], [3, 0, 2, 1]])
         sign_matrix = jnp.array(