c++ - float 不准确导致计算错误

Question

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

1E+81E+121E+16

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

问题(已更新 - 具体化):因为 float 的值 float (和 double )b^2 & 4ac 有误差，足够小(相对于这些非常大的数字 [测量的绝对误差约为 1E+18])，但作为 D = = 它们之间的差异，所以当D(从较大的值侧)到所提到的顺序值时，不准确性是(1E+18) ，它的值开始在+1E+18 .. -1E+18范围内波动(即波动范围大于实际值的[-100% .. +100%]!

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

注意:当我说“当 D 接近于零时”实际上 D 离零还很远，所以我不能直接将它分配给 zerro in这个范围的值。

我考虑过使用定点计算(这可以使我摆脱问题)。但是，不建议在报价逻辑中使用(因为它们的优化程度要低得多，而且可能会非常慢)。

我的问题:我该如何尝试解决我的问题？可能有一些针对我的情况的通用解决方案？非常感谢任何建议!

PS:所有公式都很好(我在 excel 中计算了所有公式并得到了正确的结果，当我的代码中的 float 失败时)。

PPS:我尝试用 double 代替 float (不是所有计算，但我的 a、b 和 c 现在是 double ) & 问题并没有消失。

更新:我犯了一个错误 - 混淆了 a、b 和 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

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

Answer 1

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

诀窍是在这种情况下使用替代公式，请参见此处: 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% 的错误。

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