Can be a little tedious on the fingers (actually, who am I kidding - I rarely modify stuff that it's a pleasure to dip back into a simple set of interfaces now and again).

How I like my C++ math interfaces

// ----- math.h --------------------------------------

namespace math
{
  // Simple types with exposed data members, no constructors and multiple addressing modes
  // No templates in the public API
  
	struct quatf
	{
		union
		{
			struct { float x, y, z, w; };
			vec3f v;
			float f[4];
		};
	};
	
	struct quatd
	{
		union
		{
			struct { double x, y, z, w; };
			vec3d v;
			double f[4];
		};
	};
	
	// Procedural operations on types as free functions
	quatf qfMultiply(const quatf& a, const quatf& b);
	quatd qdMultiply(const quatd& a, const quatd& b);
}


// ----- math.cpp ------------------------------------


// Hidden templated implementation with all your wonderful workings and comments
template <typename QUAT_TYPE>
QUAT_TYPE qMultiply(const QUAT_TYPE& a, const QUAT_TYPE& b)
{
	//
	// Multiplication is simply the expanded product of:
	//
	//    (p0 + i.p1 + j.p2 + k.p3)(q0 + i.q1 + j.q2 + k.q3)
	//
	// where:
	//
	//    i^2 = j^2 = k^2 = ijk = -1,
	//    ij = k, ji = -k,
	//    jk = i, kj = -i,
	//    ki = j, ik = -j
	//
	// This results in the following:
	//
	// pq = p0q0  + i.p0q1   + j.p0q2   + k.p0q3 +
	//     i.p1q0 + i.i.p1q1 + i.j.p1q2 + i.k.p1q3 +
	//     j.p2q0 + j.i.p2q1 + j.j.p2q2 + j.k.p2q3 +
	//     k.p3q0 + k.i.p3q1 + k.j.p3q2 + k.k.p3q3
	//
	// pq = p0q0   + i.p0q1 + j.p0q2 + k.p0q3 +
	//      i.p1q0 - p1q1   + k.p1q2 - j.p1q3 +
	//      j.p2q0 - k.p2q1 - p2q2   + i.p2q3 +
	//      k.p3q0 + j.p3q1 - i.p3q2 - p3.q3
	//
	// pq = p0q0   - p1q1 - p2q2 - p3q3 +
	//      i(p0q1 + p1q0 + p2q3 - p3q2)
	//      j(p0q2 - p1q3 + p2q0 + p3q1)
	//      k(p0q3 + p1q2 - p2q1 + p3q0)
	//
	// 16 muls, 12 add/sub
	//
	QUAT_TYPE q =
	{
		a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
		a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
		a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
		a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
	};

	return q;
}


// Implementations of procedural implementations
// Fully-scoped with namespaces to catch function prototype mis-matches at compile-time
// Also allows you to see namespace ownership of functions at a glance
math::quatf math::qfMultiply(const math::quatf& a, const math::quatf& b)
{
	return ::qMultiply(a, b);
}
math::quatd math::qdMultiply(const math::quatd& a, const math::quatd& b)
{
	return ::qMultiply(a, b);
}