使用针状三角形计算两个向量之间的角度

问：

我实现了一个函数（）来计算两个向量之间的角度。它利用针状三角形，并基于针状三角形的错误计算面积和角度以及这个相关问题。angle_between

该函数似乎大部分时间都工作正常，除了一个奇怪的情况，我不明白发生了什么：

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

深入研究我的函数，它是通过取负数的平方根产生的，负数似乎是该方法精度提高的结果：np.nan

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

根据 Kahan 的论文，这意味着三元组（foo、bar、baz）实际上并不代表三角形的边长。但是，考虑到我如何构造三角形，这应该是这种情况（请参阅代码中的注释）。

从这里开始，我对在哪里寻找错误的来源感到有点迷茫。有人可以向我解释发生了什么吗？

为了完整起见，以下是我的函数的完整代码：

import numpy as np
from numpy.typing import ArrayLike

def angle_between(
    vec_a: ArrayLike, vec_b: ArrayLike, *, axis: int = -1, eps=1e-10
) -> np.ndarray:
    """Computes the angle from a to b

    Notes
    -----
    Implementation is based on this post:
    https://scicomp.stackexchange.com/a/27694
    """

    vec_a = np.asarray(vec_a)[None, :]
    vec_b = np.asarray(vec_b)[None, :]

    if axis >= 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 >= len_b
    tmp = np.where(mask, len_a, len_b)
    np.putmask(len_b, ~mask, len_a)
    len_a = tmp

    mask = len_c > 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 > 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]

编辑：当使用较大的浮点数执行计算时，该问题肯定与此相关并消失：float64

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