Cli*_*ord 11 c++ embedded fixed-point sqrt square-root
我正在使用Dobb博士的文章" 使用定点运算优化数学密集型应用 "中描述的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()algorithim简要上面链接的文章第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米的距离是有利的,因此限制不是定点表示.
使用rhumb line方法是该应用的标准主体建议,因此无法更改,并且在任何情况下,在应用程序的其他地方或将来的应用程序中可能需要更准确的平方根函数.
问题:是否有可能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().然而,一些横向思维提供了更好的解决方 我的应用程序测试计算的距离与某个距离限制.后见之明的一个相当明显的解决方案是测试距离的平方与极限的平方!
Kei*_*ith 11
考虑到sqrt(ab) = sqrt(a)sqrt(b),那么你就不能陷在您的数量少,由位给定数目的转移起来的情况下,计算出的根源,转变,回落的比特数的一半得到的结果？
即
sqrt(n) = sqrt(n.2^k)/sqrt(2^k)
= sqrt(n.2^k).2^(-k/2)
例如,对于任何小于2 ^ 8的n,选择k = 28.
原来的实现显然存在一些问题。我对尝试用当前代码完成的方式修复所有这些问题感到沮丧，最终采用了不同的方法。我现在可能可以修复原来的版本，但无论如何我更喜欢我的方式。
我将输入数字视为在 Q64 中开始，这与先移 28，然后再移回 14（开方将其减半）相同。但是，如果您这样做，那么精度将限制为 1/2^14 = 6.1035e-5，因为最后 14 位将为 0。为了解决这个问题，我然后进行移位a和remainder继续填写数字再次循环。代码可以变得更高效、更简洁，但我会将其留给其他人。下面显示的精度几乎与 Q36.28 一样好。如果将输入数字的定点 sqrt 与被定点截断后的浮点 sqrt 进行比较（将其转换为定点并返回），则错误约为 2e-9（我没有在下面的代码，但需要修改一行）。这与 Q36.28 的最佳精度一致，即 1/2^28 = 3.7529e-9。
顺便说一下，原始代码中的一个大错误是从未考虑 m = 0 的项，因此永远无法设置该位。无论如何，这是代码。享受！
#include <iostream>
#include <cmath>
typedef unsigned long uint64_t;
uint64_t sqrt(uint64_t in_val)
{
const uint64_t fixed_resolution_shift = 28;
const unsigned max_shift=62;
uint64_t a_squared=1ULL<<max_shift;
unsigned b_shift=(max_shift>>1) + 1;
uint64_t a=1ULL<<(b_shift - 1);
uint64_t x=in_val;
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=1ULL<<(2*(b_shift - 1));
uint64_t two_a_b=(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((remainder)>=delta && b_shift)
{
a+=(1ULL<<(b_shift - 1));
remainder-=delta;
--b_shift;
}
}
a <<= (fixed_resolution_shift/2);
b_shift = (fixed_resolution_shift/2) + 1;
remainder <<= (fixed_resolution_shift);
while(remainder && b_shift)
{
uint64_t b_squared=1ULL<<(2*(b_shift - 1));
uint64_t two_a_b=(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((remainder)>=delta && b_shift)
{
a+=(1ULL<<(b_shift - 1));
remainder-=delta;
--b_shift;
}
}
return a;
}
double fixed2float(uint64_t x)
{
return static_cast<double>(x) * pow(2.0, -28.0);
}
uint64_t float2fixed(double f)
{
return static_cast<uint64_t>(f * pow(2, 28.0));
}
void finderror(double num)
{
double root1 = fixed2float(sqrt(float2fixed(num)));
double root2 = pow(num, 0.5);
std::cout << "input: " << num << ", fixed sqrt: " << root1 << " " << ", float sqrt: " << root2 << ", finderror: " << root2 - root1 << std::endl;
}
main()
{
finderror(0);
finderror(1e-5);
finderror(2e-5);
finderror(3e-5);
finderror(4e-5);
finderror(5e-5);
finderror(pow(2.0,1));
finderror(1ULL<<35);
}
程序的输出是
input: 0, fixed sqrt: 0 , float sqrt: 0, finderror: 0
input: 1e-05, fixed sqrt: 0.00316207 , float sqrt: 0.00316228, finderror: 2.10277e-07
input: 2e-05, fixed sqrt: 0.00447184 , float sqrt: 0.00447214, finderror: 2.97481e-07
input: 3e-05, fixed sqrt: 0.0054772 , float sqrt: 0.00547723, finderror: 2.43815e-08
input: 4e-05, fixed sqrt: 0.00632443 , float sqrt: 0.00632456, finderror: 1.26255e-07
input: 5e-05, fixed sqrt: 0.00707086 , float sqrt: 0.00707107, finderror: 2.06055e-07
input: 2, fixed sqrt: 1.41421 , float sqrt: 1.41421, finderror: 1.85149e-09
input: 3.43597e+10, fixed sqrt: 185364 , float sqrt: 185364, finderror: 2.24099e-09