Как улучшить извлечение квадратного корня из фиксированной точки для малых значений

Я использую библиотеку с фиксированной точкой Энтони Вильямса, описанную в статье доктора Добба « Оптимизация математических приложений с помощью арифметики с фиксированной точкой », чтобы вычислить расстояние между двумя географическими точками с помощью Rhumb Линейный метод .

Это работает достаточно хорошо, когда расстояние между точками значительно (больше нескольких километров), но очень мало на меньших расстояниях. В худшем случае, когда две точки равны или почти равны, в результате получается расстояние 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<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);

        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<

Обратите внимание, что m_nVal - это внутреннее значение представления с фиксированной точкой, это int64_t , а представление использует Q36.28 формат ( fixed_resolution_shift = 28). Само представление имеет достаточную точность как минимум для 8 десятичных знаков, а доля экваториальной дуги подходит для расстояний около 0,14 метра, поэтому ограничение не относится к представлению с фиксированной точкой.

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

Вопрос: Можно ли повысить точность алгоритма fixed :: sqrt () для небольших ненулевых значений, сохранив при этом его ограниченную и детерминированную сходимость?

Дополнительная информация Тестовый код, использованный для создания приведенной выше таблицы:

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