浮点数平方根的确定

Question

我该如何求浮点数的平方根？牛顿-拉夫逊法是一个好方法吗？我也没有硬件平方根。虽然我已经实现了浮点数除法，但如果可能的话，我想尽可能减少除法的数量，因为它们非常昂贵。

另外，为了减少总迭代次数，初始猜测应该是多少？

非常感谢!

- Veridian

---

是的，牛顿-拉夫逊法在找平方根方面速度很快。 - ypercubeᵀᴹ

更适合于math.stackexchange。 - Andrey Rubshtein

我只是在问问题，没有必要这么粗鲁。 - Fred Foo

Answer 1

当你使用牛顿-拉弗森方法计算平方根时，实际上是希望使用迭代来寻找倒数平方根（之后你可以通过乘以输入的值（需要注意舍入）来得到平方根）。

更准确地说：我们使用函数 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)

请注意（与求平方根的迭代不同），这个求倒数平方根的迭代不涉及除法，因此通常更加高效。

我在你的除法问题中提到过，你应该查看现有的soft-float库，而不是重新发明轮子。这个函数已经在现有的soft-float库中实现了。

编辑：提问者似乎仍然困惑，那么我们来做一个例子：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。现在进行牛顿-拉夫逊迭代，以找到 x = 1/sqrt(f)，使用起始猜测值为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是源类型和目标类型的精度）。然而，您还可以采取计算比p位多一点的sqrt(f)的策略，对该值进行平方，如果必要的话，将尾数位调整为1（这通常更便宜）。

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

您可以在opensource.apple.com上找到OS X软浮点sqrtf函数，该函数使用上述算法（恰好是我写的）。它受APSL许可证的管辖，这可能适合或不适合您的需求。

- Stephen Canon

---

如我所解释的，使用 x^-2 - n 可避免在迭代中进行除法运算，在保留二次收敛的同时，使其在典型硬件上更加高效。 - Stephen Canon

@starbox: 我相信你可以从公共代码库中获取libgcc和compiler_rt的源码。compiler_rt也可以在网上浏览：http://llvm.org/viewvc/llvm-project/compiler-rt/trunk/lib/divsf3.c?revision=144660&view=markup - Stephen Canon

@starbox: 612 = 2.390625 x 2^8，因此 sqrt(612) = sqrt(2.390625) x sqrt(2^8) = sqrt(2.390625) x 2^4。提取指数的偶数部分，你会得到一个值在 [1,4) 范围内。 - Stephen Canon

@starbox：类似于1.1 - 1/6 * f的线性近似。（1/6是最优的，你可以解出最优的第一项；我想不起来确切的值，但它几乎等于1.1） - Stephen Canon

@StephenCanon，只是为了澄清，您是如何确定第一个最优术语的？ - Veridian

Answer 2

可能仍然是查找倒数平方根最快的实现，这是我最喜欢的10行代码。

它基于牛顿逼近法，但有一些小问题。周围甚至有一个很棒的故事。

- emboss

---

链接无效。可能已移至 https://www.beyond3d.com/content/articles/8/？ - treehouse

这实际上是一篇不同的博客文章，早于2年。https://web.archive.org/web/20111126121744/https://codemaestro.com/reviews/9 是2011年的存档副本。 - Adam Hyland

Answer 3

最容易实现的（你甚至可以在计算器上实现）：

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

这与牛顿-拉弗森法完全相等：

y(new) = y - f(y)/f'(y)

其中，f(y) = y^2-x，f'(y) = 2y

代入这些值：

y(new) = y - (y^2-x)/2y = (y^2+x)/2y = (y+x/y)/2

如果除法的开销很大，您应该考虑使用移位根算法。移位算法：

假设您有两个数字a和b，使得最低有效位（等于1）大于b，而b只有一位等于1（例如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.

- ElKamina

---

这基本上是交叉法的一种变体，预计收敛速度比牛顿-拉夫逊方法慢得多。 - amit

@ElKamina 虽然我已经实现了浮点数除法，但是我希望能够减少执行次数，因为它非常昂贵。我使用的处理器没有硬件浮点数除法单元。 - Veridian

@amit 这个方法等同于牛顿-拉夫逊法。 - ElKamina

@starbox 你的乘法效率与加法或减法相比如何？也就是说，多少次加法等于乘法？这是否取决于被乘数的值？ - ElKamina

@ElKamina，你能举个例子说明如何使用移位n次根算法来处理浮点数吗？我看了那个页面，但它似乎不容易实现。顺便说一下，我有硬件移位器。 - Veridian

Answer 4

在区间[1,4)上，一个很好的平方根近似值是

def sqrt(x):
  y = x*-0.000267
  y = x*(0.004686+y)
  y = x*(-0.034810+y)
  y = x*(0.144780+y)
  y = x*(-0.387893+y)
  y = x*(0.958108+y)
  return y+0.315413

将浮点数规范化，使尾数在[1,4)范围内，然后使用上述算法对其进行处理，最后将指数除以2。不要进行任何浮点除法。

在相同的CPU时间预算下，您可能可以做得更好，但这似乎是一个很好的起点。

- sigfpe

---

如何将尾数归一化，使其范围在 [1,4] 内？请以 x = 612 为例。 - Veridian

如果您正在使用IEEE浮点数，则可以轻松地取得指数。如果它是偶数，则浮点数的形式为a*2^(2n)，其中a在[1,2)。因此，使用上述算法可得到结果sqrt(a)2^n。如果指数是奇数，则我们有a2^(2n+1)=(2a)2^(2n)。这次a在[2,4)中，我们使用上述算法给出：sqrt(2a)*2^n。 - sigfpe