@simonlynen : I think the squared norm is correct:

We want to see which pair is closer (q0, q1) or (q0, -q1) in terms of the distance in S^3 (=length of a shortest connecting curve between them (=: a 'shortest') ). Paul's solution is based on the assumption that this closer - pre-order on pairs in S^3 is identical to the pre-order given by the square norm distance of R^4. 
Or do we disagree already in that?

If not: Consider a pair (a, b) \in (S^3)^2. The shortest are always geodesics. It is a well known theorem that the geodesics of S^N are sections of great circles. This implies that all connecting geodesics and with them the shortest are contained in the plain determined by the 0, a, b. 
Having projected the whole comparison to a plain we can 2-d elemental geometry: a, b lie on a common unit circle (exactly the one great circle of S^3 containing all geodesic connections). 
Let's look at the connecting secant of (a,b). It is the strait line in R^4/our plain connecting the pair having the length L:= sqrt(squaredNorm(a - b)) = d_{R^4}(a,b). After drawing the radius lines to a, b we see that the resulting triangle yields the equation : L = 2 * sin (phi/2) where phi is the inner angle at the origin and thus exactly the length of the shorter circle section connecting a, b which is then just the distance in S^3. So we have :

(2 * sin (d_{S^3}(a,b) / 2) )^2 = squaredNorm(a - b)

As the map phi |-> (2*sin(phi / 2)^2) is strictly monotone increasing on phi\in [0, pi] we have proven the equality of the two pre - orders. (This is because any two pre - orders given by pulling back the order of the real numbers through two mappings (f1, f2) are equal when they are translatable by a strictly monotone increasing function g (f1 = g \circ f2). Because pulling back the order with f1 is according to that the same as pulling back first with g on R again and the the resulting order back with f2. But g is strictly monotone increasing on R which is identical for it and its inverse to be order preserving on R. So pulling back with g did not change the order on R!)

Having discussed that I want to suggest the following optimizations : 
Consider: squaredNorm(a,b) = (a - b)^T(a-b) = 1 - 2 a^T*b + 1.
This gives us:  squaredNorm(q0 - q1) > squaredNorm(q0 + q2) <=> 2q0^T q1 > -2q0^T q1 <=> 4 q0^Tq1 > 0 <=> q0^T q1 > 0
(this result is also elemental geometric obvious if you look at it the right way)

@furgalep : This would be faster for your check.