Я использую библиотеку Энтони Уильямса с фиксированной запятой, описанную в статье доктора Добба "Оптимизация математических приложений с арифметикой с фиксированной точкой«рассчитать расстояние между двумя географическими точками, используяМетод Rhumb Line.

Это работает достаточно хорошо, когда расстояние между точками значительно (больше нескольких километров), но очень мало на меньших расстояниях. В худшем случае, когда две точки равны или почти равны, результатом является расстояние 194 метра, в то время как мне нужна точность не менее 1 метра на расстояниях> = 1 метр.

По сравнению с реализацией с плавающей запятой двойной точности, я нашел проблему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 метра, поэтому ограничение не является представлением с фиксированной точкой.

Использование метода прямой линии является рекомендацией органа стандартизации для этого приложения, поэтому его нельзя изменить, и в любом случае более точная функция квадратного корня, вероятно, потребуется в другом месте приложения или в будущих приложениях.

Вопрос: Можно ли улучшить точность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 ;
    }
}

Заключение В свете решения и анализа Джастина Пила, а также сравнения с алгоритмом в«Заброшенное искусство арифметики с фиксированной точкой»Я адаптировал последнее следующим образом:

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 () без потери точности в несколько бит. Однако некоторое латеральное мышление дает лучшее решение. Мое приложение проверяет рассчитанное расстояние с некоторым пределом расстояния. Довольно очевидное решение в ретроспективе - проверить квадрат расстояния против квадрата предела!