Numerically stable method for angle between 3D vectors

I'm looking for a numerically stable method for computing the angle between two 3D vectors. Which of the following methods ought to be preferred?

Method 1: $$ u\times v = ||u||~||v|| \sin(\theta) \textbf{n}\\ u\cdot v = ||u||~||v|| \cos(\theta)\\ \theta = \arctan2(||u\times v||,~u\cdot v) $$

Method 2:

$$ \theta = 2~\arctan2(||u/||u|| - v/||v||~||,~||u/||u|| + v/||v||~||) $$

where method 2 is based on the fact that the sum and difference vectors of two unit (or equal length) vectors are orthogonal.

Or is there an even better method which I haven't thought of? The vectors I am considering have angles typically smaller than 2 degrees.

Short answer: Method 2 is better for small angles. Use this slight rearrangement: $$ \theta = 2~atan2(||~||v||u - ||u||v~||,~||~||v||u + ||u||v~||) $$

This formula comes from W. Kahan's advice in his paper "How Futile are Mindless Assessments of Roundoff in Floating-Point Computation?" (https://www.cs.berkeley.edu/~wkahan/Mindless.pdf), section 12 "Mangled Angles."

Before I found that paper, I also did my own analysis, comparing the double-precision results with result's using bc with 50 digits of precision. Method 2 is commonly within 1 ulp of correct, and almost always within 50 ulps, whereas Method 1 was much, much less accurate.

This might seem surprising, since Method 1 is mathematically insensitive to the magnitude of u and v, whereas they have to be normalized in Method 2. And indeed, the accuracy of the normalization is the limiting factor for Method 2 - virtually all the error comes from the fact that even after normalization, the vectors aren't exactly length 1.

However, for small angles you get catastrophic cancellation in the cross product for method 1. Specifically, all the products like $u_y v_z - u_z v_y$ end up close to 0, and I believe the cross-multiplication before the subtraction loses accuracy, compared to the direct subtraction in Method 2.