我该如何自己编写一个求幂函数?

Question

我该如何自己编写一个求幂函数?

c++mathfloating-point

54

我一直在思考如何编写一个计算幂次方的函数（例如 2³）。大多数语言都包含在标准库中，通常为 pow(double x, double y)，但我该如何自己编写它呢？

我想到了使用 for 循环，但是当我尝试进行非整数指数的幂运算时（例如 5^4.5 或负数 2^-21），我的大脑似乎进入了死循环，感觉有些疯狂 ;)

那么，我该如何编写一个计算实数幂次方的函数呢？谢谢！

哦，也许需要注意的是：我不能使用使用幂的函数（例如 exp），这让这个问题变得无用。

- user142019

4

您想基于现实世界处理器可用的特性来实现吗？大多数处理器都包含某种形式的exp(x)硬件指令，因此pow是基于exp和ln来实现的。还是要基于一些缺少这些功能的机器来实现？如果选择后者，请描述所需的指令集。 - Pete Kirkham

可能是如何计算任意幂/根？的重复问题。 - jww

14个回答

23

A^B = Log^-1(Log(A)*B)

编辑：是的，这个定义确实有用。例如，在 x86 上，它几乎直接转换为 FYL2X（Y * Log₂(X)）和 F2XM1（2^x-1）：

fyl2x
fld st(0)
frndint
fsubr st(1),st
fxch st(1)
fchs
f2xmi
fld1
faddp st(1),st
fscale
fstp st(1)

代码看起来比你想象的略长，主要是因为 F2XM1 只能处理范围在 -1.0 到 1.0 之间的数字。 fld st(0)/frndint/fsubr st(1),st 这一步骤减去整数部分，只剩下小数部分。我们对其应用 F2XM1，加回 1，然后使用 FSCALE 来处理指数的整数部分。

- Jerry Coffin

3

一个问题：log^-1(x) = exp(x) = e^x。解释：这个等式的意思是“对数的反函数等于自然指数的指数等于e的指数”。 - user142019

5

有趣的递归定义：对数的倒数也是一个幂，我认为Koning想用更基本的运算（加、减、乘、除等）来计算它。 - fortran

20

通常数学库中pow(double, double)函数的实现是基于以下公式：

pow(x,y) = pow(a, y * log_a(x))

使用这个身份，你只需要知道如何将一个数a提高到任意指数，以及如何取对数基于a。你已经将一个复杂的多变量函数转化为了两个单变量函数和一个乘法，这非常容易实现。最常选择的值是e或2——e因为e^x和log_e(1+x)具有一些非常好的数学性质，2则因为它在浮点运算中具有一些好的性质。

这种方法的缺陷是（如果您想获得完全的精度），您需要计算log_a(x)项（以及其与y的乘积），比浮点表示的x和y都要高的精度。例如，如果x和y是双精度浮点数，并且您希望获得高精度结果，则需要想出一种在更高精度格式中存储中间结果（并进行算术运算）的方法。 Intel x87格式是一个常见的选择，64位整数也是如此（虽然如果您真的想要高质量的实现，则需要执行几个96位整数计算，在某些语言中有一点痛苦）。如果您要使用这种方法，用powf(float,float)来实现会更容易处理，因为那样你只需要使用double类型的中间计算。如果您想使用这种方法，我建议从这个开始。

我概述的算法不是计算pow的唯一可能方法。它只是最适合提供满足固定先验精度限制的高速结果的方式。它在其他一些情况下不太适用，并且当然要比其他人建议的重复平方[root]算法难得多。

如果您想尝试重复平方[root]算法，请先编写一个仅使用重复平方的无符号整数幂函数。一旦你对该缩减案例的算法有了很好的把握，你将发现将其扩展到处理分数指数相当容易。

- Stephen Canon

9

有两种不同的情况需要处理：整数指数和分数指数。

对于整数指数，可以使用平方求幂算法。

def pow(base, exponent):
    if exponent == 0:
        return 1
    elif exponent < 0:
        return 1 / pow(base, -exponent)
    elif exponent % 2 == 0:
        half_pow = pow(base, exponent // 2)
        return half_pow * half_pow
    else:
        return base * pow(base, exponent - 1)

第二个“elif”语句是区别于朴素pow函数的关键。它使函数能够进行O(log n)递归调用，而不是O(n)。

对于分数指数，您可以使用恒等式a^b = C^(b*log_C(a))。方便起见，取C=2，因此a^b = 2^(b * log2(a))。这将问题简化为编写2^x和log2(x)函数。

之所以方便取C=2，是因为浮点数以基数2浮点数存储。log2(a * 2^b) = log2(a) + b。这使得编写log2函数更容易：您不需要让它对每个正数都准确，只需在区间[1,2)上准确即可。类似地，要计算2^x，您可以将2^(x的整数部分) * 2^(x的小数部分)相乘。第一部分在浮点数中存储非常简单，对于第二部分，您只需要在区间[0,1)上拥有一个2^x函数即可。

难点在于找到2^x和log2(x)的良好近似值。一个简单的方法是使用泰勒级数。

- dan04

8

根据定义：

a^b = exp(b ln(a))

其中exp(x) = 1 + x + x^2/2 + x^3/3! + x^4/4! + x^5/5! + ...，

n! = 1 * 2 * ... * n。

在实践中，您可以存储前10个值的数组1/n!，然后近似计算

exp(x) = 1 + x + x^2/2 + x^3/3! + ... + x^10/10!

因为10！是一个巨大的数字，所以1/10！非常小（2.7557319224⋅10^-7）。

- Andreas Rejbrand

2

@Pindatjuh：感谢您尝试使格式更好看。但是，您不小心（我希望如此）完全弄错了公式。当人们写x^n/n!时，他们总是指的是(x^n)/(n!)，而不是x^(n/n!)。 - Andreas Rejbrand

麦克劳林（和泰勒）级数是微积分课程中最有用的东西之一 :-) - fortran

5

Wolfram函数提供了各种用于计算幂的公式。其中一些非常容易实现。

- High Performance Mark

4

对于正整数幂，可以参考平方求幂算法和加法链求幂算法。

- Richie Cotton

2

使用三个自行实现的函数iPow(x, n)、Ln(x)和Exp(x)，我能够计算fPow(x, a)，其中x和a都是双精度浮点数。以下函数均不使用库函数，而只使用迭代。

一些关于已实现函数的解释：

(1) iPow(x, n)：x为双精度浮点数，n为整数。这是一个简单的迭代，因为n是整数。

(2) Ln(x)：此函数使用泰勒级数迭代。在迭代中使用的级数为Σ（从i = 0到n）{(1 / (2 * i + 1)) * ((x - 1) / (x + 1)) ^ (2 * n + 1)}。符号“^”表示第一个函数中实现的幂函数Pow(x, n)，它使用简单的迭代。

(3) Exp(x)：再次使用泰勒级数迭代。在迭代中使用的级数为Σ（从i = 0到n）{x^i / i!}。在这里，“^”表示幂函数，但它不是通过调用第一个Pow(x, n)函数来计算的；相反，在第三个函数中，它与阶乘同时实现，使用d *= x / i。我觉得我必须使用这个技巧，因为在这个函数中，相对于其他函数，迭代需要更多的步骤，而阶乘(i!)大多数情况下会溢出。为了确保迭代不会溢出，在此部分中，幂函数与阶乘并行迭代。通过这种方式，我克服了溢出问题。

(4) fPow(x, a)：x和a都是双精度浮点数。此函数仅调用上述实现的其他三个函数。该函数的主要思想取决于一些微积分：fPow(x, a) = Exp(a * Ln(x))。现在，我已经拥有所有带有迭代的函数iPow、Ln和Exp。

注意：我使用了一个常量MAX_DELTA_DOUBLE来决定在哪一步停止迭代。我将其设置为1.0E-15，这对于双精度浮点数似乎是合理的。因此，如果(delta < MAX_DELTA_DOUBLE)，则迭代停止。如果需要更高的精度，可以使用并减小MAX_DELTA_DOUBLE的常量值，例如减小到1.0E-18（1.0E-18将是最小值）。

这是我用过的代码，可以正常运行。

#define MAX_DELTA_DOUBLE 1.0E-15
#define EULERS_NUMBER 2.718281828459045

double MathAbs_Double (double x) {
    return ((x >= 0) ? x : -x);
}

int MathAbs_Int (int x) {
    return ((x >= 0) ? x : -x);
}

double MathPow_Double_Int(double x, int n) {
    double ret;
    if ((x == 1.0) || (n == 1)) {
        ret = x;
    } else if (n < 0) {
        ret = 1.0 / MathPow_Double_Int(x, -n);
    } else {
        ret = 1.0;
        while (n--) {
            ret *= x;
        }
    }
    return (ret);
}

double MathLn_Double(double x) {
    double ret = 0.0, d;
    if (x > 0) {
        int n = 0;
        do {
            int a = 2 * n + 1;
            d = (1.0 / a) * MathPow_Double_Int((x - 1) / (x + 1), a);
            ret += d;
            n++;
        } while (MathAbs_Double(d) > MAX_DELTA_DOUBLE);
    } else {
        printf("\nerror: x < 0 in ln(x)\n");
        exit(-1);
    }
    return (ret * 2);
}

double MathExp_Double(double x) {
    double ret;
    if (x == 1.0) {
        ret = EULERS_NUMBER;
    } else if (x < 0) {
        ret = 1.0 / MathExp_Double(-x);
    } else {
        int n = 2;
        double d;
        ret = 1.0 + x;
        do {
            d = x;
            for (int i = 2; i <= n; i++) {
                d *= x / i;
            }
            ret += d;
            n++;
        } while (d > MAX_DELTA_DOUBLE);
    }
    return (ret);
}

double MathPow_Double_Double(double x, double a) {
    double ret;
    if ((x == 1.0) || (a == 1.0)) {
        ret = x;
    } else if (a < 0) {
        ret = 1.0 / MathPow_Double_Double(x, -a);
    } else {
        ret = MathExp_Double(a * MathLn_Double(x));
    }
    return (ret);
}

- ssd

不要使用泰勒级数来近似一个区间上的函数：http://lolengine.net/blog/2011/12/21/better-function-approximations - Pascal Cuoq

我已经阅读了链接中的文档。虽然文档中的想法依赖于一些科学结果，但是我代码中的迭代最多需要15到20步（直到delta < 1.0E-15）。我没有进行彻底测试，但由于在C语言中double类型可以处理的有效数字（小数点后）最多为15，我有点相信上面的代码会按照预期完成工作。 - ssd

1

您可以像这样找到 pow 函数：

static double pows (double p_nombre, double p_puissance)
{
    double nombre   = p_nombre;
    double i=0;
    for(i=0; i < (p_puissance-1);i++){
          nombre = nombre * p_nombre;
       }
    return (nombre);
}

你可以像这样找到floor函数：

static double floors(double p_nomber)
{
    double x =  p_nomber;
    long partent = (long) x; 

    if (x<0)
    {
        return (partent-1);
    }
    else
    {
        return (partent);
    }
}

最好的祝福

- che.moor

这里实现的pows函数仅处理整数指数。 - programmerjake

1

这是一个有趣的练习。以下是一些建议，你可以按照以下顺序尝试：

使用循环。
使用递归（虽然不一定更好，但同样有趣）
通过使用分治技术大幅优化递归
使用对数

- Wouter Lievens

你知道有哪些算法存在吗？谢谢。 - user142019

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- fortran · Accepted Answer

负数幂不是问题，它们只是正数幂的倒数 (1/x)。

浮点数幂稍微有些复杂；如你所知，分数幂等价于根式 (例如，x^(1/2) == sqrt(x))，而且你也知道，将具有相同底数的幂乘起来等价于将它们的指数相加。

有了以上所有这些，你可以：

将指数分解为整数部分和有理部分。
使用循环计算整数幂 (可以通过分解因子并重复使用部分计算进行优化)。
使用任何你喜欢的算法 (像二分法或牛顿法这样的迭代逼近方法都可以工作) 计算根式。
将结果相乘。
如果指数为负数，则应用其倒数。

示例：

2^(-3.5) = (2^3 * 2^(1/2)))^-1 = 1 / (2*2*2 * sqrt(2))