在浮点中实现pow()函数的最有效方法

Question

我正在尝试实现我自己版本的 pow() 和 sqrt() 函数，因为我的自定义库没有 pow()/sqrt() 浮点支持。

有人可以帮忙吗？

Answer 1

是的，Sun 可以（我猜现在是 Oracle）：

fdlibm，“可自由分发的数学库”，具有 sqrt 和 pow，以及许多其他数学函数。

不过，它们是相当高科技的实现，当然，没有什么是像这样的“最有效”的实现。您是在寻找源代码来完成它，还是您真的不是在寻找

pow

和

sqrt

，而是实际上在寻找浮点算法编程教育？

Answer 2

当然 - 如果你有指数函数和自然对数函数，这很容易。

从

y = x^n

开始，可以取两边的自然对数：

ln(y) = n*ln(x)

然后取两边的指数就可以得到你想要的：

y = exp(n*ln(x))

如果你想要更好的东西，我知道最好的地方是Abramowitz 和 Stegun。

Answer 3

请注意，如果您的指令集有平方根或幂指令，那么使用它会更好。例如，x87 浮点指令有一条指令

fsqrt

，而 SSE2 添加的指令包括另一条指令

sqrtsd

，这可能比大多数用 C 编写的解决方案要快得多。事实上，至少 gcc 使用这两个指令在 x86 机器上进行编译时的说明。

然而，对于权力来说，事情变得有些模糊。 x87浮点指令集中有一条指令可以用来计算n*log2(n)，即

fyl2x

。另一条指令

fldl2e

将 log2(e) 存储在浮点堆栈中。您可能想看看这些。

您可能还想看看各个 C 库是如何做到这一点的。例如，

dietlibc

，只需使用

fsqrt

:

sqrt:
    fldl 4(%esp)
    fsqrt
    ret

glibc

对于硬件平方根指令不可用的机器（在

sysdeps/ieee754/flt-32/e-sqrtf.c

下）使用 Sun 的实现，并在 x86 指令集上使用

fsqrt

（尽管可以指示 gcc 改为使用

sqrtsd

指令）。 )

Answer 4

平方根是通过迭代牛顿法正确实现的。

Answer 5

double ipow(int base, int exp)
{

bool flag=0;
if(exp<0) {flag=1;exp*=-1;}
int result = 1;
while (exp)
{
    if (exp & 1)
        result *= base;
    exp >>= 1;
    base *= base;
}
if(flag==0)
return result;
else
return (1.0/result);
}
//most suitable way to implement power function for integer to power integer

Answer 6

为了计算 C 中浮点数的平方根，如果您的目标是 x86，我建议使用

fsqrt

。 您可以使用这样的 ASM 指令：

asm("fsqrt" : "+t"(myfloat));

对于海湾合作委员会或

asm {

fstp myfloat

fsqrt

fldp myfloat

}

或者类似 Visual Studio 的东西。

为了实现 pow，应该使用像 upitasoft.com/link/powLUT.h 中那样的大型 switch 语句。 它可能会导致一些缓存问题，但如果你保持这样，它应该不是问题，只需限制范围（注意，你仍然可以优化我提供的代码）。

如果你想支持浮点运算，那就更难了...... 您可以尝试使用自然对数和指数函数，例如：

float result = exp(number * log(power));

但通常它很慢和/或不精确。

希望我有帮助。

Answer 7

下面是用于处理单精度数字的标准 powf() 函数的快速模拟代码。该代码适用于 Visual Studio C++、x86 架构。它适用于整个输入数据范围。误差幅度为几个 ULP，性能（在我的系统上）比标准实现高 12 倍以上。请注意，这不是 powf() 函数的精确模拟：某些退化情况下的行为是不同的。更详细的描述以及双精度数的类似函数的代码可以在并行论坛帖子中找到：对实数求幂的最快方法是什么？。

_declspec(naked) float _vectorcall pwrf(float x, float y)
{
  static const float ct[9] =       // Table of constants
  {
    // Binary logarithm constants
    -1.0f,                         // -1
    0.576379478f,                  // b2*5^5
    0.961804628f,                  // b1*5^3
    2.88539004f,                   // b0*5 and b0 for exponentiation
    0.438455284f,                  // b3*5^7
    // Binary exponential function constants
    0.115524434f,                  // b1*2^2
    -9.21293359E-4f,               // b2*2^4
    // Exponent-correcting coefficients
    2.0f,                          // 2
    2.98023224E-8f                 // 2^-25
  };
  _asm
  {
    vmovd eax,xmm0                 // Read the binary representation of x into eax
    vcvtss2sd xmm4,xmm1,xmm1       // xmm4 = y in double format
    mov edx,offset ct              // edx - constant table address
    cmp eax,0x7F800000             // Compare x with  Inf
    jnc lb_bad_x                   // Abnormal branch if x<=-0, x=+Inf or x=+NaN
    ror eax,23                     // Shift eax so that the exponent is in al
    movzx ecx,al                   // Get the exponent of a number in ecx with offset +127
    jecxz lb_denorm                // Goto if x is denormal, incl. x=0
      lb_calc:                     // Entry point after number normalization
    setnc al                       // al=1 if the mantissa is less than 1.5, otherwise al=0
    adc ecx,-127                   // ecx = k, where k is the integer part of the logarithm
    or eax,126                     // Change the order of x so that it becomes 0.75<=x<1.5
    ror eax,9                      // In eax - the binary representation of the reduced value of x
    vmovd xmm0,eax                 // In xmm0 - reduced value of x
    vsubss xmm1,xmm0,[edx]         // xmm1 = x+1
    vaddss xmm0,xmm0,[edx]         // xmm0 = x-1
    vcvtsi2sd xmm3,xmm3,ecx        // xmm3 = k (в формате double) - добавка к двоичному логарифму
    vdivss xmm0,xmm0,xmm1          // xmm0 = (x-1)/(x+1)=t/5
    vmovss xmm1,[edx+16]           // xmm1 = 5^7*b3 - initialize sum
    vmulss xmm2,xmm0,xmm0          // xmm2 = t^2/25 - prepare polynomial argument
    vfmadd213ss xmm1,xmm2,[edx+4]  // Step 1 of calculating a polynomial using Horner's scheme
    vfmadd213ss xmm1,xmm2,[edx+8]  // Step 2 of calculating a polynomial using Horner's scheme
    vfmadd213ss xmm1,xmm2,[edx+12] // Step 3 of calculating a polynomial using Horner's scheme
    vmulss xmm0,xmm0,xmm1          // Get the fractional part of the logarithm in xmm0
    jmp lb_done                    // Continue to the general part of the function
      lb_denorm:                   // Processing of denormalized x values, incl. x=+0
    bsr ecx,eax                    // Finding the highest setted bit; zf=1 if x=0
    ror eax,cl                     // Form the mantissa of the normalized number x*2^31
    lea ecx,[ecx-31]               // 31 is added to the order, so subtract 31 from ecx
    jnz lb_calc                    // Go to calculate lb(x) if x>0
    vmulss xmm0,xmm0,[edx]         // xmm0 = -0
      lb_bad_x:                    // Branch for handling abnormal x values
    vxorpd xmm3,xmm3,xmm3          // xmm3 = 0
    jnl lb_done                    // If x=+Inf or +NaN leave the value x
    vsqrtss xmm0,xmm0,xmm0         // The root of a negative number gives the result NaN
    vrcpss xmm0,xmm0,xmm0          // In case x=-0, generate the result -Inf
      lb_done:                     // In xmm0 the logarithm value is obtained
    vcvtss2sd xmm0,xmm0,xmm0       // Convert xmm0 to double
    vaddsd xmm0,xmm0,xmm3          // Add the integer part, now xmm0 = lb(x)
    vmulsd xmm0,xmm0,xmm4          // xmm0 = y*lb(x) - binary exponent
    vmovq xmm4,[edx+28]            // Read the correction coefficients into xmm4
    vroundsd xmm1,xmm0,xmm0,0      // xmm1 = k - exponent rounded to integer in double format
    vcvtsd2si eax,xmm1             // eax = k
    dec eax                        // of=1 if edx=80000000h (i.e. overflow or NaN)
    jno exp_cont                   // Goto form 2^k if |y*lb(x)| not too big
    vucomisd xmm0,xmm4             // Test the sign of xmm0 and identify the value of NaN
    jnp exp_break                  // Next, 0 is formed (if y*lb(x)<0) or Inf (if y*lb(x)>0)
    xor eax,eax                    // When forming NaN, make the incremental exponent equal to 0
      exp_cont:                    // Continue for not too large |y*lb(x)|
    vsubsd xmm1,xmm0,xmm1          // xmm1 = x*lb(e)-k = t/2 - the value from -0.5 to 0.5
    vmovss xmm3,[edx+24]           // Initialize the sum with the minor coefficient
    vcvtsd2ss xmm1,xmm1,xmm1       // xmm1 = t/2 in float format
    add eax,127                    // In eax - incremental exponent 2^(k-1)
    jg exp_test_big                // Skip if the exponent is not too small
    vmovshdup xmm4,xmm4            // vmm4 = 2^-25
    sub eax,-26                    // Incremental exponent is k+25 for supposedly denormalized result
    jle exp_break                  // Goto if x<<0
      exp_test_big:                // Checking x>>0
    vmulss xmm2,xmm1,xmm1          // xmm2 = t^2/4 - prepare polynomial argument
    cmp eax,254                    // Check that the exponent is not too large
    jg exp_break                   // Go to set Inf in case of overflow
    vfmadd213ss xmm3,xmm2,[edx+20] // xmm3 = 4*b1+4*b2*t^2
    shl eax,23                     // Binary representation of a number 2^(k-1) or 2^(k+31)
    vfmsub213ss xmm3,xmm2,xmm1     // xmm3 = -t/2+b1*t^2+b2*t^4
    vaddss xmm0,xmm1,xmm1          // xmm0 = t
    vaddss xmm3,xmm3,[edx+12]      // xmm3 = b0-t/2+b1*t^2+b2*t^4 = f(t)-t/2
    vmovd xmm1,eax                 // xmm1 = 2^(k-1) или 2^(k+25)
    vdivss xmm0,xmm0,xmm3          // xmm0 = t/(f(t)-t/2)
    vfmadd213ss xmm0,xmm1,xmm1     // xmm0 = 2^x with the exponent shifted by -1 or 25
    vmulss xmm0,xmm0,xmm4          // xmm0 = 2^x
    ret                            // Return
      exp_break:                   // Next, 0 is formed (if y*lb(x)<0) or Inf (if y*lb(x)>0)
    vxorps xmm0,xmm0,xmm0          // Initialize result to zero
    jbe exp_end                    // Go to the end if the result is really close to 0
    vrcpss xmm0,xmm0,xmm0          // In the case of y*lb(x)>0, place infinity in xmm0 (Inf)
      exp_end:                     // The result in xmm0 is ready
    ret                            // Return
  }
}

Answer 8

我能想到的执行 pow() 的最快方法就是沿着这些思路（注意，这非常复杂）：

//raise x^y
double pow(double x, int y) {
    int power;
    map<int, double> powers;
    for (power = 1; power < y; power *= 2, x *= x)
        powers.insert(power, x);

    while (power > y) {
        //figure out how to get there
        map<int, double>::iterator p = powers.lower_bound(power - y);
        //p is an iterator that points to the biggest power we have that doesn't go over power - y
        power -= p->first;
        x /= p->second;
    }

    return x;
}

我不知道如何实现十进制幂。我最好的猜测是使用对数。

编辑：我正在尝试使用对数解决方案（基于 y），而不是您建议的线性解决方案。让我解决这个问题并编辑它，因为我知道它有效。

编辑2：呵呵，我的错。 power *= 2 而不是 power++