x86中的三角函数指令错误是什么？

Question

我在哪里可以找到关于x86处理器上三角函数指令的误差范围信息，例如 fsincos？

- mdakin

---

我期望它符合IEEE 754的要求，是1 ulp。 - lhf

IEEE-754对三角函数并没有任何要求（即使有要求，要求也不会是1 ulp）。由IEEE-754标准化的运算通常需要被正确舍入，这大致相当于0.5 ulp的容差。 - Stephen Canon

相关链接：https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/ - 英特尔在之前的fsin文档中低估了1.3千万亿的误差边界。 - Peter Cordes

Answer 1

您所提问的问题很少是有趣的，并且最可能您真正想知道的是不同的东西。因此，让我先回答不同的问题：

“如何计算三角函数以达到一定的精度？”只需使用更长的数据类型。对于x86，如果您需要双倍精度的结果，请进行80位扩展双倍运算，这样就安全了。

“如何获得平台无关的精度？” 您需要一个专门的软件解决方案，例如 MPFR。

尽管如此，让我回到您最初的问题。简短的答案：对于小操作数，通常应该在1 ulp之内。 对于较大的操作数，情况会变得更糟。 唯一确定的方法是自己进行测试，就像这个人所做的那样。 处理器供应商没有可靠的信息。

- pentadecagon

---

谢谢。实际上我对超越指令的错误很感兴趣，但我也对如何获得更高的精度很感兴趣，你也回答了这个问题。 - mdakin

Answer 2

对于英特尔处理器，内置的超越函数指令的精度记录在Intel® 64和IA-32体系结构软件开发人员手册，第1卷第8.3.10节“超越指令精度”中：

Pentium处理器和后来的IA-32处理器，超越函数的最坏情况误差小于1 ulp（舍入到最近的（偶数）），在其他模式下小于1.5 ulps。

需要注意的是，1 ulp的误差界限适用于80位扩展精度格式，因为所有超越函数指令都提供扩展精度结果。关于三角函数指令FSIN、FCOS、FSCINCOS、FPTAN相对于数学参考的精度损失问题，由于使用66位机器PI进行参数约减，Intel已经确认。提供以下指导：

无论目标精度（单精度、双精度或双扩展精度），对于FSIN，将参数减小到绝对值大约为3π/4以下是安全的，对于FCOS、FSINCOS和FPTAN，将参数减小到绝对值小于约3π/8以下是安全的。例如，准确度测量表明，当|x|＜2.82时，FSIN的双扩展精度结果不会出现大于0.72ulp的误差[...]同样，当|x|＜1.31时，FCOS的双扩展精度结果不会出现大于0.82 ulp的误差[...]

此外，需要注意的是，对数函数指令FYL2X和FYL2XP1的1ulp误差界限仅在y=1时成立（这在英特尔早期的一些文档中并不清楚）:

FYL2X和FYL2XP1指令是两个操作数指令，并且仅在y等于1时保证在1 ulp范围内。当y不等于1时，最大ulp误差始终在1.35以内。

使用多精度库，可以轻松地对英特尔的声明进行测试。为了收集以下数据，我使用Richard Brent的MP库作为参考，并在指定的间隔内运行了2³¹个随机测试用例：

Intel Xeon CPU E3-1270 v2 "IvyBridge", Intel64 Family 6 Model 58 Stepping 9, GenuineIntel

2xm1 [-1,1]        max. ulp = 0.898306 at x = -1.8920e-001 (BFFC C1BED062 C071D472)
sin [-2.82,+2.82]  max. ulp = 0.706783 at x =  5.1323e-001 (3FFE 8362D6B1 FC93DFA0)
cos [-1.41,+1.41]  max. ulp = 0.821634 at x = -1.3201e+000 (BFFF A8F8486E 591A59D7)
tan [-1.41,+1.41]  max. ulp = 0.990388 at x =  1.3179e+000 (3FFF A8B0CAB9 0039C790)
atan [-1,1]        max. ulp = 0.747328 at x =  1.2252e-002 (3FF8 C8BB9E06 B9EB4DF8), y =  3.9204e-001 (3FFD C8B8DC94 AA6655B4)
y2lx [0.5,2.0]     max. ulp = 0.994396 at x =  1.0218e+000 (3FFF 82C95B56 8A70EB2D), y =  1.0000e+000 (3FFF 80000000 00000000)
yl2x [1.0,1.2]     max. ulp = 1.202769 at x =  1.0915e+000 (3FFF 8BB70F1B C5F7E103), y = -9.8934e-001 (BFFE FD453A23 AC926478)
yl2xp1 [-0.7,1.44] max. ulp = 0.990469 at x =  2.1709e-002 (3FF9 B1D61A98 BF349080), y =  1.0000e+000 (3FFF 80000000 00000000)
yl2xp1 [-1, 1]     max. ulp = 1.206979 at x =  9.1169e-002 (3FFB BAB69127 C1D5C158), y = -9.9281e-001 (BFFE FE28A91F 132F0C35)

虽然这种非穷尽测试不能“证明”误差边界，但找到的最大误差似乎证实了英特尔的文档。

我没有现代的 AMD 处理器进行测试，但有一个旧的 32 位 Athlon CPU 的测试数据。完全披露：我设计了用于 32 位 Athlon 处理器中的超越函数指令的算法。我的精度目标是所有指令的小于 1 ulp；但是上面提到的三角函数的 66 位机器 PI 的参数减少也适用于此。

Athlon XP-2100 "Palomino", x86 Family 6 Model 6 Stepping 2, AuthenticAMD

2xm1 [-1,1]        max. ulp = 0.720006 at x =  5.6271e-001 (3FFE 900D9E90 A533535D)
sin [-2.82, +2.82] max. ulp = 0.663069 at x = -2.8200e+000 (C000 B47A7BB2 305631FE)
cos [-1.41, +1.41] max. ulp = 0.671089 at x = -1.3189e+000 (BFFF A8D0CF9E DC0BCA43)
tan [-1.41, +1.41] max. ulp = 0.783821 at x = -1.3225e+000 (BFFF A947067E E3F4C39C)
atan [-1,1]        max. ulp = 0.665893 at x =  5.5333e-001 (3FFE 8DA6B606 C58B206A) y =  5.5169e-001 (3FFE 8D3B9DC8 5EA87546)
yl2x [0.4,2.5]     max. ulp = 0.716276 at x =  6.9826e-001 (3FFE B2C128C3 0EF1EC00) y = -1.2062e-001 (BFFB F7064049 BC362838)
yl2xp1 [-1,4]      max. ulp = 0.691403 at x =  1.9090e-001 (3FFC C37C0397 F8184934) y = -2.4796e-001 (BFFC FDE93CA9 980BF78C)

AMD64架构程序员手册第1卷中的6.4.5.1节“超越函数结果的精度”文档记录了误差范围如下：

x87计算以双扩展精度格式进行，因此超越函数为每种浮点数据类型提供的结果精确到最后一位（ulp）的误差不超过1个单位。

- njuffa

Answer 3

你可以阅读Intel® 64 和 IA-32 架构开发人员手册：卷 1第8.3.10节关于超越指令精度的内容。其中有一个精确的公式，但也有更易理解的陈述：

使用 Pentium 处理器及其后续的 IA-32 处理器，当四舍五入到最近的偶数时，超越函数的最坏情况误差小于 1 ulp，而在其他模式下舍入时小于 1.5 ulps。

- MvG

---

在考虑三角函数时，特别需要牢记的是，它们的精度边界是根据使用66位近似值pi的参考函数计算的（请参见同一文档中的8.3.8）。如果将结果与数学上精确的函数进行比较（大多数人可能会天真地想要这样做），误差可能会比1 ulp大得多（一旦超出函数的基本域，误差就会非常快地增长）。 - Stephen Canon

自从Bruce Dawson指出范围缩减导致fsin输入接近+Pi时发生灾难性取消后，英特尔已经更正了该文档：英特尔低估了1.3千万亿的误差界限。 - Peter Cordes