x87 FPU 计算 e 驱动 x，也许带有泰勒级数？

Question

我正在尝试在 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 获取总和和被加数的指数。如果 (exp(sum) - exp(summand)) > 23，被加数将不再影响单精度中的位表示（因为单精度中的尾数有 23 位）--> 所以退出。2.：将被加数与0比较，如果是0显然也不会影响结果了。

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

Answer 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缩放。

一个替代实现——基本上是由 fuz 建议的，让人想起 MSVC 为expC 标准库中的函数生成的代码——如下：

; 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，而不是使用fprem1 除后获得余数。

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

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/?

上面使用的符号是“操作？延迟？互惠吞吐量”。

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指令的一种可能方法——如果你愿意为了速度而牺牲代码大小——将是一种基于查找表的方法。如果你搜索，你可以找到几篇关于此的论文，尽管其中大多数都在付费墙后面. :-(这是一个可公开访问的参考。事实上，这可能是用 C 编写的优化数学库为实现exp. Ben Voigt 用 C# 编写的类似代码所做的. 与往常一样，基本问题就在这里，您是针对大小还是速度进行优化？换句话说，你是否经常调用这个函数并且需要它尽可能快？或者您是否不经常调用它并且只需要它是精确的，而不会在二进制文件中消耗太多空间，并且可能会通过将其从缓存中逐出而减慢热路径上的其他代码。）

但是fpremvs.frndint呢？

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

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

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

如果您不想费心对其进行基准测试，或者您想要在所有 CPU 上运行良好的通用东西，要么实现是好的。两人都不打算放火烧世界在性能方面，但正如我之前所说的，既会显著快于一个循环，试图计算泰勒级数。那里的开销会扼杀你的表现。