这个反编译的 f2xm1/fscale 序列是用来做什么的?

what's this decompiled f2xm1/fscale sequence meant to do?

我正在尝试对一些最初用 C/C++ 编写的反编译代码进行逆向工程,即我怀疑下面的 FPU 相关代码序列可能来自一些简单的 C 代码“double”在生成的汇编代码中处理看起来更复杂。在此之前,已经执行了一些浮点乘法,结果在 ST0(对应于 d1)中。我已经阅读了有关底层 FPU 操作在技术上的作用的文档,但相应代码序列的意图对我来说仍然不是很明显。

  d1 = (float10)1.442695040888963 *(float10)0.6931471805599453 * (float10)DOUBLE_00430088 * (float10)param_1[0x58];
  d2 = ROUND(d1);
  d1 = (float10)f2xm1(d1 - d2);
  d1 = (float10)fscale((float10)1 + d1,d2);

对原始 d1 结果的预期更改是什么,即 C 代码对原始 d1“double”做了什么?

PS:低于实际的 x86 代码(以防 Ghidra 在反编译时误解某些内容):


                     * Push ST(i) onto the FPU register stack, i.e.               *
                     * duplicate ST0                                              *
                     **************************************************************
0040aeef d9 c0           FLD        ST0
                     **************************************************************
                     * Rounds the source value in the ST(0) register              *
                     * to the nearest integral value, depending on                *
                     * the current rounding mode (lets suppose some               *
                     * "floor" mode was used?)                                    *
                     **************************************************************
0040aef1 d9 fc           FRNDINT
                     **************************************************************
                     * Exchange the contents of ST(0) and ST(1).                  *
                     **************************************************************
0040aef3 d9 c9           FXCH
                     **************************************************************
                     * get fractional part?                                       *
                     **************************************************************
0040aef5 d8 e1           FSUB       d2[0],d1[0]
                     **************************************************************
                     * Computes the exponential value of 2 to the power           *
                     * of the source operand minus 1.                             *
                     **************************************************************
0040aef7 d9 f0           F2XM1
                     **************************************************************
                     * Push +1.0 onto the FPU register stack.                     *
                     **************************************************************
0040aef9 d9 e8           FLD1
                     **************************************************************
                     * Add ST(0) to ST(1), store result in ST(1),                 *
                     * and pop the register stack.                                *
                     **************************************************************
0040aefb de c1           FADDP
                     **************************************************************
                     * Scale ST(0) by ST(1).  This instruction provides           *
                     * rapid multiplication or division by integral               *
                     * powers of 2.                                               *
                     **************************************************************
0040aefd d9 fd           FSCALE
                     **************************************************************
                     * The FSTP instruction copies the value in the ST(0)         *
                     * register to the destination operand and then pops          *
                     * the register stack.                                        *
                     **************************************************************
0040aeff dd d9           FSTP       d1[0]

似乎是 pow(x,y) 实现的一些变体(参见 How can I write a power function myself?)。 Ghidra 在反编译代码视图中把它弄得一团糟。

在调试器中跟踪结果,执行的功能确实是:

pow((float10)DOUBLE_00430088, (float10)param_1[0x58])

没有完整的上下文很难完全确定,但这里的计算似乎是通过 F2XM1exp(x) 的天真计算。请注意,F2XM1 的参数在早期的 x87 实现中被限制为 [-0.5, 0.5],在后来的实现中被限制为 [-1, 1]。用 FRNDINT 计算参数的整数部分并从参数中减去它会产生适合 F2XM1 使用的小数部分。该代码假定默认的舍入模式,舍入到最近或偶数,是有效的。

下面是用于简单 exp(x) 计算的完整带注释的指令序列,已编程为尽可能接近反汇编代码。我使用了 Intel C/C++ 编译器的内联汇编工具,所以它使用了 Intel 语法。在注释中,插入符号 ^ 表示求幂。这个程序的输出是:

x=2.5000000000000000e+000 exp(x)=1.2182493960703475e+001 lib=1.2182493960703473e+001
#include <stdio.h>
#include <math.h>

int main (void)
{
    double r, x = 2.5; 

    __asm fld qword ptr[x]; // x
    __asm fldl2e;           // log2(e)=1.442695040888963, x
    __asm fmulp st(1), st;  // x*log2(e)
    __asm fld st(0);        // x*log2(e), x*log2(e)
    __asm frndint;          // rint(x*log2(e)), x*log2(e)
    __asm fxch st(1);       // x*log2(e), rint(x*log2(e))
    __asm fsub st, st(1);   // frac(x*log2(e)), rint(x*log2(e))
    __asm f2xm1;            // 2^frac(x*log2(e))-1, rint(x*log2(e))
    __asm fld1;             // 1, 2^frac(x*log2(e))-1, rint(x*log2(e))
    __asm faddp st(1),st;   // 2^frac(x*log2(e)), rint(x*log2(e))
    __asm fscale;           // exp(x)=2^frac(x*log2(e))*2^rint(x*log2(e)), rint(x*log2(e)
    __asm fstp qword ptr[r];// rint(x*log2(e)
    __asm fstp st(0);       // <empty>
    
    printf ("x=%23.16e exp(x)=%23.16e lib=%23.16e\n", x, r, exp(x));
    return 0;
}