I'm working on an Inverse Kinematic system using Jacobian, but I am struggling with creating a Jacobian matirx populated by partial derivatives for/with quaternions. At the moment I am using Euler angles using the current formula:

where V is local x,y or z axis transformed into world space. Here is a snippet of the code for that:

    glm::vec3 forward = transforms[parent].forward();
    glm::vec3 forwardCross = glm::cross(bone, forward);
    ikt.jacobi(index + 0, j + 2) = forwardCross.x;
    ikt.jacobi(index + 1, j + 2) = forwardCross.y;
    ikt.jacobi(index + 2, j + 2) = forwardCross.z;

And after doing Damped Least Squares I add the Theta changes through the eulers.

glm::vec3 oldEul = transforms[child].getEulerRotation();
glm::vec3 eulDelta = { dTheta(jointNum + 0), dTheta(jointNum + 1), dTheta(jointNum + 2) };
transforms[child].setEulerRotation(oldEul + eulDelta);

This method is mostly functional but falls apart approaching singularities/gimble lock.

Does anyone know a method to instead use quaternions to avoid this?