Estoy usando la biblioteca de puntos fijos de Anthony Williams descrita en el artículo del Dr. Dobb " Optimización de aplicaciones intensivas en matemáticas con aritmética de punto fijo "para calcular la distancia entre dos puntos geográficos utilizando laétodo @Rhumb Line.

Esto funciona lo suficientemente bien cuando la distancia entre los puntos es significativa (mayor que unos pocos kilómetros), pero es muy pobre a distancias más pequeñas. El peor de los casos es cuando los dos puntos son iguales o casi iguales, el resultado es una distancia de 194 metros, mientras que necesito una precisión de al menos 1 metro a distancias> = 1 metro.

n comparación con una implementación de punto flotante de doble precisión, he localizado el problema en lafixed::sqrt()unción @, que funciona mal en valores pequeños:

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

Corrección del resultado parafixed::sqrt(0) es trivial al tratarlo como un caso especial, pero eso no resolverá el problema para distancias pequeñas que no sean cero, donde el error comienza en 194 metros y converge hacia cero a medida que aumenta la distancia. Probablemente necesito al menos un mejoramiento en el orden de precisión en la precisión hacia cero.

Losfixed::sqrt() Algoritmo se explica brevemente en la página 4 del artículo vinculado anteriormente, pero estoy luchando por seguirlo y mucho menos determinar si es posible mejorarlo. El código para la función se reproduce a continuación:

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

Tenga en cuenta quem_nVal es el valor de representación de punto fijo interno, es unint64_t y la representación usa Q36.28 formato fixed_resolution_shift = 28). La representación en sí tiene suficiente precisión para al menos 8 lugares decimales, y como una fracción del arco ecuatorial es buena para distancias de alrededor de 0,14 metros, por lo que la limitación no es la representación de punto fij

El uso del método de línea de rumbo es una recomendación del cuerpo de estándares para esta aplicación, por lo que no se puede cambiar, y en cualquier caso es probable que se requiera una función de raíz cuadrada más precisa en otra parte de la aplicación o en futuras aplicaciones.

Pregunta ¿Es posible mejorar la precisión de lafixed::sqrt() algoritmo para valores pequeños que no son cero mientras se mantiene su convergencia limitada y determinista?

Información Adiciona El código de prueba utilizado para generar la tabla anterior:

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

Conclusió A la luz de la solución y el análisis de Justin Peel, y la comparación con el algoritmo en "El arte descuidado de la aritmética de punto fijo", He adaptado este último de la siguiente manera:

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

Aunque esto proporciona una precisión mucho mayor, la mejora que necesitaba no se logra. El formato Q36.28 solo proporciona la precisión que necesito, pero no es posible realizar un sqrt () sin perder algunos bits de precisión. Sin embargo, algo de pensamiento lateral proporciona una mejor solución. Mi aplicación prueba la distancia calculada contra algún límite de distancia. ¡La solución bastante obvia en retrospectiva es probar el cuadrado de la distancia contra el cuadrado del límite!