c++ - 如何提高小值的定点平方根

Question

我正在使用 Dr Dobb 的文章“ Optimizing Math-Intensive Applications with Fixed-Point Arithmetic ”中描述的 Anthony Williams 的定点库来使用 Rhumb Line method 计算两个地理点之间的距离。

当点之间的距离很大(大于几公里)时，此方法效果很好，但在较小的距离时效果很差。最坏的情况是当两点相等或接近相等时，结果是 194 米的距离，而我需要在距离 >= 1 米时至少有 1 米的精度。

通过与 double 浮点实现的比较，我将问题定位到 fixed::sqrt() 函数，该函数在小值时表现不佳:

x       std::sqrt(x)    fixed::sqrt(x)  error
----------------------------------------------------
0       0               3.05176e-005    3.05176e-005
1e-005  0.00316228      0.00316334      1.06005e-006
2e-005  0.00447214      0.00447226      1.19752e-007
3e-005  0.00547723      0.0054779       6.72248e-007
4e-005  0.00632456      0.00632477      2.12746e-007
5e-005  0.00707107      0.0070715       4.27244e-007
6e-005  0.00774597      0.0077467       7.2978e-007
7e-005  0.0083666       0.00836658      1.54875e-008
8e-005  0.00894427      0.00894427      1.085e-009

将 fixed::sqrt(0) 的结果作为特例来修正是微不足道的，但这并不能解决小的非零距离的问题，错误开始于194 米，随着距离的增加向零收敛。我可能至少需要将精度提高一个数量级才能达到零。

fixed::sqrt() 算法在上面链接的文章的第 4 页有简要说明，但我很难遵循它，更不用说确定是否可以改进它了。该函数的代码转载如下:

fixed fixed::sqrt() const
{
    unsigned const max_shift=62;
    uint64_t a_squared=1LL<<max_shift;
    unsigned b_shift=(max_shift+fixed_resolution_shift)/2;
    uint64_t a=1LL<<b_shift;

    uint64_t x=m_nVal;

    while(b_shift && a_squared>x)
    {
        a>>=1;
        a_squared>>=2;
        --b_shift;
    }

    uint64_t remainder=x-a_squared;
    --b_shift;

    while(remainder && b_shift)
    {
        uint64_t b_squared=1LL<<(2*b_shift-fixed_resolution_shift);
        int const two_a_b_shift=b_shift+1-fixed_resolution_shift;
        uint64_t two_a_b=(two_a_b_shift>0)?(a<<two_a_b_shift):(a>>-two_a_b_shift);

        while(b_shift && remainder<(b_squared+two_a_b))
        {
            b_squared>>=2;
            two_a_b>>=1;
            --b_shift;
        }
        uint64_t const delta=b_squared+two_a_b;
        if((2*remainder)>delta)
        {
            a+=(1LL<<b_shift);
            remainder-=delta;
            if(b_shift)
            {
                --b_shift;
            }
        }
    }
    return fixed(internal(),a);
}

注意m_nVal是内部定点表示值，它是一个int64_t，表示使用Q36.28格式(fixed_resolution_shift = 28) .表示本身具有至少 8 位小数的足够精度，并且作为赤道弧的一小部分适用于大约 0.14 米的距离，因此限制不是定点表示。

使用恒向线法是标准机构针对此应用程序的建议，因此无法更改，无论如何，在应用程序的其他地方或 future 的应用程序中可能需要更精确的平方根函数。

问题:是否有可能提高 fixed::sqrt() 算法对小的非零值的准确性，同时仍然保持其有界和确定性收敛？

附加信息用于生成上表的测试代码:

#include <cmath>
#include <iostream>
#include "fixed.hpp"

int main()
{
    double error = 1.0 ;
    for( double x = 0.0; error > 1e-8; x += 1e-5 )
    {
        double fixed_root = sqrt(fixed(x)).as_double() ;
        double std_root = std::sqrt(x) ;
        error = std::fabs(fixed_root - std_root) ;
        std::cout << x << '\t' << std_root << '\t' << fixed_root << '\t' << error << std::endl ;
    }
}

结论根据Justin Peel的解决方案和分析，并与"The Neglected Art of Fixed Point Arithmetic"中的算法进行比较，我对后者进行了如下调整:

fixed fixed::sqrt() const
{
    uint64_t a = 0 ;            // root accumulator
    uint64_t remHi = 0 ;        // high part of partial remainder
    uint64_t remLo = m_nVal ;   // low part of partial remainder
    uint64_t testDiv ;
    int count = 31 + (fixed_resolution_shift >> 1); // Loop counter
    do 
    {
        // get 2 bits of arg
        remHi = (remHi << 2) | (remLo >> 62); remLo <<= 2 ;

        // Get ready for the next bit in the root
        a <<= 1;   

        // Test radical
        testDiv = (a << 1) + 1;    
        if (remHi >= testDiv) 
        {
            remHi -= testDiv;
            a += 1;
        }

    } while (count-- != 0);

    return fixed(internal(),a);
}

虽然这提供了更高的精度，但我需要的改进并没有实现。仅 Q36.28 格式就可以提供我需要的精度，但不可能在不损失几位精度的情况下执行 sqrt()。然而，一些横向思维提供了更好的解决方案。我的应用程序根据某个距离限制测试计算出的距离。事后看来，相当明显的解决方案是测试距离的平方与极限的平方!

Answer 1

鉴于 sqrt(ab) = sqrt(a)sqrt(b)，那么你不能只捕获你的数字较小的情况并将其向上移动给定的位数, 计算根并将其向下移动一半的位数以获得结果？

即

 sqrt(n) = sqrt(n.2^k)/sqrt(2^k)
         = sqrt(n.2^k).2^(-k/2)

例如对于小于 2^8 的任何 n，选择 k = 28。