在浮点中实现 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 and 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 将 Sun 的实现用于硬件平方根指令不可用的机器（在 sysdeps/ieee754/flt-32/e-sqrtf.c 下），并在 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));

对于 GCC 或

asm {

fstp myfloat

fsqrt

fldp myfloat

}

或类似 Visual Studio 的东西。

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

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

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

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

希望我能帮上忙。

Answer 7

我能想到执行 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++