粒子物理学中不变质量的数值稳定计算？

Question

在粒子物理学中，我们必须计算很多invariant mass，这是针对二体衰变

当动量 (p1, p2) 有时与质量 (m1, m2) 相比非常大（高达 1000 倍或更多）时。在那种情况下，在计算机上用浮点数进行计算时，最后两项之间会发生很大的抵消。

对于任何输入，可以使用什么样的数字技巧来准确计算这个值？

问题是关于适当的数字技巧来提高浮点数计算的准确性，因此解决方案应该与语言无关。出于演示目的，首选 Python 中的实现。重新表述问题并增加基本运算量的解决方案是可以接受的，但建议使用其他数字类型（如小数或多精度浮点数）的解决方案则不可接受。

注意：原始问题以 Python 表达式的形式提出了一个简化的 1D 维度问题，但问题是针对动量以 3D 维度给出的一般情况。问题以这种方式重新表述。

Answer 1

如果你把例如表达式中的 m1 = 1e-4、m2 = 1e-4、p1 = 1 和 p2 = 1，双精度计算得到 4e-8，单精度计算得到 0.0。我假设，您的问题是关于如何通过单精度计算获得 4e-8。

您可以对上面的表达式进行泰勒展开（大约 m1 = 0 和 m2 = 0）。

e ~ e|(m1=0,m2=0) + de/dm1|(m1=0,m2=0) * m1 + de/dm2|(m1=0,m2=0) * m2 + ...

如果我计算正确，第零和一阶项为 0，二阶展开将为

e ~ (p1+p2)/p1 * m1**2 + (p1+p2)/p2 * m2**2

即使使用单精度计算，这也恰好产生 4e-8。如果需要，您当然可以在扩展中做更多的项，直到达到单个浮点数的精度限制。

编辑
如果 mi 并不总是比 pi 小很多，您可以进一步修改方程式以获得

$m_1^2 + m_2^2 + 2p_1p_2 \left[ \sqrt{\frac{m_1^2m_2^2}{p_1^2p_2^2}+\frac{m_1^2}{p_1^2}+\frac{m_2^2}{p_2^2}+1}-1\right]$

复杂的部分现在是方括号中的部分。对于范围广泛的 x 值，它本质上是 sqrt(x+1)-1。如果 x 很小，我们可以使用平方根的泰勒展开（例如 here）。如果 x 值较大，则公式工作得很好，因为 1 的加法和减法不再因浮点精度而丢弃 x 的值。因此 x 的一个阈值必须选择低于一个切换到泰勒展开。

Answer 2

通过 Whosebug 上列出的一些技巧和 Jakob Stark 在他的回答中描述的转换，可以将方程重写为不再遭受灾难性抵消的形式。

原题求的是一维的解法，有简单的解法，但实际中需要3维的公式，解法比较复杂。参见 this notebook for a full derivation。

Python 中 3D 数值稳定计算的示例实现：

import numpy as np

# numerically stable implementation
@np.vectorize
def msq2(px1, py1, pz1, px2, py2, pz2, m1, m2):
    p1_sq = px1 ** 2 + py1 ** 2 + pz1 ** 2
    p2_sq = px2 ** 2 + py2 ** 2 + pz2 ** 2
    m1_sq = m1 ** 2
    m2_sq = m2 ** 2
    x1 = m1_sq / p1_sq
    x2 = m2_sq / p2_sq
    x = x1 + x2 + x1 * x2
    a = angle(px1, py1, pz1, px2, py2, pz2)
    cos_a = np.cos(a)
    if cos_a >= 0:
        y1 = (x + np.sin(a) ** 2) / (np.sqrt(x + 1) + cos_a) 
    else:
        y1 = -cos_a + np.sqrt(x + 1) 
    y2 = 2 * np.sqrt(p1_sq * p2_sq)
    return m1_sq + m2_sq + y1 * y2

# numerically stable calculation of angle
def angle(x1, y1, z1, x2, y2, z2):
    # cross product
    cx = y1 * z2 - y2 * z1
    cy = x1 * z2 - x2 * z1
    cz = x1 * y2 - x2 * y1
    
    # norm of cross product
    c = np.sqrt(cx * cx + cy * cy + cz * cz)
    
    # dot product
    d = x1 * x2 + y1 * y2 + z1 * z2
    
    return np.arctan2(c, d)

数值稳定的实现永远不会产生负结果，这是天真的实现中经常出现的问题，即使是双精度。

让我们将数值稳定的函数与简单的实现进行比较。

# naive implementation
def msq1(px1, py1, pz1, px2, py2, pz2, m1, m2):
    p1_sq = px1 ** 2 + py1 ** 2 + pz1 ** 2
    p2_sq = px2 ** 2 + py2 ** 2 + pz2 ** 2
    m1_sq = m1 ** 2
    m2_sq = m2 ** 2
    
    # energies of particles 1 and 2
    e1 = np.sqrt(p1_sq + m1_sq)
    e2 = np.sqrt(p2_sq + m2_sq)

    # dangerous cancelation in third term
    return m1_sq + m2_sq + 2 * (e1 * e2 - (px1 * px2 + py1 * py2 + pz1 * pz2))

对于下图，momenta p1 和 p2 是从 1 到 1e5 中随机选取的，值 m1 和 m2 是从 1e-5 到 1e5 中随机选取的。所有实现都以单精度获取输入值。两种情况下的参考都是使用带有 100 位小数的朴素公式使用 mpmath 计算的。

朴素的实现会失去一些输入的所有准确性，而数值稳定的实现则不会。

粒子物理学中不变质量的数值稳定计算？

Numerically stable calculation of invariant mass in particle physics?

numeric