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

Inaccuracy of floating point numbers causes calculations error

本文关键字：计算 错误 不准确 浮点数      更新时间：2023-10-16

我的代码求解二次方程（在游戏逻辑tick中）以求解任务 - 在空间中沿着可移动对象的轨道找到卫星刻度偏移。而且我遇到了判别（D）计算中的错误。我会提醒：D = b^2 - 4ac。由于它是大对象的轨道，我的 a， b＆amp;c是订单数：

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

因此，b^2是关于1E+24的订单数量，＆amp;4ac也是1E+24。但是这个方程的根数要少得多，" cos它们只是场景上的坐标。因此根是1E+3 ... 1E+4。

问题（更新 - 具体化）：，因为浮子值（＆amp; doubles） b^2＆amp;4ac具有不准确性，足够小（相对于这些非常大的数字（测量的绝对不准确率）的顺序约为 1E+18]），但如D == 差异中的时，当D是从较大的值侧）到上述不准确的订单值为（ 1E+18），其值开始在约 +1E+18 .. -1E+18的范围内波动（即波动范围比[-100％.. 100％.. 100％]的实际价值！

>

显然，这种波动会导致错误（甚至是错误的定向）滴答偏移。我的卫星开始摆动（这很糟糕）。

注意：当我说"当D接近零时"时，D实际上距离零还足够远，所以我不能仅在此值范围内将其分配给Zerro。

我已经考虑使用定点计算（这可能使我摆脱问题）。但是，不建议在tick逻辑中使用（"它们的优化程度要少得多，并且可能非常慢）。

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

ps：所有公式都很好（我在excel＆amp;中都取得了适当的结果，当我的代码中的漂浮失败时）。

pps：我尝试了装有浮子的双打（并非全部计算，而是我的a，b＆amp; c现在是双打）＆amp;问题没有消失。

更新：我犯了一个错误 - 混乱的a订单顺序，b＆amp;c。因此，" b^2是关于1E+16的订单数，4ac是关于1E+28的订单数字是错误的。现在，它将其固定在1E+24上。（我已经写给已经写的评论是可以理解的）

更新＃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

双打看起来还不错，但不是，'cos在这里我只给出了一部分计算 - 在这里我从同一a，b＆amp;开始。C值，但还计算了它们在我的代码中的实际值。并包含不准确性，即使双打也会产生问题。

使用标准二次公式可以给出"灾难性的取消"，其中2个减去2个数量相同的幅度会损失精度。

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

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

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

对于〜1138和1361的解决方案。现在，如果计算相对导数。我可以使用forwarddiff.jl软件包在朱莉娅（Julia）进行自动差异：

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％。

这里唯一的解决方案是重新制定问题，以使其对输入值不太敏感。

请参阅我对这个问题的答案：ADA

中的二次方程式

诀窍是始终使用

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

作为第一个根，并使用

x1 * x2 = c / a

找到第二个。这样，您将自己对冲4ac＆lt;＆lt;B^2和-b SQRT（Delta）表现出灾难性的取消。

如果您所谓的问题是B^2和4ac具有相同的幅度，那么与B相比，Delta实际上很小，并且您没有圆形问题，您也许应该重新确定问题（这两种解决方案都非常接近-b/2a）。

c 具有标准的数学库函数 fma()，它提供了一种简单的方法，可以通过强大的计算在给定的浮点类型中在给定的浮点类型中尽可能准确地计算二次方程的根源判别d =√（b 2 -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()映射实现的融合乘数ADD操作（FMA），可在大多数现代处理器体系结构上提供单个硬件指令。当FMA在添加之前计算完整的，不可能的双宽度产品时，它用于准确计算产品的误差。

正如西蒙·伯恩（Simon Byrne）在他的回答中提到的那样，手头的具体问题是不明显的，准确的计算无法解决此问题，只有对基础数学的重新制定。

相关文章：