Я использую библиотеку Энтони Уильямса с фиксированными точками, описанную в статье доктора Добба " Оптимизация математических приложений с арифметикой с фиксированной точкой ", чтобы вычислить расстояние между двумя географическими точкамииспользование метода 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 метра, поэтому ограничение не является представлением с фиксированной точкой.
Использование прямого удараМетод line является рекомендацией органа по стандартизации для этого приложения, поэтому его нельзя изменить, и в любом случае более точная функция квадратного корня может потребоваться в другом месте приложения или в будущих приложениях.
Вопрос: Можно ли повысить точность алгоритма 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 () без потери точности в несколько бит.Однако некоторое латеральное мышление дает лучшее решение.Мое приложение проверяет рассчитанное расстояние с некоторым пределом расстояния.Довольно очевидное решение в ретроспективе - проверить квадрат расстояния против квадрата предела!