x87 FPU 计算 e 驱动 x，也许是泰勒级数？答案

【问题标题】：x87 FPU computing e powered x, maybe with a Taylor series?x87 FPU 计算 e 驱动 x，也许是泰勒级数？
【发布时间】：2017-12-10 22:43:21
【问题描述】：

我正在尝试在 x87 中计算函数 e^x（x 是单精度浮点数）。为了实现这一点，我使用泰勒级数（提醒一下：e^x = 1 + (x^1)/1! + (x^2)/2! + ... + (x^n)/n ！）

当我使用 x87 时，我可以以扩展精度（80 位而不是单精度 32 位）计算所有值。

到目前为止，我的认识是：我在 FPU 的两个单独的寄存器中拥有和和和（加上其他一些不太重要的东西）。我有一个循环，不断更新总和并将其添加到总和中。所以在模糊的伪代码中：

loop:
;update my summand
n = n + 1
summand = summand / n
summand = summand * x
;add the summand to the total sum
sum = sum + summand

我现在的问题是循环的退出条件。我想以这样一种方式设计它，即一旦将总和添加到总和不会影响单精度总和的值，它就会退出循环，即使我正在弄清楚实现这种退出条件的艰难方法非常复杂，会占用大量指令 -> 计算时间。

到目前为止，我唯一的好主意是： 1.：通过FXTRACT获取sum和summand的指数。如果 (exp(sum) - exp(summand)) > 23，则 summand 将不再影响单精度中的位表示（因为单精度中的尾数有 23 位）-> 所以退出。 2.：和0比较，如果是0，显然也不会影响结果了。

我的问题是，对于现有条件，是否有人会比我有更有效的想法？

【问题讨论】：

我知道这不是你的问题，但你为什么要使用泰勒级数？
现在了解破坏条件的难度...姑且称之为初学者错误吧，这是我第一个真正的 x87 项目 ;-)
如果你很懒，你可以用fyl2x 和f2xmi 来实现exp，不过它们是一些最慢的指令。我可能会缩小范围并使用 Padé 近似值。
但是，如果您致力于做泰勒的事情，您也可以将您的号码停止更改作为退出条件。您可以按照您的建议手动测试，但如果添加数字没有任何区别，您可以简单地进行添加，然后检测它没有任何区别。
@MaruanElMahgary 就是这样......感觉比现在更糟。 fstp 和 fld 是 1-µop 并且（相对）快，而像 fxtract 这样的特殊指令是微编码的并且速度很慢，所以如果这是替代方案，那么“奇怪的存储和重新加载”开始看起来相当不错。比所有这些更简单的替代方法就是没有动态退出条件。

标签： assembly x86 exponential exponentiation x87

【解决方案1】：

也许使用 x87 浮点单元计算 e^x 的最有效方法是以下指令序列：

; Computes e^x via the formula 2^(x * log2(e))
fldl2e                  ; st(0) = log2(e)        <-- load log2(e)
fmul    [x]             ; st(0) = x * log2(e)
fld1                    ; st(0) = 1              <-- load 1
                        ; st(1) = x * log2(e)
fld     st(1)           ; st(0) = x * log2(e)    <-- make copy of intermediate result
                        ; st(1) = 1
                        ; st(2) = x * log2(e)
fprem                   ; st(0) = partial remainder((x * log2(e)) / 1)  <-- call this "rem"
                        ; st(1) = 1
                        ; st(2) = x * log2(e)
f2xm1                   ; st(0) = 2^(rem) - 1
                        ; st(1) = 1
                        ; st(2) = x * log2(e)
faddp   st(1), st(0)    ; st(0) = 2^(rem) - 1 + 1 = 2^(rem)
                        ; st(1) = x * log2(e)
fscale                  ; st(0) = 2^(rem) * 2^(trunc(x * log2(e)))
                        ; st(1) = x * log2(e)
fstp    st(1)

结果留在st(0)。

如果您已经在浮点堆栈上输入了x，那么您可以稍微修改指令序列：

; st(0) == x
fldl2e
fmulp   st(1), st(0)
fld1
fld     st(1)
fprem
f2xm1
faddp   st(1), st(0)
fscale
fstp    st(1)

这将 x 乘以 log₂e，找到除以 1 的余数，将 2 取幂为该值的幂（f2xm1 也减去 1，然后我们再加 1），最后将其缩放 x × log₂e。

另一种实现——本质上是suggested by fuz，让人想起 MSVC 为 C 标准库中的 exp 函数生成的代码——如下：

; st(0) == x
fldl2e
fmulp    st(1), st(0)
fld      st(0)
frndint
fsub     st(1), st(0)
fxch                         ; pairable with FSUB on Pentium (P5)
f2xm1
fld1
faddp    st(1), st(0)
fscale
fstp     st(1)

主要区别在于使用frndint 和fsub 来获得-1.0 到+1.0 范围内的值，这是f2xm1 所要求的，而不是使用fprem 来获得余数除以 1 后。

为了了解这些指令的相对成本，我们将从Agner Fog's instruction tables 中提取数据。一种 ”？”表示对应的数据不可用。

Instruction      AMD K7      AMD K8      AMD K10    Bulldozer    Piledriver    Ryzen
-------------------------------------------------------------------------------------------
FLD/FLD1/FLDL2E  [all very fast, 1-cycle instructions, with a reciprocal throughput of 1]
FADD(P)/FSUB(P)  1/4/1       1/4/1       1/4/1      1/5-6/1      1/5-6/1       1/5/1
FMUL(P)          1/4/1       1/4/1       1/4/1      1/5-6/1      1/5-6/1       1/5/1
FPREM            1/7-10/8    1/7-10/8    1/?/7      1/19-62/?    1/17-60/?     2/?/12-50
FRNDINT          5/10/3      5/10/3      6/?/37     1/4/1        1/4/1         1/4/3
FSCALE           5/8/?       5/8/?       5/9/29     8/52/?       8/44/5        8/9/4
F2XM1            8/27/?      53/126/?    8/65/30?   10/64-71/?   10/60-73/?    10/50/?

上面使用的符号是“ops／latency／reciprocal throughput”。

Instruction      8087     287      387      486      P5         P6        PM      Nehalem
-------------------------------------------------------------------------------------------
FLD              17-22    17-22    14       4        1/0        1/?       1/1      1/1
FLD1             15-21    15-21    24       4        2/0        2/?       2/?      2/?
FLDL2E           15-21    15-21    40       8        5/2        2/?       2/?      2/?
FADD(P)/FSUB(P)  70-100   70-100   23-34    8-20     3/2        1/3       1/3      1/3
FMUL(P)          90-145   90-145   29-57    16       3/2        1/5       1/5      1/5
FPREM            15-190   15-190   74-155   70-138   16-64 (2)  23/?      26/37    25/14
FRNDINT          16-50    16-50    66-80    21-30    9-20 (0)   30/?      15/19    17/22
FSCALE           32-38    32-38    67-86    30-32    20-32 (5)  56/?      28/43    24/12
F2XM1            310-630  310-630  211-476  140-279  13-57 (2)  17-48/66  ~20/?    19/58

对于 P5 之前的版本，该值只是循环计数。对于 P5，符号是循环计数，括号中表示与其他 FP 指令重叠。对于 P6 及更高版本，表示法是 µops/latency。

很明显，f2xm1 是一条缓慢而昂贵的指令，但两种实现都使用它并且很难避免。尽管它很慢，但它仍然方式比将其实现为循环快。

（如果您愿意为了速度而牺牲代码大小，那么一种可能的方法是绕过缓慢的f2xm1 指令，那就是基于查找表的方法。如果您搜索，您可以找到几篇关于此发表的论文，尽管他们中的大多数都在付费墙后面。:-( Here's one publicly-accessible reference。事实上，这可能是用 C 编写的优化数学库为实现 exp 所做的事情。Ben Voigt wrote similar code in C#。与往常一样，基本问题在这里，您是否在优化大小或速度？换句话说，您是否经常调用此函数并且需要它尽可能快？或者您是否不经常调用它并且只需要精确，不会在二进制文件中占用太多空间，也不会通过从缓存中逐出热路径上的其他代码来减慢它的速度。）

但是fprem 与frndint 呢？

嗯，在 AMD CPU 上，fprem 的解码操作数比frndint 少，尽管在现代 AMD CPU 上frndint 比fprem 更快。再说一次，我想说这无关紧要，因为在那些 CPU 上，您将编写 SSE2/AVX 代码，而不是使用传统的 x87 FPU。（这些指令在后来的型号英特尔微架构上也会变慢，比如 Sandy Bridge、Haswell 等，但同样的论点也适用于那里，所以我在编译数据时完全忽略了这些。）真正重要的是这段代码如何执行较旧的 CPU，在较旧的 CPU 上，使用 fprem 的版本显然是赢家。

然而，在 Intel CPU 上，情况正好相反：frndint 通常比 fprem 快，但在 P6（Pentium Pro、Pentium II 和 Pentium III）上除外。

但是，我们只是在谈论frndint 与fprem——我们实际上是在谈论fprndint+fsub 与fprem。 fsub 是一个相对较快的指令，但是如果不执行它并对其计时，就很难预测该代码将具有什么性能。循环计数只能告诉我们这么多，而且fprem 序列总体上更短（指令更少，更重要的是，字节更少——18 对 22），这可能意味着显着的速度提升。

如果您不想费心对其进行基准测试，或者您想要在所有 CPU 上运行良好的通用程序，两种实现都很好。它们都不会在性能方面引起轰动，但正如我之前所说，它们都将明显比尝试计算泰勒级数的循环快。那里的开销会影响你的表现。

【讨论】：

我怀疑即使在 P5 上，仅使用简单的 x87 指令（如乘法和加法）的实现也可能是最快的。我不知道准确性；我知道 x87 三角函数不是很好，尤其是当范围缩小导致问题时，但是如果 e 的内部值优于 pi 的 x87 值 only accurate to 66 bits，则 IDK。我认为即使在 SSE 之前的 x87 天，优秀的数学库实现通常也使用基本操作完成所有工作。
如果这听起来很随意，那就是。我只是假设 e^x 有一些快速收敛的方法。我们可以用2^int(x) 做任何事情吗，我们可以通过将x 的整数部分填充到IEEE 二进制浮点数的指数中来获得？可能不是; ln(x) 和 log10(x) 的逆向适用，因为 ln(x) = log2(x) * ln(2)，那么您只需要在尾数域上对 log2(x) 进行多项式逼近。但是e^x / 2^x 是(e/2)^x，这不是一个常数。
你可能是正确的关于乘法和加法更快，@Peter。如果不进行测试，我实际上并不确定。但是与旧处理器上的循环和指令获取一样慢，我怀疑它在那里是否成立。也许在 P5 上，在 P6 上可能更是如此，FPU 最终（部分）流水线化了。然而，好的数学库实现过去只使用基本操作并不是真的。真正优秀的优化编译器，例如 Watcom，实际上确实生成了类似上述序列的内联复杂 x87 FPU 指令。事实上，我“最高效”的实现基本上是 Watcom 过去的......
…用于exp 功能，它是现代版本的 MSVC 在禁用 SSE 指令时继续使用的功能。我不确定这个谱系在 MSVC 上有多远，但我怀疑已经有一段时间了，因为最近没有关于 x87 代码生成的开发！我认为e的内部值是准确的；我没有听到任何其他建议，pi 的问题只是一个错误，简单明了。 trig insns 在 287 中是新的，可能只是仓促发货而没有经过全面测试。不过，做点什么来替换f2xm1 会很好。
...claims that, as of the Pentium III (and possibly earlier), a LUT-based approach would give no considerable speedup: "自从引入 Pentium III 模型以来，即使是基于查找表的 exp() 替换的最快实现，与优化的系统给出了这个函数的代码。”我希望我能在某个地方找到原文，但我不想仅仅出于好奇就订购这本书！