问题描述
我一直在用 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。