1

我的代码求解二次方程（在游戏逻辑刻度中）以解决任务 - 沿着空间中可移动物体的轨道找到卫星刻度偏移。而且我在判别（更远D）计算中遇到了错误。我会提醒：D = b^2 - 4ac。因为它是大物体的轨道，所以 mya和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==差异其中，当Dis（从较大的值一侧）到所提到的不准确的顺序值是1E+18（+1E+18 .. -1E+18价值！

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

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

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

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

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

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

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

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

4

3 回答 3

5

使用标准的二次公式可以给出“灾难性的抵消”，其中两个相同数量的数字相减会导致精度损失。

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

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

于 2016-12-12T08:01:05.350 回答
2

请参阅我对这个问题的回答：Ada 中的二次方程

诀窍是始终使用

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

作为第一个根，并使用

x1 * x2 = c / a

找到第二个。这样，您就可以避免 4ac << b^2 和 -b + sqrt(delta) 出现灾难性取消的情况。

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

于 2016-12-12T11:29:29.480 回答
1

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

/*
  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 在加法之前计算完整的、未舍入的双倍宽度乘积，因此它用于准确计算乘积的误差。

正如西蒙伯恩在他的回答中提到的那样，手头的具体问题是病态的，准确的计算无法解决这个问题，只有重新制定基础数学可以。

于 2016-12-12T17:49:15.083 回答