Btw. this all leads after some trigonometric transformations also to the following quite efficient unit quaternion distance calculation, as I developed recently for the orientation spline estimator evaluation. It would have lead directly to this dot product check, after the same isotonic function argument:

function [ dist ] = unitLengthQuaternionRowwiseDist( qs1, qs2 )
dist = 2 * acos(min(1, max(-1, dot(qs1, qs2))));
end

This is much more efficient than the way of involving a full quaternion log on the fraction q0/q1. The clamping to [-1, 1] is of course only for numerical reasons. It was guarantied by the Cauchy-Scharz inequality if we had proper unit length quaternions.
(can be easily, directly inferred from the logarithm formula).