我正在使用Dobb博士的文章「Optimizing Math-Intensive Applications with Fixed-Point Arithmetic」中描述的Anthony Williams的定點庫來計算使用Rhumb Line method的兩個地理點之間的距離。如何提高小數值的定點平方根
當點之間的距離顯着（大於幾公里）時，這種方法的效果不錯，但在較小的距離上很差。最差的情況是當兩點相等或接近相等時，結果是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 ;
}
}
結論 在賈斯汀剝離溶液和分析的光，和對算法在"The Neglected Art of Fixed Point Arithmetic"，我後者改編如下：
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（）。然而一些側面思考提供了更好的解決方案我的應用程序測試距離限制距離的計算距離。事後相當明顯的解決辦法是測試距離的極限平方的距離的平方！
克利福德的輸出 - URL的文章都用水管沖刷（要求一個博聞世網登錄，您很可能登錄和沒有經歷過這一點）。我試圖找到其他地方的文章，但沒有 - 谷歌緩存似乎是最好的選擇。感謝您首先參考文章。 – Dan 2012-01-04 03:36:50
@丹，我去了原來的鏈接，沒有問題。我從來沒有使用UBM，我不認爲我登錄到DDJ，所以我不知道爲什麼有問題。 – 2012-01-04 04:36:51
我已經使用了幾年的代碼;我可能在登錄要求之前下載了該庫。從[安東尼自己的網站]得到它（http://www.justsoftwaresolutions.co.uk/news/optimizing-applications-with-fixed-point-arithmetic.html） – Clifford 2012-01-04 16:06:17