如何提高小值的定点平方根 [英] How to improve fixed point square-root for small values
问题描述
我使用Anthony Williams的Dobb博士文章中描述的定点库。 使用定点算法优化数学密集型应用 使用 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 amp; 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&& amp; amp;<（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 * remaining）> 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方法是此应用程序的标准机构建议，因此不能更改，在任何情况下，在应用程序或未来应用程序中可能需要更准确的平方根函数。 / p>
问题：可以提高 fixed :: sqrt（）算法对小的非零值，同时仍保持其有界和确定性收敛？
其他信息
测试代码用于生成上面的表：
#include< cmath>
#include< iostream>
#includefixed.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'<<错误< std :: endl;
}
}
结论
根据Justin Peel的解决方案和分析，并与定点算术的被忽视的艺术，我修改了后者如下：
fixed fixed :: sqrt（）const
{
uint64_t a = 0; // root accumulator
uint64_t remHi = 0; //高部分的部分剩余
uint64_t remLo = m_nVal; // low part of partial remainder
uint64_t testDiv;
int count = 31 +（fixed_resolution_shift>> 1）; //循环计数器
do
{
//获得2位arg
remHi =（remHi << 2）| （remLo>> 62）; remLo <= 2;
//为根中的下一位准备好
a<< = 1;
//测试激进
testDiv =（a<< 1）+ 1;
if（remHi> = testDiv）
{
remHi - = testDiv;
a + = 1;
}
} while（count--！= 0）;
return fixed（internal（），a）;
}
虽然这给出了更高的精度，但我需要的改进。单独的Q36.28格式提供了我需要的精度，但是不可能执行sqrt（）而不损失几个位的精度。然而，一些横向思维提供了更好的解决方案。我的应用程序测试计算距离与一些距离限制。事后看起来相当明显的解决方案是测试距离的平方与极限的平方！
有一些问题。我试图修复他们所有的代码当前的方式，结束了以一种不同的方法，变得沮丧。
我将输入的数字视为在Q64中开始，这与移动28相同。然后再向后移动14（sqrt减半）。然而，如果你只是这样做，那么准确性被限制为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 amp; a_squared> x）
{
a>> = 1;
a_squared>> = 2;
--b_shift;
}
uint64_t remainder = x-a_squared;
--b_shift;
while（remaining& b_shift）
{
uint64_t b_squared = 1ULL<<（2 *（b_shift - 1））;
uint64_t two_a_b =（a<< b_shift）;
while（b_shift&& amp; amp;<（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（remaining& b_shift）
{
uint64_t b_squared = 1ULL<<（2 *（b_shift - 1））;
uint64_t two_a_b =（a<< b_shift）;
while（b_shift&& amp; amp;<（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）;
}
，程序的输出为
输入：0，固定sqrt：0，float sqrt：0，finderror：0
输入：1e-05，fixed sqrt：0.00316207，float sqrt： 0.00316228，finderror：2.10277e-07
输入：2e-05，固定sqrt：0.00447184，float sqrt：0.00447214，finderror：2.97481e-07
输入：3e-05，fixed sqrt：0.0054772，float sqrt：0.00547723，finderror：2.43815e-08
输入：4e-05，固定sqrt：0.00632443，float sqrt：0.00632456，finderror：1.26255e-07
输入：5e-05，fixed sqrt：0.00707086 ，float sqrt：0.00707107，finderror：2.06055e-07
输入：2，固定sqrt：1.41421，float sqrt：1.41421，finderror：1.85149e-09
输入：3.43597e + 10， 185364，float sqrt：185364，finderror：2.24099e-09
I am using Anthony Williams' fixed point library described in the Dr Dobb's article "Optimizing Math-Intensive Applications with Fixed-Point Arithmetic" to calculate the distance between two geographical points using the Rhumb Line method.
This works well enough when the distance between the points is significant (greater than a few kilometers), but is very poor at smaller distances. The worst case being when the two points are equal or near equal, the result is a distance of 194 meters, while I need precision of at least 1 metre at distances >= 1 metre.
By comparison with a double precision floating-point implementation, I have located the problem to the fixed::sqrt() function, which performs poorly at small values:
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
Correcting the result for fixed::sqrt(0) is trivial by treating it as a special case, but that will not solve the problem for small non-zero distances, where the error starts at 194 metres and converges toward zero with increasing distance. I probably need at least an order of maginitude improvement in precision toward zero.
The fixed::sqrt() algorithim is briefly explained on page 4 of the article linked above, but I am struggling to follow it let alone determine whether it is possible to improve it. The code for the function is reproduced below:
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);
}
Note that m_nVal is the internal fixed point representation value, it is an int64_t and the representation uses Q36.28 format (fixed_resolution_shift = 28). The representation itself has enough precision for at least 8 decimal places, and as a fraction of equatorial arc is good for distances of around 0.14 metres, so the limitation is not the fixed-point representation.
Use of the rhumb line method is a standards body recommendation for this application so cannot be changed, and in any case a more accurate square-root function is likely to be required elsewhere in the application or in future applications.
Question: Is it possible to improve the accuracy of the fixed::sqrt() algorithm for small non-zero values while still maintaining its bounded and deterministic convergence?
Additional Information The test code used to generate the table above:
#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 ;
}
}
Conclusion In the light of Justin Peel's solution and analysis, and comparison with the algorithm in "The Neglected Art of Fixed Point Arithmetic", I have adapted the latter as follows:
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);
}
While this gives far greater precision, the improvement I needed is not to be achieved. The Q36.28 format alone just about provides the precision I need, but it is not possible to perform a sqrt() without loss of a few bits of precision. However some lateral thinking provides a better solution. My application tests the calculated distance against some distance limit. The rather obvious solution in hindsight is to test the square of the distance against the square of the limit!
The original implementation obviously has some problems. I became frustrated with trying to fix them all with the way the code is currently done and ended up going at it with a different approach. I could probably fix the original now, but I like my way better anyway.
I treat the input number as being in Q64 to start which is the same as shifting by 28 and then shifting back by 14 afterwards (the sqrt halves it). However, if you just do that, then the accuracy is limited to 1/2^14 = 6.1035e-5 because the last 14 bits will be 0. To remedy this, I then shift a and remainder correctly and to keep filling in digits I do the loop again. The code can be made more efficient and cleaner, but I'll leave that to someone else. The accuracy shown below is pretty much as good as you can get with Q36.28. If you compare the fixed point sqrt with the floating point sqrt of the input number after it has been truncated by fixed point(convert it to fixed point and back), then the errors are around 2e-9(I didn't do this in the code below, but it requires one line of change). This is right in line with the best accuracy for Q36.28 which is 1/2^28 = 3.7529e-9.
By the way, one big mistake in the original code is that the term where m = 0 is never considered so that bit can never be set. Anyway, here is the code. Enjoy!
#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);
}
with the output of the program being
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
这篇关于如何提高小值的定点平方根的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
