【发布时间】:2010-09-14 23:08:44
【问题描述】:
我正在尝试确定我的一种算法的渐近运行时间,该算法使用指数,但我不确定如何以编程方式计算指数。
我正在专门寻找用于双精度浮点数的 pow() 算法。
【问题讨论】:
-
第一次编辑后,这个问题还不清楚。您提到“使用指数的算法”和“用于双精度浮点数的算法”。算法做什么?多个两个这样的数字?计算具有双精度 x-y-z 的 N 个点的 3D 三角剖分?
我正在尝试确定我的一种算法的渐近运行时间,该算法使用指数,但我不确定如何以编程方式计算指数。
我正在专门寻找用于双精度浮点数的 pow() 算法。
【问题讨论】:
我有机会了解 fdlibm 的实现。 cmets 描述了所使用的算法:
* 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** 计算。并且其中任何一个都没有循环。因此,尽管浮点基元很复杂,但它看起来像是一个渐近常数时间算法。
欢迎浮点专家(我不是其中之一)发表评论。 :-)
【讨论】:
**是什么?
x ** y 是pow(x, y),除了2 ** n 的情况,其中n 是一个整数,您可以使用表查找。输入是x(一个浮点数,分为指数n和尾数f)和y。
除非他们发现了更好的方法,否则我相信三角函数、对数函数和指数函数(例如指数增长和衰减)的近似值通常使用算术规则和 泰勒级数 扩展以产生精确到要求精度的近似结果。 (有关幂级数、泰勒级数和麦克劳林级数展开函数的详细信息,请参阅任何微积分书籍。)请注意,我已经有一段时间没有这样做了,所以我不能告诉你,例如,确切地如何计算您需要包含的系列中的项数保证误差小到可以在双精度计算中忽略不计。
例如,e^x 的 Taylor/Maclaurin 级数展开是这样的:
+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 + 分数) * (2^exponent), 1
x * log2(e) = log2(1 + 分数) + 指数, 0
指数 = floor(x * log2(e))
1 + 分数 = 2^(x * log2(e) - 指数) = e^((x * log2(e) - 指数) * ln2) = e^(x - 指数 * ln2), 0
【讨论】: