An alternative formulation of the popular slerp function is:

(http://en.wikipedia.org/wiki/Slerp#Quaternion_Slerp)

my question is how is this actually computed in practice? Is it done by first finding the polar form and then multiplying the cosine and sine? I believe this is somewhat similar to how you would do it for complex numbers.

Managed to figure it out, also thanks to boost for their exp function

	// Raise a quaternion to a real power
	template <typename S>
	inline quaternion<S> pow(const quaternion<S> &q, const S ex)
	{
		return exp(ex * log(q) );
	}

	// Raise a unit quaternion to a real power
	template <typename S>
	inline quaternion<S> unit_pow(const quaternion<S> &q, const S ex)
	{
		static const S C = S(1) - std::numeric_limits<S>::epsilon();

		// Check for the case of an identity quaternion.
		// This will protect against divide by zero
		if ( std::abs(real(q)) > C ) 
			return q;

		S a = acos( real(q) );
		S u = a*ex;

		return quaternion<S>(std::cos(u), 
			imag(q) * (std::sin(u) / std::sin(a)));
	}

	// logarithm of a quaternion.
	// when q is unit this gives { 0 + acos(r) * v/|v| }
	template <typename S>
	inline quaternion<S> log ( const quaternion<S> &q )
	{	
		S N = std::sqrt(norm(q));

		return quaternion<S>( std::log(N), 
			std::acos(real(q) / N) * normalize(imag(q)) );	
	}

	// exponential of a quaternion
	template <typename S>
	inline quaternion<S> exp( const quaternion<S> &q )
	{
		S z = length(imag(q));
		S w = std::sin(z) / z;

		return std::exp(real(q)) * 
			quaternion<S>(std::cos(z), w * imag(q));
	}