H.glsl

vec4
qconj(vec4 x)
{
	return vec4(-x.xyz, x.w);
}

vec4
qmul(vec4 x, vec4 y)
{
	return vec4(
		x.w * y.xyz + y.w * x.xyz + cross(x.xyz, y.xyz),
		x.w * y.w - dot(x.xyz, y.xyz));
}

vec4
qinv(vec4 x)
{
	return qconj(x) / dot(x, x);
}

/*
Inverse of a unit quaternion
*/
vec4
qinvu(vec4 x)
{
	return qconj(x);
}

vec4
qlog(vec4 x)
{
	float nrm = length(x);
	float vnrm = length(x.xyz);
	if (vnrm == 0.0) {
		return vec4(vec3(0.0), log(nrm));
	} else {
		return vec4((acos(x.w / nrm) / vnrm) * x.xyz, log(nrm));
	}
}

/*
Logarithm of a unit quaternion
*/
vec4
qlogu(vec4 x)
{
	float vnrm = length(x.xyz);
	if (vnrm == 0.0) {
		return vec4(0.0);
	} else {
		return vec4((acos(x.w) / vnrm) * x.xyz, 0.0);
	}
}

vec4
qexp(vec4 x)
{
	float vnrm = length(x.xyz);
	if (vnrm == 0.0) {
		return vec4(vec3(0.0), 1.0);
	} else {
		float wexp = exp(x.w);
		return vec4((wexp * sin(vnrm) / vnrm) * x.xyz, wexp * cos(vnrm));
	}
}

/*
Exponential of an imaginary quaternion
*/
vec4
qexpv(vec4 x)
{
	float vnrm = length(x.xyz);
	if (vnrm == 0.0) {
		return vec4(vec3(0.0), 1.0);
	} else {
		return vec4((sin(vnrm) / vnrm) * x, cos(vnrm));
	}
}

vec4
qpow(vec4 x, float y)
{
	return qexp(y * qlog(x));
}

/*
Rotation
*/
vec3
qvrot(vec4 x, vec3 y)
{
	return y + 2.0 * cross(x.xyz, x.w * y + cross(x.xyz, y));
}

vec4[2]
dmul(vec4 x[2], vec4 y[2])
{
	return vec4[2](qmul(x[0], y[0]), qmul(x[0], y[1]) + qmul(x[1], y[0]));
}

/*
Inverse of a unit dual quaternion
*/
vec4[2]
dinvu(vec4 x[2])
{
	return vec4[2](qconj(x[0]), qconj(x[1]));
}

vec4[2]
dlog(vec4 x[2])
{
	vec4 rlog = qlog(x[0]);
	return vec4[2](rlog, qmul(qinvu(rlog), x[1]) / dot(x[0], x[0]));
}

/*
Logarithm of a unit dual quaternion
*/
vec4[2]
dlogu(vec4 x[2])
{
	vec4 rlog = qlogu(x[0]);
	return vec4[2](rlog, qmul(qinvu(rlog), x[1]));
}

vec4[2]
dexp(vec4 x[2])
{
	vec4 rexp = qexp(x[0]);
	return vec4[2](rexp, qmul(rexp, x[1]));
}

vec4[2]
dpow(vec4 x[2], float y)
{
	vec4 xlog[2] = dlog(x);
	return dexp(vec4[2](y * xlog[0], y * xlog[1]));
}

vec4[2]
sclerp(vec4 x[2], vec4 y[2], float t)
{
	return dmul(dpow(dmul(y, dinvu(x)), t), x);
}

/*
Rotation and displacement
*/
vec3
dvdisp(vec4 x[2], vec3 y)
{
	return qvrot(x[0], y)
		+ 2.0 * (x[0].w * x[1].xyz - x[1].w * x[0].xyz + cross(x[0].xyz, x[1].xyz));
}