After some more investigation, most other libraries I've looked at handle this issue (of very similar quats, or non-unit quats) in the slerp() method. And they just accept that normalize() may not return truly unit quats.

See http://cgkit.svn.sourceforge.net/viewvc/cgkit/cgkit/trunk/supportlib/include/quat.h?revision=341&content-type=text%2Fplain (Line 187)

acos(ca) results in a complex number for values of ca > 1 or ca < -1. I've rewritten my slerp function as follows to handle this better:

#This is a Python version (but the same would work for C)

def slerp(t,q0,q1,shortest=True):

neg_q1 = False #Does q1 have to be negated (so that the shortest path is taken)?

ca = q0.dot(q1)

if (shortest and ca < 0):

ca = -ca

neg_q1 = True

if ( (1.0 - ca) > 1e-12): #epsilion?

#standard case

o = math.acos(ca)

so = math.sin(o)

a = math.sin(o*(1.0-t)) / so

b = math.sin(o*t) / so

else:

#quats are very close, do linear interpolation

a = 1.0 - t

b = t

#do the interpolation

if (neg_q1):

return q0*a - q1*b

else:

return q0*a + q1*b