void mul(Quat4d & q1, Quat4d & q2, Quat4d & q3) {
q3.w = q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z;
q3.x = q1.x * q2.w + q1.w * q2.x - q1.z * q2.y + q1.y * q2.z;
q3.y = q1.y * q2.w + q1.w * q2.y - q1.x * q2.z + q1.z * q2.x;
q3.z = q1.z * q2.w + q1.w * q2.z - q1.y * q2.x + q1.x * q2.y;
}
void mul(Point3D & v, Quat4d & q2, Quat4d & qout) {
qout.w = -v.x * q2.x - v.y * q2.y - v.z * q2.z;
qout.x = v.x * q2.w - v.z * q2.y + v.y * q2.z;
qout.y = v.y * q2.w - v.x * q2.z + v.z * q2.x;
qout.z = v.z * q2.w - v.y * q2.x + v.x * q2.y;
}
void rotate(Quat4d & q1, Point3D & v, Point3D & vout) {
w2 = -v.x * q2.x - v.y * q2.y - v.z * q2.z;
x2 = -v.x * q2.w + v.z * q2.y - v.y * q2.z;
y2 = -v.y * q2.w + v.x * q2.z - v.z * q2.x;
z2 = -v.z * q2.w + v.y * q2.x - v.x * q2.y;
vout.x = q1.x * w2 + q1.w * x2 - q1.z * y2 + q1.y * z2;
vout.y = q1.y * w2 + q1.w * y2 - q1.x * z2 + q1.z * x2;
vout.z = q1.z * w2 + q1.w * z2 - q1.y * x2 + q1.x * y2;
}
Uses an equivalent formulation but with cross products