为什么(1-x)(1 + x)更喜欢(1-x ^ 2)?

我正在查看arcsin的运行时库实现,它是通过计算来实现的:

 ArcTan(X, Sqrt(1 - X*X)) 

然而,计算1 - X*X的代码实际上评估了(1-X)*(1+X) 。 有更好的理由select后者吗? 我怀疑后者减less了X接近于零的舍入误差,但是我不能解释为什么会如此。

ArcTan(X, Sqrt(1-X*X))相对于X的导数为1/Sqrt(1-X*X) 。 这是无限的| X | 因此,当X接近1或-1时,任何评估错误都会对结果产生巨大的影响。 因此,在这些情况下评估最小化错误是至关重要的。

当X接近1时, 1-X评估没有错误(在IEEE 754或任何良好的浮点系统中,因为结果的规模使得其最低有效位至less与最低有效位一样低1或X,所以确切的math结果在可用的有效位之外没有位)。 由于1-X是精确的,通过考虑ArcTan(X, Sqrt((1-X)*(1+X+e))对e的导数,考虑1+X误差的影响,其中e是当x接近1且e很小时,导数为-1/10(以Maple为导数,用1代替x,得到-1/(sqrt(4+2e)*(5+2e) ,然后用0代替,得到-1/10。)因此, 1+X的误差并不重要。

因此,评估expression式为ArcTan(X, Sqrt((1-X)*(1+X))是评估它的好方法。

这种情况在X接近-1时是对称的。 ( 1+X没有错误, 1-X不重要。)

相反,如果我们考虑X*X的误差,那么ArcTan(X, Sqrt(1-X*X+e))对e的导数是当X接近1时,约为-1 /(2 * sqrt (e)*(1 + e)),所以当e很小的时候是很大的。 所以当X接近1时,评估X*X一个小错误会导致结果中的大错误。


问Pascal Cuoq指出,当评估一个函数f(x)时,我们通常希望最小化最终结果中的相对误差。 而且,正如我所指出的那样,计算过程中出现的错误通常是由于浮点四舍五入导致的中间结果的相对误差。 我可以在上面忽略这个,因为当X接近1时我认为是函数,所以所考虑的中间值(1 + X和X * X)和最终值都有大约1的值,所以将值这些大小不会有什么重大改变。

不过,为了完整起见,我更仔细地检查了这种情况。 在Maple中,我写了g := arctan(x, sqrt((1-x*x*(1+e0))*(1+e1))*(1+e2)) ,因此允许相对误差e0,e1 ,和e2分别计算x*x1-x*xsqrt ,我写成h:= arctan(x, sqrt((1-x)*(1+x)*(1+e0))*(1+e2)) 。 注意,在这种情况下,e0将1-x1+x的三个误差和它们的相乘结合起来; 完全误差项可以是(1+ea)*(1+eb)*(1+ec) ,但是这对于e0来说有效地为1+e0并且具有更大的可能范围。

然后,我检查了这些函数的导数(e0,e1和e2除以abs(f(x)),其中f是理想函数arctan(x, sqrt(1-x*x)) 。 例如,在Maple中,我检查了diff(g, e0) / abs(f(x)) 。 我没有对这些进行全面的分析评估。 我检查了一些x值接近0和​​接近1的值,以及e0,e1和e2的值在它们的一个极限值-2-54

对于0附近的x,这些值都是大约1或更小。 也就是说,计算中的任何相对误差都会导致结果中相似的相对误差,或者更小。

对于接近1的x,e1和e2的导数值很小,大约在10 -8或更小。 然而,这两种方法的e0的导数值是非常不同的。 对于1-x*x方法,值约为2•10 7 (使用x = 1-2 -53 )。 对于(1-x)*(1+x)方法,其值约为5•10 -9

总之,两种方法在x = 0附近差别不大,但后一种方法在x = 1附近明显更好。

我写了下面的程序,以获得单精度的一些实证结果。

 #include <float.h> #include <math.h> #include <stdio.h> long double d1, d2, rel1, rel2; float i1, i2; int main() { float f; for (f = nextafterf(0, 2); f <= 1; f = nextafterf(f, 2)) { long double o = 1.0L - ((long double)f * f); float r1 = (1 - f) * (1 + f); float r2 = 1 - f * f; long double c1 = fabsl(o - r1); long double c2 = fabsl(o - r2); if (c1 > d1) d1 = c1; if (c2 > d2) d2 = c2; if (c1 / o > rel1) rel1 = c1 / o, i1 = f; if (c2 / o > rel2) rel2 = c2 / o, i2 = f; } printf("(1-x)(1+x) abs:%Le relative:%Le\n", d1, rel1); printf("1-x*x abs:%Le relative:%Le\n\n", d2, rel2); printf("input1: %a 1-x:%a 1+x:%a (1-x)(1+x):%ao:%a\n", i1, 1-i1, 1+i1, (1-i1)*(1+i1), (double)(1 - ((long double)i1 * i1))); printf("input2: %ax*x:%a 1-x*x:%ao:%a\n", i2, i2*i2, 1 - i2*i2, (double)(1 - ((long double)i2 * i2))); } 

几句话:

  • 我使用了80位long double来计算元数据。 这还不足以准确地表示x所有值的错误,但是我担心程序会变得过于慢而精度更高。
  • 参考值o被计算为1.0L - ((long double)f * f) 。 这总是最接近真正结果的long double数,因为(long double)f * f是精确的(看来,forms1 - x*x有时可能会更好:))。

我得到了以下结果:

 (1-x)(1+x) abs:8.940394e-08 relative:9.447410e-08 1-x*x abs:4.470348e-08 relative:8.631498e-05 input1: 0x1.6a046ep-1 1-x:0x1.2bf724p-2 1+x:0x1.b50238p+0 (1-x)(1+x):0x1.0007bep-1 o:0x1.0007bc6a305ep-1 input2: 0x1.ffe96p-1 x*x:0x1.ffd2cp-1 1-x*x:0x1.6ap-12 o:0x1.69f8007p-12 

根据这些结果, 1 - x*x具有更好的绝对精度, (1-x)*(1+x)具有更好的相对精度。 浮点都是关于相对精度(整个系统的devise是为了能够相对准确地表示小数值和大数值),所以后一种forms是优选的。

编辑:计算最终的错误更有意义,如埃里克的答案所示。 如ArcTan(X, Sqrt(1 - X*X))这样的expression式中的子expression式可能被选中的原因并不在于其总体上更精确,而是因为它在最重要的地方是准确的。 将下面的行添加到循环的主体中:

  long double a = atan2l(f, sqrtl(o)); float a1 = atan2f(f, sqrtf(r1)); float a2 = atan2f(f, sqrtf(r2)); long double e1 = fabsl(a - a1); long double e2 = fabsl(a - a2); if (e1 / a > ae1) ae1 = e1 / a, i1 = f; if (e2 / a > ae2) ae2 = e2 / a, i2 = f; 

使用atan2l(f, sqrtf(r1))可能会atan2l(f, sqrtf(r1))意义,因为我没有与系统完全相同的ArcTan函数。 无论如何,有了这些注意事项,对于完整的expression式,(1-x)(1 + x)版本的区间[-1 … 1]的最大相对误差为1.4e-07,对于1的版本为5.5e-7 -x 2版本。