腾讯云

文章/答案/技术大牛

发布

社区首页 >问答首页 >确定浮点平方根

问确定浮点平方根
EN

Stack Overflow用户

提问于 2012-02-10 22:04:55

回答 4查看 10.7K关注 0票数 7

如何确定浮点数的平方根？牛顿-拉夫森方法是一种好方法吗？我也没有硬件平方根。我也没有硬件除法(但我实现了浮点除法)。

如果可能的话，我更愿意尽可能地减少分割的数量，因为它们是如此昂贵。

另外，减少总迭代次数的初始猜测应该是什么？

非常感谢!

math

floating-point

computer-science

algorithm

回答 4

Stack Overflow用户

回答已采纳

发布于 2012-02-10 23:30:04

当您使用Newton-Raphson来计算平方根时，您实际上希望使用迭代来找到平方根的倒数(之后，您可以简单地乘以输入--需要注意舍入--以产生平方根)。

更准确地说:我们使用函数f(x) = x^-2 - n。显然，如果是f(x) = 0，那么就是x = 1/sqrt(n)。这就产生了牛顿迭代：

x_(i+1) = x_i - f(x_i)/f'(x_i)
        = x_i - (x_i^-2 - n)/(-2x_i^-3)
        = x_i + (x_i - nx_i^3)/2
        = x_i*(3/2 - 1/2 nx_i^2)

请注意(与平方根的迭代不同)，这种平方根倒数的迭代不涉及除法，因此它通常效率更高。

我在你关于divide的问题中提到，你应该看看现有的软浮点库，而不是重新发明轮子。这条建议在这里也适用。这个函数已经在现有的软浮点库中实现了。

编辑:提问者似乎仍然很困惑，所以让我们来做一个例子：sqrt(612)。612是1.1953125 x 2^9 (或者b1.0011001 x 2^9，如果您更喜欢二进制)。取出指数(9)的偶数部分，将输入写为f * 2^(2m)，其中m是整数，f在[1，4]范围内。然后我们将拥有：

sqrt(n) = sqrt(f * 2^2m) = sqrt(f)*2^m

将这种简化应用到我们的示例中，可以得到f = 1.1953125 * 2 = 2.390625 (b10.011001)和m = 4。现在进行牛顿-拉普森迭代，使用起始猜测0.5 (正如我在评论中提到的，此猜测对所有f都收敛，但使用线性近似作为初始猜测可以做得更好)：

x_0 = 0.5
x_1 = x_0*(3/2 - 1/2 * 2.390625 * x_0^2)
    = 0.6005859...
x_2 = x_1*(3/2 - 1/2 * 2.390625 * x_1^2)
    = 0.6419342...
x_3 = 0.6467077...
x_4 = 0.6467616...

因此，即使有一个(相对较差的)初始猜测，我们也能快速收敛到1/sqrt(f) = 0.6467616600226026的真实值。

现在我们简单地将最终结果组合在一起：

sqrt(f) = x_n * f = 1.5461646...
sqrt(n) = sqrt(f) * 2^m = 24.738633...

并检查: sqrt(612) = 24.738633...

显然，如果你想要正确的舍入，需要仔细的分析，以确保你在计算的每个阶段都有足够的精度。这需要仔细的记账，但这不是火箭科学。您只需保持仔细的误差界限，并通过算法传播它们。

如果您希望在不显式检查残差的情况下校正舍入，则需要将sqrt(f)计算为2p +2位的精度(其中p是源和目标类型的精度)。但是，您也可以采用计算sqrt(f)的策略，使其比p位多一点，对该值进行平方，并在必要时逐个调整尾随位(这通常更便宜)。

sqrt很好，因为它是一个一元函数，这使得在商用硬件上对单精度进行详尽的测试是可行的。

您可以在opensource.apple.com上找到OS X软浮点sqrtf函数，它使用上述算法(碰巧是我写的)。它是根据APSL许可的，这可能适合也可能不适合您的需求。

票数 11

Stack Overflow用户

发布于 2012-02-11 05:28:04

可能(仍然)是查找inverse square root和我最喜欢的10行代码的最快实现。

它是基于牛顿近似的，但有一些奇怪的地方。甚至有一个围绕这个的great story。

票数 4

Stack Overflow用户

发布于 2012-02-10 22:11:14

最容易实现的(你甚至可以在计算器中实现)：

def sqrt(x, TOL=0.000001):
    y=1.0
    while( abs(x/y -y) > TOL ):
        y= (y+x/y)/2.0
    return y

这完全等同于牛顿·拉普森：

Y(新)=y- f(y)/f'(y)

f(y) = y^2-x和f'(y) = 2y

替换这些值：

Y(新)=y- (y^2-x)/2y = (y^2+x)/2y = (y+x/y)/2

如果除法成本很高，你应该考虑：http://en.wikipedia.org/wiki/Shifting_nth-root_algorithm。

移位算法：

让我们假设你有两个数字a和b，使得最低有效位(等于1)大于b，而b只有一个比特等于(例如：a=1000和b=10)。设s(b) = log_2(b) (它只是b中值为1的位元的位置)。

假设我们已经知道a^2的值。现在(a+b)^2 = a^2 + 2ab + b^2。a^2已经已知，2ab:将a移位s(b)+1，b^2:将b移位s(b)。

算法：

Initialize a such that a has only one bit equal to one and a^2<= n < (2*a)^2. 
Let q=s(a).    
b=a
sqra = a*a

For i = q-1 to -10 (or whatever significance you want):
    b=b/2
    sqrab = sqra + 2ab + b^2
    if sqrab > n:
        continue
    sqra = sqrab
    a=a+b

n=612
a=10000 (16)

sqra = 256

Iteration 1:
    b=01000 (8) 
    sqrab = (a+b)^2 = 24^2 = 576
    sqrab < n => a=a+b = 24

Iteration 2:
    b = 4
    sqrab = (a+b)^2 = 28^2 = 784
    sqrab > n => a=a

Iteration 3:
    b = 2
    sqrab = (a+b)^2 = 26^2 = 676
    sqrab > n => a=a

Iteration 4:
    b = 1
    sqrab = (a+b)^2 = 25^2 = 625
    sqrab > n => a=a

Iteration 5:
    b = 0.5
    sqrab = (a+b)^2 = 24.5^2 = 600.25
    sqrab < n => a=a+b = 24.5

Iteration 6:
    b = 0.25
    sqrab = (a+b)^2 = 24.75^2 = 612.5625
    sqrab < n => a=a


Iteration 7:
    b = 0.125
    sqrab = (a+b)^2 = 24.625^2 = 606.390625
    sqrab < n => a=a+b = 24.625

and so on.