首页 文章

确定浮点平方根

提问于
浏览 1934
6

如何确定浮点数的平方根? 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) . 6121.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.390625b10.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许可的,可能适合或不适合您的需求 .

  • 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时间预算,您可能会做得更好,但这似乎是一个很好的起点 .

  • 2

    最容易实现(您甚至可以在计算器中实现):

    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.
    
  • 1

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

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

相关问题