我正在使用 Dobb 博士的文章“ Optimizing Math-Intensive Applications with Fixed-Point Arithmetic ”中描述的 Anthony Williams 的定点库，使用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()。然而，一些横向思维提供了更好的解决方案。我的应用程序根据某个距离限制测试计算出的距离。事后看来，相当明显的解决方案是测试距离的平方与极限的平方！