35

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

4

3 回答 3

67

只需将其计算为2^(y*log2(x)).

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

于 2011-01-09T09:26:48.657 回答
16

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

代码： 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
于 2011-01-09T18:18:59.567 回答
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()，但我仍然找不到任何东西......

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

于 2019-07-27T03:43:36.817 回答