well you could check long double for it
just make something lie this
rewrite dot4 so it will get quaterions (like old did but i will operate on long doubles and will return long double)
ofc function returns long double
then return 2.0 * acosl( fabsl( DotDouble4( qa, qb ) ) ); }
this could use the best precision for floats
anyway try to rotate something by in ex 90 degrees and try to find that angle.