我正在使用Anthony Williams的固定点库，该库在Dr Dobb's文章“使用固定点算术优化数学密集型应用程序”中进行了描述，使用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米的距离，因此限制不在于固定点表示。
在此应用中，使用航向线方法是标准机构的建议，因此无法更改，而且在任何情况下，可能需要在应用程序的其他地方或未来的应用程序中使用更准确的平方根函数。
问题：是否可以提高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 ;
}
}
结论 根据Justin Peel的解决方案和分析，以及与“固定点算术的被忽视艺术”中的算法进行比较，我已将后者改编为以下形式：
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)，我的结果都与表格中的x相差约10^-5到10^-6。当您填写表格时，您是否使用了超过28位的小数位？您是如何转换成/从fixed的？您的浮点类型大小是多少（顺便问一下，它是double还是long double）？ - Alexey Frunze