Я использую библиотеку фиксированной точки Энтони Уильямса, описанную в статье доктора Добба "Оптимизация приложений, интенсивно использующих математические вычисления, с помощью арифметики с фиксированной точкой" для расчета расстояния между двумя географическими точками с использованием метод локсодромии.
Это работает достаточно хорошо, когда расстояния между точками значительны (более нескольких километров), но очень плохо на меньших расстояниях. В худшем случае, когда две точки равны или почти равны, результатом будет расстояние 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() без потери нескольких битов точности. Однако некоторое нестандартное мышление обеспечивает лучшее решение. Мое приложение проверяет расчетное расстояние на соответствие некоторому пределу расстояния. Оглядываясь назад, довольно очевидное решение состоит в том, чтобы сравнить квадрат расстояния с квадратом предела!
fixed::sqrt(), показанные в таблице? Какой компилятор + ОС вы использовали? Я не получаю одинаковые числа ни для чего, кроме квадратного корня из 0. Ни с gcc (DJGPP/DOS), ни с Open Watcom (Windows). Все мои результаты отличаются примерно от 10^-5 до 10^-6 дляxиз таблицы, а не на 10^-7 или 10^-9. Вы использовали более 28 дробных бит при заполнении таблицы? Как вы конвертировали в/изfixedи каков размер вашего типа с плавающей запятой (кстати, это было double или long double)? - person Alexey Frunze schedule 05.01.2012fabs( std::sqrt(x) - sqrt(fixed(x)))- person Clifford schedule 05.01.2012fixed::as_double().as_double()выполняетx / (double)1<<fixed_resolution_shift- person Clifford schedule 05.01.2012