平方根计算算法

Question

我一直在用 C 语言实现控制软件，其中一种控制算法需要平方根计算。我一直在寻找合适的平方根计算算法，它具有恒定的执行时间，而与基值无关。此要求排除了标准库中的 sqrt 函数。

就我的平台而言，我一直在使用基于浮点 32 位 ARM Cortex A9 的机器。就我的应用程序中的radicand范围而言，算法是以物理单位计算的，所以我希望遵循范围<0, 400>。至于所需的误差，我认为大约 1% 的误差就足够了。谁能推荐一个适合我目的的平方根计算算法？

Answer 1

我决定使用以下方法。我选择了牛顿法，然后我通过实验设置了固定的迭代次数，以使整个radicand范围内的误差，即<0,400>不超过规定值。我已经结束了六次迭代。至于值为 0 的基数，我决定不进行任何计算就返回 0。

Answer 2

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, vdup_n_f32(0.0f)));
    float32x2_t a2est = vrsqrte_f32(a2);
    a2est = vand_u32(is_legal, a2est);
    return vmul_f32(a2, a2est);
}

这个实现给出了sqrtest2(1) == 0.998 和sqrtest2(400) == 19.97（在带有arm64 的MacBook M1 上测试）。由于无分支且无 LUT，这可能具有恒定的执行时间，假设所有指令都以恒定数量的周期执行。

Answer 3

我最初的方法是使用泰勒级数求平方根，并在多个固定点上预先计算系数。这会将计算减少为减法和乘法。

查找表将是一个二维数组，例如：

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 in [0 - 1] -> 距离 0.1，x in [1 - 2] -> 距离 0.25，x in [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