How to improve fixed point square-root for small values

Question

I am using Anthony Williams' fixed point library described in the Dr Dobb's article "Optimizing Math-Intensive Applications with Fixed-Point Arithmetic" to calculate the distance between two geographical points using the Rhumb Line method.

This works well enough when the distance between the points is significant (greater than a few kilometers), but is very poor at smaller distances. The worst case being when the two points are equal or near equal, the result is a distance of 194 meters, while I need precision of at least 1 metre at distances >= 1 metre.

By comparison with a double precision floating-point implementation, I have located the problem to the fixed::sqrt() function, which performs poorly at small values:

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

Correcting the result for fixed::sqrt(0) is trivial by treating it as a special case, but that will not solve the problem for small non-zero distances, where the error starts at 194 metres and converges toward zero with increasing distance. I probably need at least an order of maginitude improvement in precision toward zero.

The fixed::sqrt() algorithim is briefly explained on page 4 of the article linked above, but I am struggling to follow it let alone determine whether it is possible to improve it. The code for the function is reproduced below:

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);
}

Note that m_nVal is the internal fixed point representation value, it is an int64_t and the representation uses Q36.28 format (fixed_resolution_shift = 28). The representation itself has enough precision for at least 8 decimal places, and as a fraction of equatorial arc is good for distances of around 0.14 metres, so the limitation is not the fixed-point representation.

Use of the rhumb line method is a standards body recommendation for this application so cannot be changed, and in any case a more accurate square-root function is likely to be required elsewhere in the application or in future applications.

Question: Is it possible to improve the accuracy of the fixed::sqrt() algorithm for small non-zero values while still maintaining its bounded and deterministic convergence?

Additional Information The test code used to generate the table above:

#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 ;
    }
}

Conclusion In the light of Justin Peel's solution and analysis, and comparison with the algorithm in "The Neglected Art of Fixed Point Arithmetic", I have adapted the latter as follows:

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);
}

While this gives far greater precision, the improvement I needed is not to be achieved. The Q36.28 format alone just about provides the precision I need, but it is not possible to perform a sqrt() without loss of a few bits of precision. However some lateral thinking provides a better solution. My application tests the calculated distance against some distance limit. The rather obvious solution in hindsight is to test the square of the distance against the square of the limit!

Answer 1

Given that sqrt(ab) = sqrt(a)sqrt(b), then can't you just trap the case where your number is small and shift it up by a given number of bits, compute the root and shift that back down by half the number of bits to get the result?

I.e.

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

E.g. Choose k = 28 for any n less than 2^8.

Answer 2

The original implementation obviously has some problems. I became frustrated with trying to fix them all with the way the code is currently done and ended up going at it with a different approach. I could probably fix the original now, but I like my way better anyway.

I treat the input number as being in Q64 to start which is the same as shifting by 28 and then shifting back by 14 afterwards (the sqrt halves it). However, if you just do that, then the accuracy is limited to 1/2^14 = 6.1035e-5 because the last 14 bits will be 0. To remedy this, I then shift a and remainder correctly and to keep filling in digits I do the loop again. The code can be made more efficient and cleaner, but I'll leave that to someone else. The accuracy shown below is pretty much as good as you can get with Q36.28. If you compare the fixed point sqrt with the floating point sqrt of the input number after it has been truncated by fixed point(convert it to fixed point and back), then the errors are around 2e-9(I didn't do this in the code below, but it requires one line of change). This is right in line with the best accuracy for Q36.28 which is 1/2^28 = 3.7529e-9.

By the way, one big mistake in the original code is that the term where m = 0 is never considered so that bit can never be set. Anyway, here is the code. Enjoy!

#include <iostream>
#include <cmath>

typedef unsigned long uint64_t;

uint64_t sqrt(uint64_t in_val)
{
    const uint64_t fixed_resolution_shift = 28;
    const unsigned max_shift=62;
    uint64_t a_squared=1ULL<<max_shift;
    unsigned b_shift=(max_shift>>1) + 1;
    uint64_t a=1ULL<<(b_shift - 1);

    uint64_t x=in_val;

    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=1ULL<<(2*(b_shift - 1));
        uint64_t two_a_b=(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((remainder)>=delta && b_shift)
        {
            a+=(1ULL<<(b_shift - 1));
            remainder-=delta;
            --b_shift;
        }
    }
    a <<= (fixed_resolution_shift/2);
    b_shift = (fixed_resolution_shift/2) + 1;
    remainder <<= (fixed_resolution_shift);

    while(remainder && b_shift)
    {
        uint64_t b_squared=1ULL<<(2*(b_shift - 1));
        uint64_t two_a_b=(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((remainder)>=delta && b_shift)
        {
            a+=(1ULL<<(b_shift - 1));
            remainder-=delta;
            --b_shift;
        }
    }

    return a;
}

double fixed2float(uint64_t x)
{
    return static_cast<double>(x) * pow(2.0, -28.0);
}

uint64_t float2fixed(double f)
{
    return static_cast<uint64_t>(f * pow(2, 28.0));
}

void finderror(double num)
{
    double root1 = fixed2float(sqrt(float2fixed(num)));
    double root2 = pow(num, 0.5);
    std::cout << "input: " << num << ", fixed sqrt: " << root1 << " " << ", float sqrt: " << root2 << ", finderror: " << root2 - root1 << std::endl;
}

main()
{
    finderror(0);
    finderror(1e-5);
    finderror(2e-5);
    finderror(3e-5);
    finderror(4e-5);
    finderror(5e-5);
    finderror(pow(2.0,1));
    finderror(1ULL<<35);
}

with the output of the program being

input: 0, fixed sqrt: 0 , float sqrt: 0, finderror: 0
input: 1e-05, fixed sqrt: 0.00316207 , float sqrt: 0.00316228, finderror: 2.10277e-07
input: 2e-05, fixed sqrt: 0.00447184 , float sqrt: 0.00447214, finderror: 2.97481e-07
input: 3e-05, fixed sqrt: 0.0054772 , float sqrt: 0.00547723, finderror: 2.43815e-08
input: 4e-05, fixed sqrt: 0.00632443 , float sqrt: 0.00632456, finderror: 1.26255e-07
input: 5e-05, fixed sqrt: 0.00707086 , float sqrt: 0.00707107, finderror: 2.06055e-07
input: 2, fixed sqrt: 1.41421 , float sqrt: 1.41421, finderror: 1.85149e-09
input: 3.43597e+10, fixed sqrt: 185364 , float sqrt: 185364, finderror: 2.24099e-09