如何解决平方根计算算法
我一直在用 C 语言实现控制软件,其中一种控制算法需要平方根计算。我一直在寻找合适的平方根计算算法,无论基数如何,它的执行时间都是恒定的。此要求排除了标准库中的 sqrt
函数。
就我的平台而言,我一直在使用基于浮点 32 位 ARM Cortex A9 的机器。至于我的应用程序中的基数范围,算法是以物理单位计算的,所以我希望遵循范围 <0,400>
。至于所需的误差,我认为大约 1% 的误差就足够了。任何人都可以向我推荐适合我的目的的平方根计算算法吗?
解决方法
我最初的方法是使用泰勒级数求平方根,并在多个固定点处预先计算系数。这会将计算减少为减法和乘法。
查找表将是一个二维数组,如:
point | C0 | C1 | C2 | C3 | C4 | ...
-----------------------------------------
0.5 | f00 | f01 | f02 | f03 | f04 |
-----------------------------------------
1.0 | f10 | f11 | f12 | f13 | f14 |
-----------------------------------------
1.5 | f20 | f21 | f22 | f23 | f24 |
-----------------------------------------
....
所以在计算 sqrt(x) 时,使用最接近 x 的点的表格行。
示例:
sqrt(1.1) (i.e. use point 1.0 coeffients)
f10 +
f11 * (1.1 - 1.0) +
f12 * (1.1 - 1.0) ^ 2 +
f13 * (1.1 - 1.0) ^ 3 +
f14 * (1.1 - 1.0) ^ 4
上表建议您预先计算系数的点之间的固定距离(即每个点之间的 0.5)。但是,由于平方根的性质,您可能会发现点之间的距离对于不同的 x
范围会有所不同。例如 x
在 [0 - 1] -> 距离 0.1,x
在 [1 - 2] -> 距离 0.25,x
在 [2 - 10] -> 距离 0.5 和以此类推。
另一件事是获得所需精度所需的项数。在这里您可能还会发现 x
的不同范围可能需要不同数量的系数。
所有这些都可以在普通计算机上轻松预先计算(例如使用 excel)。
注意:对于非常接近于零的值,这种方法不好。也许牛顿法会是更好的选择。
泰勒系列:https://en.wikipedia.org/wiki/Taylor_series
牛顿法:https://en.wikipedia.org/wiki/Newton%27s_method
也相关:https://math.stackexchange.com/questions/291168/algorithms-for-approximating-sqrt2
,Arm v7 指令集为两个同时逼近的倒数平方根计算 vrsqrte_f32
和四个逼近的 vrsqrteq_f32
提供了快速指令。 (标量变体 vrsqrtes_f32
仅在 Arm64 v8.2 上可用)。
那么结果可以简单地通过x * vrsqrte_f32(x);
计算,它在整个正值x范围内具有优于0.33%的相对准确度。见https://www.mdpi.com/2079-3197/9/2/21/pdf
ARM NEON 指令 FRSQRTE 给出了 8.25 位正确的结果。
在 x==0
vrsqrtes_f32(x) == Inf,所以 x*vrsqrtes_f32(x) 将是 NaN。
如果x==0的值不可避免,那么最优的两条指令序列需要稍微调整一下:
float sqrtest(float a) {
// need to "transfer" or "convert" the scalar input
// to a vector of two
// - optimally we would not need an instruction for that
// but we would just let the processor calculate the instruction
// for all the lanes in the register
float32x2_t a2 = vdup_n_f32(a);
// next we create a mask that is all ones for the legal
// domain of 1/sqrt(x)
auto is_legal = vreinterpret_f32_u32(vcgt_f32(a2,vdup_n_f32(0.0f)));
// calculate two reciprocal estimates in parallel
float32x2_t a2est = vrsqrte_f32(a2);
// we need to mask the result,so that effectively
// all non-legal values of a2est are zeroed
a2est = vand_u32(is_legal,a2est);
// x * 1/sqrt(x) == sqrt(x)
a2 = vmul_f32(a2,a2est);
// finally we get only the zero lane of the result
// discarding the other half
return vget_lane_f32(a2,0);
}
当然,这种方法的吞吐量几乎是
的两倍void sqrtest2(float &a,float &b) {
float32x2_t a2 = vset_lane_f32(b,vdup_n_f32(a),1);
float32x2_t is_legal = vreinterpret_f32_u32(vcgt_f32(a2,vdup_n_f32(0.0f)));
float32x2_t a2est = vrsqrte_f32(a2);
a2est = vand_u32(is_legal,a2est);
a2 = vmul_f32(a2,a2est);
a = vget_lane_f32(a2,0);
b = vget_lane_f32(a2,1);
}
如果您可以直接使用 float32x2_t
或 float32x4_t
输入和输出,那就更好了。
float32x2_t sqrtest2(float32x2_t a2) {
float32x2_t is_legal = vreinterpret_f32_u32(vcgt_f32(a2,a2est);
return vmul_f32(a2,a2est);
}
这个实现给出了 sqrtest2(1) == 0.998
和 sqrtest2(400) == 19.97
(在带有 arm64 的 MacBook M1 上测试)。由于无分支和 LUT 自由,假设所有指令在恒定数量的周期内执行,这可能具有恒定的执行时间。
我决定使用以下方法。我选择了牛顿法,然后我通过实验设置了固定的迭代次数,以便在整个基数范围内的误差,即 public void sortByID()
{
DefaultListModel tmp = new DefaultListModel();
for(int x = 0; x < listModel.size(); x++)
{
String[] a = listModel.get(x).toString().split(",");
for(int y = 0; y < listModel.size(); y++)
{
String[] b = listModel.get(y).toString().split(",");
if(a[2].compareTo(b[2]) > 0 && a[1].equals(b[1]) != true)
{
tmp.add(0,listModel.get(y));
listModel.set(y,listModel.get(x));
listModel.set(x,tmp.get(0));
}
}
}
}
不超过规定值。我已经结束了六次迭代。至于值为 0 的被数,我决定不做任何计算就返回 0。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。