浮点数不准确导致计算错误

Question

我的代码求解二次方程（在游戏逻辑刻度中）以解决任务 - 沿着空间中可移动物体的轨道找到卫星刻度偏移。 而且我在判别式（更远的D）计算中遇到了错误。我会提醒：D = b^2 - 4ac。 因为它是大物体的轨道，所以我的a,b & c 是这样的序数：

1E+8 1E+12 1E+16

据此，b^2 是关于1E+24 的订单数，而4ac 也是关于1E+24 的订单数。 但是这个方程根的数字要少得多，因为它们只是场景中的坐标。所以根大约是1E+3 ... 1E+4。

问题（已更新-具体化）：因为浮点数（和双精度）b^2 和4ac 的值的浮动有不准确，这是足够小的（相对于这些非常大的数字 [测量的绝对不准确度大约有 1E+18]），但是因为 D == 它们之间的差异，所以当 D 是（从更大的值方面) 到上面提到的不准确的订单值是 (1E+18)，它的值开始在 +1E+18 .. -1E+18 左右的范围内波动（即波动范围大于实际值的 [-100% .. +100%]！

显然，这种波动会导致错误的（甚至是错误的方向）刻度偏移。我的卫星开始摇晃（这很糟糕））。

注意：当我说“当D 接近零时”实际上D 离零还很远，所以我不能在这个值范围内将它分配给零。

我考虑过使用定点计算（这可以让我摆脱问题）。但是，不建议在滴答逻辑中使用（因为它们的优化程度要低得多，而且可能会很慢）。

我的问题：如何尝试解决我的问题？我的情况可能有一些常见的解决方案吗？非常感谢您的建议！

PS：所有公式都很好（当我的代码中的浮点数失败时，我在 excel 中计算了所有公式并得到了正确的结果）。

PPS：我尝试了双精度浮点数（不是所有计算，但我的 a、b 和 c 现在是双精度数）并且问题没有消失。

更新：我犯了一个错误-混淆了a、b和c的顺序。所以“b^2 是关于1E+16 的订单数，&4ac 是关于1E+28”是错误的。现在它都固定为1E+24。 （我已经写到这个已经写的cmets可以理解了）

更新#2：“问题”部分已具体化。

更新#3：值的真实案例（供参考）： 注意：这里作为“准确值”，我将在 Excel 中手动计算的值标记出来。

a == 1.43963872E+8
b == 3.24884062357827E+12
c == 1.83291898112689E+16

//floats:
b^2 == 1.05549641E+25
4ac == 1.05549641E+25
D == 0.0
root:
y = -1.12835273E+4

//doubles:
b^2 == 1.0554965397412443E+25
4ac == 1.0554964543412880E+25
D == 8.5399956328598733E+17
roots:
y1 == -1.1280317962726038E+4
y2 == -1.1286737079932651E+4

//accurate values:
b^2 == 1.05549653974124E+25
4ac == 1.05549645434129E+25
D == 8.53999563285987E+17
roots: 
y1 == -1.128031796E+4 
y2 == -1.128673708E+4

双打看起来不错，但不是，因为这里我只给出了部分计算 - 这里我从相同的 a、b 和 c 值开始，但它们在我的代码中的实际值也被计算出来。并且包含不准确性，即使使用双打也会产生问题。

Answer 1

使用标准的二次公式会产生“灾难性的抵消”，其中两个相同大小的数字相减会导致精度损失。

诀窍是在这种情况下使用替代公式，请参见此处： https://math.stackexchange.com/a/311397

更新：我误读了您的问题。我认为问题更可能是您的结果对输入数字的敏感性。让我们选择，说

a = 4e8
b = -1e12
c = 6.2e14

解决方案是 ~1138 和 1361。现在如果你计算相对导数。我可以在 Julia 中通过使用 ForwardDiff.jl 包的自动微分来做到这一点：

julia> import ForwardDiff.Dual

julia> function p(a,b,c)
    D = sqrt(b^2-4*a*c)
    (-b+D)/(2a), (-b-D)/(2a)
end

julia> p(a,Dual(b,b),c)
(Dual(1361.803398874989,15225.424859373757),Dual(1138.196601125011,-12725.424859373757))

julia> p(Dual(a,a),b,c)
(Dual(1361.803398874989,-8293.614129124373),Dual(1138.196601125011,5793.614129124373))

julia> p(a,b,Dual(c,c))
(Dual(1361.803398874989,-6931.8107302493845),Dual(1138.196601125011,6931.8107302493845))

这里的结果是两个解和它们的缩放导数（即 (df/dx)*x）。请注意，它们都是 O(10000) 的量级，所以如果输入有 0.000001% 的误差，输出也会有 0.1% 的误差。

这里唯一的解决方案是重新表述您的问题，使其对输入值不那么敏感。

Answer 2

查看我对这个问题的回答：Quadratic equation in Ada

诀窍是始终使用

x1 = (-b - sign(b) * sqrt(b^2 - 4ac)) / 2a

作为第一个根，并使用

x1 * x2 = c / a

找到第二个。这样，您就可以避免 4ac

如果您声称的问题是 b^2 和 4ac 具有相同的幅度，那么与 b 相比，delta 实际上很小，并且您没有舍入问题，您也许应该重新调整您的问题（两种解决方案都非常接近 -b/ 2a)。

Answer 3

C++ 有一个标准数学库函数fma()，它提供了一种简单的方法，通过对判别式 d = √ 的稳健计算，在给定的浮点类型内尽可能准确地计算二次方程的根(b² - 4ac):

/*
  Compute a*b-c*d with error < 1.5 ulp

  Claude-Pierre Jeannerod, Nicolas Louvet, and Jean-Michel Muller, 
  "Further Analysis of Kahan's Algorithm for the Accurate Computation of 2x2 Determinants". 
  Mathematics of Computation, Vol. 82, No. 284, Oct. 2013, pp. 2245-2264
*/
T diff_of_products (T a, T b, T c, T d)
{
    T w = d * c;
    T e = fma (-d, c, w);
    T f = fma (a, b, -w);
    return f + e;
}

/* George E. Forsythe, "How Do You Solve a Quadratic Equation"
   Stanford University Technical Report No. CS40 (June 16, 1966)
*/ 
T a, b, c;
T d = diff_of_products (b, b, 2*a, 2*c);
T x1 = 2*c / (-b - sqrt (d));
T x2 = 2*c / (-b + sqrt (d));

fma() 实现的融合乘加运算 (FMA) 映射到大多数现代处理器架构上的单个硬件指令。由于 FMA 在加法之前计算完整的、未舍入的双倍宽度乘积，因此它用于准确计算乘积的误差。

正如 Simon Byrne 在 his answer 中提到的那样，手头的具体问题是病态的，准确的计算无法解决这个问题，只有重新制定基础数学可以。