我一直想知道如何自己制作一个计算功率的函数（例如 23）。在大多数语言中，这些都包含在标准库中，主要是 pow(double x, double y)，但我如何自己编写呢？
我在想for loops，但它认为我的大脑陷入了循环（当我想用非整数指数做幂时，比如 54.5 或负数 2-21) 我疯了;)
那么，我该如何编写一个计算实数幂的函数呢？谢谢
哦，也许需要注意的是：我不能使用使用权力的函数（例如exp），这将使它最终无用。
【问题讨论】：
标签： c++ math floating-point
负能量不是问题，它们只是正能量的倒数 (1/x)。
浮点幂只是稍微复杂一点；如您所知，分数幂等价于根（例如x^(1/2) == sqrt(x)），并且您还知道，具有相同底数的幂的乘法等价于将它们的指数相加。
通过以上所有内容，您可以：
例子：
2^(-3.5) = (2^3 * 2^(1/2)))^-1 = 1 / (2*2*2 * sqrt(2))
【讨论】：
AB = 对数-1(对数(A)*B)
编辑：是的，这个定义确实提供了一些有用的东西。例如，在 x86 上，它几乎直接转换为 FYL2X (Y * Log2(X)) 和 F2XM1 (2x-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 来处理取幂的整数部分。
【讨论】：
log^-1(x) = exp(x) = e^x
通常数学库中 pow(double, double) 函数的实现是基于身份的：
pow(x,y) = pow(a, y * log_a(x))
使用此恒等式，您只需要知道如何将单个数字 a 提高到任意指数，以及如何取对数底 a。您已经有效地将一个复杂的多变量函数变成了一个单变量和一个乘法的两个函数，这很容易实现。 a 最常用的值是e 或2 -- e 因为e^x 和log_e(1+x) 有一些非常好的数学属性，而2 因为它有一些很好的实现属性在浮点运算中。
这样做的好处是（如果您想获得完全准确），您需要计算log_a(x) 项（及其与y 的乘积），使其准确度高于@ 的浮点表示987654334@ 和y。例如，如果 x 和 y 是双精度数，并且您想要获得高精度的结果，则需要想出一些方法以更高精度的格式存储中间结果（并进行算术运算）。 Intel x87 格式是一个常见的选择，64 位整数也是如此（尽管如果你真的想要一个高质量的实现，你需要做几个 96 位整数计算，这在某些情况下有点痛苦语言）。如果您实现powf(float,float)，则处理这个问题要容易得多，因为您可以使用double 进行中间计算。如果您想使用这种方法，我建议您从它开始。
我概述的算法并不是计算pow 的唯一可能方法。它只是最适合提供满足固定先验精度界限的高速结果。它在其他一些情况下不太适合，并且肯定比其他人建议的重复平方[root]-ing 算法更难实现。
如果您想尝试重复平方[root] 算法，请先编写一个仅使用重复平方的无符号整数幂函数。一旦您很好地掌握了该简化情况的算法，您会发现将其扩展为处理小数指数是相当简单的。
【讨论】：
有两种不同的情况需要处理：整数指数和小数指数。
对于整数指数，您可以使用平方取幂。
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) 的良好近似值。一个简单的方法是使用Taylor series。
【讨论】：
根据定义：
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.
实际上，您可以存储1/n! 的前 10 个值的数组，然后近似
exp(x) = 1 + x + x^2/2 + x^3/3! + ... + x^10/10!
因为 10！是一个巨大的数字，所以 1/10！非常小（2.7557319224⋅10^-7）。
【讨论】：
Wolfram functions 提供了多种计算幂的公式。其中一些非常容易实现。
【讨论】：
对于正整数幂，请查看 exponentiation by squaring 和 addition-chain exponentiation。
【讨论】：
使用三个自行实现的函数iPow(x, n)、Ln(x) 和Exp(x)，我能够计算fPow(x, a)、x 和一个双打。下面的函数都没有使用库函数，而只是迭代。
关于实现的功能的一些解释：
(1)iPow(x, n)：x 是double，n 是int。这是一个简单的迭代，因为 n 是一个整数。
(2) Ln(x)：此函数使用泰勒级数迭代。迭代中使用的系列是Σ (from int i = 0 to n) {(1 / (2 * i + 1)) * ((x - 1) / (x + 1)) ^ (2 * n + 1)}。符号^表示在第一个函数中实现的幂函数Pow(x, n)，它使用简单的迭代。
(3) Exp(x)：此函数再次使用泰勒级数迭代。迭代中使用的系列是Σ (from int i = 0 to 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 的所有功能并进行了迭代。
n.b. 我使用了constant MAX_DELTA_DOUBLE 来决定在哪一步停止迭代。我已将其设置为1.0E-15，这对于双打来说似乎是合理的。因此，如果(delta < MAX_DELTA_DOUBLE)，则迭代停止如果您需要更高的精度，您可以使用long 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);
}
【讨论】：
delta < 1.0E-15）。我没有彻底测试过，但是由于 c 中 double 类型可以处理的有效数字的数量最多是 15，我有点相信上面的代码可以完成它的工作假设。
这是一个有趣的练习。以下是一些建议，您应该按此顺序尝试：
【讨论】：
您可以像这样找到 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);
}
}
最好的问候
【讨论】：
有效计算正整数幂的更好算法是重复对基数求平方，同时跟踪额外的余数被乘数。这是一个 Python 示例解决方案，应该相对容易理解并翻译成您的首选语言：
def power(base, exponent):
remaining_multiplicand = 1
result = base
while exponent > 1:
remainder = exponent % 2
if remainder > 0:
remaining_multiplicand = remaining_multiplicand * result
exponent = (exponent - remainder) / 2
result = result * result
return result * remaining_multiplicand
要让它处理负指数，你所要做的就是计算正版本并将 1 除以结果，所以这应该是对上面代码的简单修改。分数指数要困难得多，因为它实际上意味着计算底的 n 次根，其中n = 1/abs(exponent % 1) 并将结果乘以整数部分幂计算的结果：
power(base, exponent - (exponent % 1))
您可以使用牛顿法将根计算到所需的准确度。查看wikipedia article on the algorithm。
【讨论】：
我正在使用定点长算法，并且我的 pow 是基于 log2/exp2 的。数字包括：
int sig = { -1; +1 }签名DWORD a[A+B]号码A 是数字的整数部分的 DWORDs 的数字B 是小数部分的 DWORDs 数我的简化解决方案是这样的：
//---------------------------------------------------------------------------
longnum exp2 (const longnum &x)
{
int i,j;
longnum c,d;
c.one();
if (x.iszero()) return c;
i=x.bits()-1;
for(d=2,j=_longnum_bits_b;j<=i;j++,d*=d)
if (x.bitget(j))
c*=d;
for(i=0,j=_longnum_bits_b-1;i<_longnum_bits_b;j--,i++)
if (x.bitget(j))
c*=_longnum_log2[i];
if (x.sig<0) {d.one(); c=d/c;}
return c;
}
//---------------------------------------------------------------------------
longnum log2 (const longnum &x)
{
int i,j;
longnum c,d,dd,e,xx;
c.zero(); d.one(); e.zero(); xx=x;
if (xx.iszero()) return c; //**** error: log2(0) = infinity
if (xx.sig<0) return c; //**** error: log2(negative x) ... no result possible
if (d.geq(x,d)==0) {xx=d/xx; xx.sig=-1;}
i=xx.bits()-1;
e.bitset(i); i-=_longnum_bits_b;
for (;i>0;i--,e>>=1) // integer part
{
dd=d*e;
j=dd.geq(dd,xx);
if (j==1) continue; // dd> xx
c+=i; d=dd;
if (j==2) break; // dd==xx
}
for (i=0;i<_longnum_bits_b;i++) // fractional part
{
dd=d*_longnum_log2[i];
j=dd.geq(dd,xx);
if (j==1) continue; // dd> xx
c.bitset(_longnum_bits_b-i-1); d=dd;
if (j==2) break; // dd==xx
}
c.sig=xx.sig;
c.iszero();
return c;
}
//---------------------------------------------------------------------------
longnum pow (const longnum &x,const longnum &y)
{
//x^y = exp2(y*log2(x))
int ssig=+1; longnum c; c=x;
if (y.iszero()) {c.one(); return c;} // ?^0=1
if (c.iszero()) return c; // 0^?=0
if (c.sig<0)
{
c.overflow(); c.sig=+1;
if (y.isreal()) {c.zero(); return c;} //**** error: negative x ^ noninteger y
if (y.bitget(_longnum_bits_b)) ssig=-1;
}
c=exp2(log2(c)*y); c.sig=ssig; c.iszero();
return c;
}
//---------------------------------------------------------------------------
地点：
_longnum_bits_a = A*32
_longnum_bits_b = B*32
_longnum_log2[i] = 2 ^ (1/(2^i)) ... precomputed sqrt table
_longnum_log2[0]=sqrt(2)
_longnum_log2[1]=sqrt[tab[0])
_longnum_log2[i]=sqrt(tab[i-1])
longnum::zero() sets *this=0
longnum::one() sets *this=+1
bool longnum::iszero() returns (*this==0)
bool longnum::isnonzero() returns (*this!=0)
bool longnum::isreal() returns (true if fractional part !=0)
bool longnum::isinteger() returns (true if fractional part ==0)
int longnum::bits() return num of used bits in number counted from LSB
longnum::bitget()/bitset()/bitres()/bitxor() are bit access
longnum.overflow() rounds number if there was a overflow X.FFFFFFFFFF...FFFFFFFFF??h -> (X+1).0000000000000...000000000h
int longnum::geq(x,y) is comparition |x|,|y| returns 0,1,2 for (<,>,==)
你只需要理解这段代码就是二进制形式的数字由 2 的幂和组成，当你需要计算 2^num 时，它可以重写为这样
2^(b(-n)*2^(-n) + ... + b(+m)*2^(+m))其中n 是小数位，m 是整数位。二进制形式的2 的乘法/除法是简单的位移，所以如果你把它们放在一起，你会得到类似于我的exp2 的代码。 log2 基于 binaru 搜索...将结果位从 MSB 更改为 LSB，直到它与搜索值匹配（与快速 sqrt 计算非常相似的算法）。希望这有助于澄清事情......
【讨论】：
其他答案中给出了很多方法。这是我认为在积分幂的情况下可能有用的东西。
在 nx 的整数幂 x 的情况下，直接的方法将采用 x-1 乘法。为了优化这一点，我们可以使用动态规划并重用较早的乘法结果来避免所有 x 乘法。例如，在 59 中，我们可以说，批次 3，即计算 53 一次，得到 125，然后立方 125使用相同的逻辑，过程中只进行 4 次乘法运算，而不是直接的 8 次乘法运算。
问题是批次 b 的理想大小是多少，以便乘法次数最少。因此，让我们为此写出方程式。如果 f(x,b) 是表示使用上述方法计算 nx 所需的乘法次数的函数，则
解释：一批 p 个数字的乘积将进行 p-1 次乘法运算。如果我们将 x 次乘法划分为 b 个批次，则每个批次内需要 (x/b)-1 次乘法，所有 b 批次都需要 b-1 次乘法。
现在我们可以计算这个函数关于 b 的一阶导数并将其等同于 0 以获得最少乘法次数的 b。
现在将 b 的值放回函数 f(x,b) 中，以获得最少的乘法次数：
对于所有正数 x，该值小于直接乘法的乘积。
【讨论】：
也许你可以使用泰勒级数展开。函数的泰勒级数是用函数在单个点的导数表示的项的无限总和。对于最常见的函数，函数及其泰勒级数之和在该点附近相等。 Taylor 的系列以 1715 年推出它们的 Brook Taylor 命名。
【讨论】：