如何：x86 中的 pow(real, real)

Posted 2023-02-17

【中文标题】如何：x86 中的 pow(real, real)

【英文标题】：How to: pow(real, real) in x86

【发布时间】：2011-01-09 09:18:11

【问题描述】：

我正在寻找在 x86 汇编中实现 pow(real, real)。我也想了解算法是如何工作的。

【问题讨论】：

glibc 对pow() 函数的实现是in sysdeps/ieee754/dbl-64/e_pow.c。它使用一些 FP 位模式的整数检查，以及一些 FP 乘法和加法，但不使用任何特殊的 x87 指令。对于 x86-64，它被编译成 __ieee754_pow_sse2() (by this code that #includes it)。无论如何，x87 并不是在现代 CPU 上执行此操作的最佳方式。

我假设 glibc 的代码比 x87 更准确或更快。可能两者都有，但可能更准确（正确四舍五入到最接近的值）。但是，它不使用循环，并且单步执行指令，pow(1.175, 33.75) 并没有那么多。 FYL2X 在现代 CPU 上是一条非常慢的指令（约 100 个周期），因此击败它应该不难。

相关：Optimizations for pow() with const non-integer exponent? 有一个快速近似版本（使用 SIMD 内在函数）。有关提供 pow 的 SIMD 数学库，另请参阅 Where is Clang's '_mm256_pow_ps' intrinsic?。

【参考方案1】：

只需计算为2^(y*log2(x))。

有一个 x86 指令 FYL2X 来计算 y*log2(x) 和一个 x86 指令 F2XM1 来做幂运算。 F2XM1 需要 [-1,1] 范围内的参数，因此您必须在两者之间添加一些代码来提取整数部分和余数，对余数取幂，使用 FSCALE 以适当的 2 次方缩放结果。

【讨论】：

我知道这是一个旧线程，但这里有一个实现：madwizard.org

正好有这里代码：fld power 987654324 fyl2x 987654326 fld st(1) 987654328 f2xm1 987654330 fscale 987654332 fstp st 987654334 @

【参考方案2】：

好的，我按照你的建议在 x86 中实现了power(double a, double b, double * result);。

代码：http://pastebin.com/VWfE9CZT

%define a               QWORD [ebp+8]
%define b               QWORD [ebp+16]
%define result          DWORD [ebp+24]
%define ctrlWord            WORD [ebp-2]
%define tmp             DWORD [ebp-6]

segment .text
    global power

power:
    push ebp
    mov ebp, esp
    sub esp, 6
    push ebx

    fstcw ctrlWord
    or ctrlWord, 110000000000b
    fldcw ctrlWord

    fld b
    fld a
    fyl2x

    fist tmp

    fild tmp
    fsub
    f2xm1
    fld1
    fadd
    fild tmp
    fxch
    fscale

    mov ebx, result
    fst QWORD [ebx]

    pop ebx
    mov esp, ebp
    pop ebp
    ret

【讨论】：

我可以建议您继续并在您的答案中包含该代码吗？

您应该 sub esp, 8 保持对齐以便推送 ebx。您还可以交换 tmp 和 ControlWord，例如%define tmp DWORD [ebp-4]，所以它是对齐的。

返回 double 而不是输出 arg 会更有意义，因此您可以将值保留在 st0 中。或者，如果您坚持使用指针，请将指针加载到 EAX、ECX 或 EDX 中，这样您根本不必保存/恢复 EBX。此外，您应该在完成后恢复原始舍入模式。 （例如，将原件保存在寄存器中，然后存储并fldcw 它）。这使它设置为截断（接近零），而不是默认的舍入到最近。 efg2.com/Lab/Library/Delphi/MathFunctions/FPUControlWord.Txt.

更新：frndint 在很多 CPU 上都很慢 (agner.org/optimize)，所以 fist/fild 实际上更好。如果要使用截断进行转换，可以使用SSE3 fisttp（如果可用），而不是更改 FP 舍入模式并恢复它。或者只是使用 SSE2 cvttsd2si 进行截断转换。

【参考方案3】：

这是我使用“The Svin”的主要算法的函数。我将它包装在 __fastcall 和 __declspec(naked) 装饰中，并添加了代码以确保 base/x 为正数。如果 x 为负，则 FPU 将完全失败。您需要检查“x”符号位，并考虑“y”的奇数/偶数位，并在完成后应用符号！让我知道你对任何随机读者的想法。如果可能的话，使用 x87 FPU 代码寻找更好的版本。它与 Microsoft VC++ 2005 一起编译/工作，我出于各种原因一直坚持使用它。

与 ANSI pow(x,y) 的兼容性：非常好！更快，可预测的结果，负值被处理，无效输入没有错误反馈。但是，如果您知道 'y' 始终可以是 INT/LONG，请不要使用此版本；我发布了 Agner Fog 的版本，并进行了一些调整，以避免非常慢的 FSCALE，请搜索我的个人资料！他是那些有限情况下最快的 x87/FPU 方式！

extern double __fastcall fs_Power(double x, double y);

// Main Source: The Svin
// pow(x,y) is equivalent to exp(y * ln(x))
// Version: 1.00

__declspec(naked) double __fastcall fs_Power(double x, double y)  __asm 
    LEA   EAX, [ESP+12]         ;// Save 'y' index in EAX
    FLD   QWORD PTR [EAX]       ;// Load 'y' (exponent) (works positive OR negative!)
    FIST  DWORD PTR [EAX]       ;// Round 'y' back to INT form to test for odd/even bit
    MOVZX EAX, WORD PTR [EAX-1] ;// Get x's left sign bit AND y's right odd/even bit!
    FLD   QWORD PTR [ESP+4]     ;// Load 'x' (base) (make positive next!)
    FABS            ;// 'x' MUST be positive, BUT check sign/odd bits pre-exit!
    AND   AX, 0180h ;// AND off all bits except right 'y' odd bit AND left 'x' sign bit!
    FYL2X       ;// 'y' * log2 'x' - (ST(0) = ST(1) * log2 ST(0)), pop
    FLD1        ;// Load 1.0f: 2 uses, mantissa extract, add 1.0 back post-F2XM1
    FLD   ST(1) ;// Duplicate current result
    FPREM1      ;// Extract mantissa via partial ST0/ST1 remainder with 80387+ IEEE cmd
    F2XM1       ;// Compute (2 ^ ST(0) - 1)
    FADDP ST(1), ST ;// ADD 1.0f back! We want (2 ^ X), NOT (2 ^ X - 1)!
    FSCALE      ;// ST(0) = ST(0) * 2 ^ ST(1) (Scale by factor of 2)
    FFREE ST(1) ;// Maintain FPU stack balance
;// Final task, make result negative if needed!
    CMP   AX, 0180h    ;// Combo-test: Is 'y' odd bit AND 'x' sign bit set?
    JNE   EXIT_RETURN  ;// If positive, exit; if not, add '-' sign!
        FCHS           ;// 'x' is negative, 'y' is ~odd, final result = negative! :)
EXIT_RETURN:
;// For __fastcall/__declspec(naked), gotta clean stack here (2 x 8-byte doubles)!
    RET   16     ;// Return & pop 16 bytes off stack

好的，为了结束这个实验，我使用 RDTSC CPU 时间戳/时钟计数器指令进行了基准测试。我遵循了使用“SetPriorityClass(GetCurrentProcess(), HIGH_PRIORITY_CLASS);”将进程设置为高优先级的建议。我关闭了所有其他应用程序。

结果：我们的复古 x87 FPU 数学函数“fs_Power(x,y)”比 MSCRT2005 pow(x,y) 版本快 50-60%，后者使用标记为“_pow_pentium4:”的相当长的 SSE 代码分支，如果它检测到 64 位 >Pentium4+ CPU。所以啊啊啊！！ :-)

注意：(1) CRT pow() 有一个大约 33 微秒的初始化分支，它在这个测试中显示了 46,000 个。在 1200 到 3000 次循环之后，它以正常平均值运行。我们手工制作的 x87 FPU 美感运行一致，第一次调用时没有初始化惩罚！

(2) 虽然 CRT pow() 输掉了所有测试，但它确实在一个区域中获胜：如果您输入了狂野的、巨大的、超出范围/溢出的值，它很快就会返回错误。由于大多数应用程序在典型/正常使用时不需要错误检查，因此无关紧要。

https://i.postimg.cc/QNbB7ZVz/FPUv-SSEMath-Power-Proc-Test.png

第二次测试（我不得不再次运行它以在图像捕捉后复制/粘贴文本）：

 x86 fs_Power(2, 32): CPU Cycles (RDTSC): 1248
MSCRT SSE pow(2, 32): CPU Cycles (RDTSC): 50112

 x86 fs_Power(-5, 256): CPU Cycles (RDTSC): 1120
MSCRT SSE pow(-5, 256): CPU Cycles (RDTSC): 2560

 x86 fs_Power(-35, 24): CPU Cycles (RDTSC): 1120
MSCRT SSE pow(-35, 24): CPU Cycles (RDTSC): 2528

 x86 fs_Power(64, -9): CPU Cycles (RDTSC): 1120
MSCRT SSE pow(64, -9): CPU Cycles (RDTSC): 1280

 x86 fs_Power(-45.5, 7): CPU Cycles (RDTSC): 1312
MSCRT SSE pow(-45.5, 7): CPU Cycles (RDTSC): 1632

 x86 fs_Power(72, -16): CPU Cycles (RDTSC): 1120
MSCRT SSE pow(72, -16): CPU Cycles (RDTSC): 1632

 x86 fs_Power(7, 127): CPU Cycles (RDTSC): 1056
MSCRT SSE pow(7, 127): CPU Cycles (RDTSC): 2016

 x86 fs_Power(6, 38): CPU Cycles (RDTSC): 1024
MSCRT SSE pow(6, 38): CPU Cycles (RDTSC): 2048

 x86 fs_Power(9, 200): CPU Cycles (RDTSC): 1152
MSCRT SSE pow(9, 200): CPU Cycles (RDTSC): 7168

 x86 fs_Power(3, 100): CPU Cycles (RDTSC): 1984
MSCRT SSE pow(3, 100): CPU Cycles (RDTSC): 2784

任何现实世界的应用程序？是的！ Pow(x,y) 被大量用于帮助将 CD 的 WAVE 格式编码/解码为 OGG，反之亦然！当您对整整 60 分钟的 WAVE 数据进行编码时，节省时间的回报将非常显着！ OGG/libvorbis 中使用了许多数学函数，如 acos()、cos()、sin()、atan()、sqrt()、ldexp()（非常重要）等。所以像这样的微调版本，不用麻烦/不需要错误检查，可以节省大量时间！！

我的实验是为 NSIS 安装程序系统构建 OGG 解码器的结果，这导致我将算法所需的所有数学“C”库函数替换为您在上面看到的内容。好吧，几乎，我需要 x86 中的 acos()，但我仍然找不到任何东西......

问候，并希望这对任何喜欢修补的人有用！

【讨论】：

为什么不直接使用fabs 指令呢？它要快得多（现代 AMD 和 Intel 上的 1 个周期延迟），并且不会在您对整个事物进行 qword 加载之前通过写入 double 的一半来导致存储转发停止。 （您仍然可以在最后对内存中的原始值进行分支）。此外，您不需要保存/恢复 EDX：它在所有标准调用约定中都被调用破坏了。

你也可以避免fxch st1。弹出堆栈时使用fstp st(1) 保持st0 = st0。

fprem 很慢。如果您只是使用它来获取整数和小数部分，请使用 frndint 并减去。 （嗯，根据agner.org/optimize frndint 在 Intel 上也很慢，但在 Ryzen 上很快。奇怪，因为 SSE/AVX roundpd 很快）

(1) 是的，关于 FRNDINT，我一直在研究 Agner Fog，他说它很慢。 (2) 老实说，在我将浮点数与 INT 命令加载到内存中并关闭符号位之后，还没有对 FABS 进行定时器测试。但是考虑到 Agner，我认为它可能会更快，而不是 100%。 ;) (3) 我保留了 EDX，因为我更愿意知道当我从呼叫返回时所有寄存器都被保留。我记得在 VC++ 中构建 DLL，在 VB/VBA 中调用它们，但如果没有恢复 EDI，你会得到一个糟糕的 DLL 调用约定。我更多地内联 ASM，尝试使用所有寄存器，所以我不会破坏任何东西。

是的，顶部的无条件fabs 是我要做的；直到 FP 值准备好很久之后，才可能出现错误预测的分支，因此 FP 延迟可以隐藏分支未命中惩罚。你不需要双字加载，只需要字节。 test [mem], imm8 / jz 小于 mov eax, [esp+4] / test eax,eax / jns。尽管在慢速微编码 FP 指令之前执行 mov eax, [esp+4] 可能有助于减少最后 test/jns 的分支错误预测检测延迟，以防微编码指令占用前端。