时间戳

首页 > 解决方案 > 平方根计算算法

algorithm - 平方根计算算法

问题描述

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

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

标签： algorithmmathembeddedsqrt

解决方案

vrsqrte_f32Arm v7 指令集为两个同时逼近和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==0vrsqrtes_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，这可能具有恒定的执行时间，假设所有指令都以恒定数量的周期执行。

推荐阅读