确定浮点平方根-Java 学习之路

如何确定浮点数的平方根？ Newton-Raphson方法是一种好方法吗？我也没有硬件平方根 . 我也没有硬件鸿沟（但我实现了浮点除法） .

如果可能的话，我宁愿尽可能减少分歧的数量，因为它们太贵了 .

另外，应该是减少迭代总数的初始猜测???

非常感谢！

4 回答

10
当你使用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)
```
请注意（与平方根的迭代不同），倒数平方根的迭代不涉及任何划分，因此通常更有效 .

我在关于鸿沟的问题中提到过，你应该看看现有的软浮动库，而不是重新发明轮子 . 这个建议也适用于此 . 此功能已在现有的软浮动库中实现 .

编辑：提问者似乎仍然感到困惑，所以让我们举个例子： 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 . 现在做一个newton-raphson迭代来找到 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是源和目标类型的精度） . 但是，您也可以采用计算sqrt（f）的策略略大于p位，将该值平方，并在必要时将尾随位调整为1（通常更便宜） .

sqrt很不错，因为它是一个一元函数，它可以对商用硬件上的单精度进行详尽的测试 .

你可以在opensource.apple.com上找到OS X soft-float sqrtf 函数，该函数使用上面描述的算法（我写的，它发生了） . 它是根据APSL许可的，可能适合或不适合您的需求 .
回复于 2024-04-20T10:07:35+08:00
1

可能（仍然）是查找inverse square root和我最喜欢的10行代码的最快实现 .

它's based on Newton Approximation, but with a few quirks. There'甚至是great story左右 .

回复于 2024-04-20T10:07:35+08:00

最容易实现（您甚至可以在计算器中实现）：

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

这与newton raphson完全相同：

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的位置） .

假设我们已经知道了^ 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.

回复于 2024-04-20T10:07:35+08:00

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时间预算，您可能会做得更好，但这似乎是一个很好的起点 .
回复于 2024-04-20T10:07:35+08:00

确定浮点平方根

4 回答

相关问题