Utilizzo la libreria a punti fissi di Anthony Williams descritta nell'articolo del dott. Dobb "Optimizing Math-Intensive Applications with Fixed-Point Arithmetic" per calcolare la distanza tra due punti geografici utilizzando Rhumb Line method.Come migliorare la radice quadrata a punti fissi per valori piccoli
Questo funziona abbastanza bene quando la distanza tra i punti è significativa (superiore a pochi chilometri), ma è molto scarsa a distanze minori. Nel peggiore dei casi, quando i due punti sono uguali o quasi uguali, il risultato è una distanza di 194 metri, mentre ho bisogno di una precisione di almeno 1 metro a distanze> = 1 metro.
In confronto con una doppia precisione implementazione in virgola mobile, ho individuato il problema alla funzione fixed::sqrt(), che funziona male a piccoli valori:
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
Correzione risultato fixed::sqrt(0) è banale ravvisandovi un caso speciale, ma che non risolverà il problema per piccole distanze diverse da zero, dove l'errore inizia a 194 metri e converge verso lo zero con l'aumentare della distanza. Probabilmente ho bisogno di almeno un ordine di miglioramento della magnitudine in precisione verso lo zero.
L'algoritmo fixed::sqrt() è brevemente spiegato a pagina 4 dell'articolo correlato sopra, ma sto lottando per seguirlo e tanto meno determinare se sia possibile migliorarlo. Il codice per la funzione è riportato di seguito:
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);
}
noti che m_nVal è il valore rappresentazione punto fisso interno, è un int64_t e la rappresentazione utilizza Q36.28 formato (fixed_resolution_shift = 28). La rappresentazione stessa ha una precisione sufficiente per almeno 8 posizioni decimali e una frazione dell'arco equatoriale è buona per distanze di circa 0,14 metri, quindi la limitazione non è la rappresentazione in virgola fissa.
L'uso del metodo della lossodromia è un consiglio per il corpo standard per questa applicazione, pertanto non può essere modificato e, in ogni caso, è probabile che sia necessaria una funzione radice quadrata più accurata altrove nell'applicazione o nelle applicazioni future.
Domanda: È possibile migliorare la precisione dell'algoritmo fixed::sqrt() per piccoli valori diversi da zero pur mantenendo la convergenza limitata e deterministica?
Informazioni aggiuntive Il codice di test utilizzato per generare la tabella di cui sopra:
#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 ;
}
}
Conclusione Alla luce della soluzione ed analisi di Justin Peel, e confronto con l'algoritmo in "The Neglected Art of Fixed Point Arithmetic", ho adattato quest'ultimo come segue:
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);
}
Mentre questo dà maggiore precisione, il miglioramento di cui avevo bisogno non è stato raggiunto. Il formato Q36.28 fornisce solo la precisione di cui ho bisogno, ma non è possibile eseguire un sqrt() senza perdita di alcuni bit di precisione. Tuttavia alcuni pensieri laterali offrono una soluzione migliore. La mia applicazione verifica la distanza calcolata rispetto ad un limite di distanza. La soluzione piuttosto ovvia col senno di poi è testare il quadrato della distanza contro il quadrato del limite!
Clifford - l'URL per l'articolo è stato hosed (per chiedere un accesso UBM TechWeb, si erano probabilmente effettuato l'accesso e non ha esperienza questo). Ho cercato di trovare l'articolo altrove ma non è stato così: la cache di Google sembrava essere l'alternativa migliore. Grazie per aver fatto riferimento all'articolo in primo luogo. – Dan
@Dan, sono andato al link originale e non ho avuto problemi. Non ho mai usato UBM e non penso di aver effettuato l'accesso a DDJ, quindi non so perché c'è un problema. –
Ho usato il codice per un paio di anni; Ho scaricato la libreria probabilmente prima del requisito di accesso. Scaricalo dal [sito web di Anthony] (http://www.justsoftwaresolutions.co.uk/news/optimizing-applications-with-fixed-point-arithmetic.html) – Clifford