Dobb博士の記事「固定小数点演算による数学集約型アプリケーションの最適化」で説明されているAnthonyWilliamsの固定小数点ライブラリを使用して、 RhumbLineメソッドを使用して2つの地理的ポイント間の距離を計算しています。
これは、ポイント間の距離が大きい場合（数キロメートルを超える場合）は十分に機能しますが、距離が短い場合は非常に不十分です。最悪の場合、2つのポイントが等しいかほぼ等しい場合、結果は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メートルから始まり、距離が長くなるにつれてゼロに向かって収束する、ゼロ以外の小さな距離の問題は解決されません。おそらく、ゼロに向けて精度を少なくとも1桁改善する必要があります。
アルゴリズムは上記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（）を実行することはできません。ただし、いくつかの水平思考はより良い解決策を提供します。私のアプリケーションは、計算された距離をある距離制限に対してテストします。後から考えると、かなり明白な解決策は、距離の2乗を限界の2乗に対してテストすることです。