title: Calculating the angle between two vectors using a needle-like triangle
body: <p>I implemented a function (<code>angle_between</code>) to calculate the angle between two vectors. It makes use of needle-like triangles and is based on <a href="https://people.eecs.berkeley.edu/%7Ewkahan/Triangle.pdf" rel="nofollow noreferrer">Miscalculating Area and Angles of a Needle-like Triangle</a> and <a href="https://scicomp.stackexchange.com/questions/27689/numerically-stable-way-of-computing-angles-between-vectors">this related question</a>.</p>
  <p>The function appears to work fine most of the time, except for one weird case where I don't understand what is happening:</p>
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  vectorA = np.array([0.008741225033460295, 1.1102230246251565e-16], dtype=np.float64)
  vectorB = np.array([1, 0], dtype=np.float64)
  angle_between(vectorA, vectorB)  # is np.nan
  </code></pre>
  <p>Digging into my function, the <code>np.nan</code> is produced by taking the square root of a negative number, and the negative number seems to be the result of the increased accuracy of the method:</p>
  <pre class="lang-py prettyprint-override"><code>foo = 1.0                  # np.linalg.norm(vectorA)
  bar = 0.008741225033460295 # np.linalg.norm(vectorB)
  baz = 0.9912587749665397   # np.linalg.norm(vectorA- vectorB)
  
  # algebraically equivalent ... numerically not so much
  order1 = baz - (foo - bar)
  order2 = bar - (foo - baz)
  
  assert order1 == 0
  assert order2 == -1.3877787807814457e-17
  </code></pre>
  <p>According to Kahan's paper, this means that the triplet (foo, bar, baz) actually doesn't represent the side lengths of a triangle. However, this should - in fact - be the case given how I constructed the triangle (see the comments in the code).</p>
  <p>From here, I feel a bit lost as to where to look for the source of the error. Could somebody explain to me what is happening?</p>
  <hr />
  <p>For completeness, here is the full code of my function:</p>
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  from numpy.typing import ArrayLike
  
  def angle_between(
      vec_a: ArrayLike, vec_b: ArrayLike, *, axis: int = -1, eps=1e-10
  ) -&gt; np.ndarray:
      &quot;&quot;&quot;Computes the angle from a to b
  
      Notes
      -----
      Implementation is based on this post:
      https://scicomp.stackexchange.com/a/27694
      &quot;&quot;&quot;
  
      vec_a = np.asarray(vec_a)[None, :]
      vec_b = np.asarray(vec_b)[None, :]
  
      if axis &gt;= 0:
          axis += 1
  
      len_c = np.linalg.norm(vec_a - vec_b, axis=axis)
      len_a = np.linalg.norm(vec_a, axis=axis)
      len_b = np.linalg.norm(vec_b, axis=axis)
  
      mask = len_a &gt;= len_b
      tmp = np.where(mask, len_a, len_b)
      np.putmask(len_b, ~mask, len_a)
      len_a = tmp
  
      mask = len_c &gt; len_b
      mu = np.where(mask, len_b - (len_a - len_c), len_c - (len_a - len_b))
  
      numerator = ((len_a - len_b) + len_c) * mu
      denominator = (len_a + (len_b + len_c)) * ((len_a - len_c) + len_b)
  
      mask = denominator &gt; eps
      angle = np.divide(numerator, denominator, where=mask)
      np.sqrt(angle, out=angle)
      np.arctan(angle, out=angle)
      angle *= 2
      np.putmask(angle, ~mask, np.pi)
      return angle[0]
  </code></pre>
  <p><strong>Edit:</strong> The problem is definitely related to <code>float64</code> and disappears when performing the computation with larger floats:</p>
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  
  vectorA = np.array([0.008741225033460295, 1.1102230246251565e-16], dtype=np.float128)
  vectorB = np.array([1, 0], dtype=np.float128)
  assert angle_between(vectorA, vectorB) == 0
  </code></pre>