where r unc n+1 is the position after the unconstrained time step n + 1, M is the diagonal matrix of particle masses, B n is the gradient matrix of the constraint equations at time step n, and the vector d contains the prescribed length of the constraints. Eq 1 acts by removing the projection of B n r unc n+1 -d onto the constraints from an unconstrained update. A key step in enforcing the constraints with eq 1 is the inversion of B n M -1 B T n , a square matrix of dimensions K × K with K being the number of distance constraints. To reduce the computational cost, one rewrites the inverse in the form