我正在尝试确定我的一种算法的渐近运行时间,该算法使用指数,但我不确定如何以编程方式计算指数。
我专门寻找用于双精度浮点数的 pow() 算法。
我有机会了解 fdlibm 的实现。评论描述了所使用的算法:
* n
* Method: Let x = 2 * (1+f)
* 1. Compute and return log2(x) in two pieces:
* log2(x) = w1 + w2,
* where w1 has 53-24 = 29 bit trailing zeros.
* 2. Perform y*log2(x) = n+y' by simulating muti-precision
* arithmetic, where |y'|<=0.5.
* 3. Return x**y = 2**n*exp(y'*log2)
随后列出所有已处理的特殊情况(0、1、inf、nan)。
在所有特殊情况处理之后,代码中最密集的部分涉及
log2
和 2**
计算。并且其中任何一个都没有循环。因此,尽管浮点原语很复杂,但它看起来像是一个渐进恒定时间算法。
欢迎浮点专家(我不是其中之一)发表评论。 :-)
除非他们发现了更好的方法,否则我相信三角函数、对数函数和指数函数的近似值(例如指数增长和衰减)通常是使用算术规则和 泰勒级数 展开来计算的,以产生近似结果精确到所要求的精度内。 (有关幂级数、泰勒级数和麦克劳林级数展开函数的详细信息,请参阅任何微积分书籍。)请注意,我已经有一段时间没有做这些了,所以我无法告诉您,例如,确切地如何计算您需要包含的级数中的项数保证误差足够小,在双精度计算中可以忽略不计。
例如,e^x 的泰勒/麦克劳林级数展开式是这样的:
+inf [ x^k ] x^2 x^3 x^4 x^5
e^x = SUM [ --- ] = 1 + x + --- + ----- + ------- + --------- + ....
k=0 [ k! ] 2*1 3*2*1 4*3*2*1 5*4*3*2*1
如果采用所有项(k 从 0 到无穷大),则此扩展是精确且完整的(没有错误)。
但是,如果您不将所有项都取为无穷大,而是在 5 项或 50 项或其他项之后停止,您将产生一个与实际 e^x 函数值相差一个余数的近似结果:相当容易计算。
指数的好消息是它收敛得很好,并且其多项式展开式的项很容易迭代编码,所以你可能(重复,可能 - 记住,已经有一段时间了)甚至不需要预先-计算需要多少项来保证误差小于精度,因为您可以测试每次迭代的贡献大小,并在它变得足够接近零时停止。实际上,我不知道这个策略是否可行——我必须尝试一下。有一些重要的细节我早已忘记了。例如:机器精度、机器误差和舍入误差等
另外,请注意,如果您不使用 e^x,而是使用另一个基数(如 2^x 或 10^x)进行增长/衰减,则近似多项式函数会发生变化。
对于整数指数,将 a 提升到 b 的通常方法是这样的:
result = 1
while b > 0
if b is odd
result *= a
b -= 1
b /= 2
a = a * a
指数的大小一般是对数。该算法基于不变式“a^b*result = a0^b0”,其中a0和b0是a和b的初始值。
对于负数或非整数指数,需要对数和近似以及数值分析。运行时间将取决于所使用的算法以及库的调整精度。
编辑:由于似乎有一些兴趣,这里有一个没有额外乘法的版本。
result = 1
while b > 0
while b is even
a = a * a
b = b / 2
result = result * a
b = b - 1
您可以使用 exp(n*ln(x)) 来计算 xn。 x 和 n 都可以是双精度浮点数。自然对数和指数函数可以使用泰勒级数计算。您可以在这里找到公式:http://en.wikipedia.org/wiki/Taylor_series
如果我正在编写一个针对 Intel 的 pow 函数,我将返回 exp2(log2(x) * y)。英特尔的 log2 微代码肯定比我能够编写的任何代码都要快,即使我还记得我第一年的微积分和研究生院数值分析。
e^x = (1 + fraction) * (2^exponent), 1 <= 1 + fraction < 2
x * log2(e) = log2(1 + fraction) + exponent, 0 <= log2(1 + fraction) < 1
exponent = floor(x * log2(e))
1 + fraction = 2^(x * log2(e) - exponent) = e^((x * log2(e) - exponent) * ln2) = e^(x - exponent * ln2), 0 <= x - exponent * ln2 < ln2