vrecpeq_f32内在的参考实现?

Question

有vrecpeq_f32ARM NEON Intrinsic

vrecpeq_f32的官方解释:https://developer.arm.com/architectures/instruction-sets/intrinsics/#f:@navigationhierarchiessimdisa=[Neon]&q=vrecpeq_f32 .

浮点倒数估计。该指令查找源simdfp寄存器中每个向量元素的近似倒数估计，将结果放入一个向量中，并将该向量写入目标simdfp寄存器。

然而，这对我来说仍然不准确。只是想知道我们是否可以在C/c++中编写一个参考实现，保持与vrecpeq_f32完全相同的结果?

我试着打电话给vrecpeq_f32并得到结果:

float32x4_t v1 = {1, 2, 3, 4};
float32x4_t v_out = vrecpeq_f32(v1);//0.99805, 0.49902, 0.33301, 0.24951

好奇为什么1的倒数是0.99805而不是1.0。

注:我对如何使用NEON固有特性和一些技巧来获得更好的精度互反结果不感兴趣，例如一个或多个牛顿-拉夫森迭代。

Answer 1

ARM文档提供了详细描述正在执行的精确算法的伪代码。查找使用定点RecipEstimate的FPRecipEstimate。

这可能看起来像很多代码，但其中很大一部分是用来处理各种边缘情况、操作模式和元素大小的。

只是想知道我们是否可以在C/c++中编写一个参考实现，保持与vrecpeq_f32完全相同的结果?

当然!毕竟，这可以归结为位操作，所以没有理由认为它不可行。将其转换为c++，同时删除大多数边缘情况处理以及扩展精度模式，如下所示:(参见godbolt)

免责声明这不是该函数的完整实现，仅足以探索精度行为，假设有限的规范化输入，没有特殊情况。不要在代码库中删除它，期望它通常与指令匹配。

#include <iostream>
#include <cstring>
#include <iomanip>
// Convenience struct to deal with encoding and decoding ieee754 floats
struct float_parts {
explicit float_parts(float v);
explicit operator float() const;
std::uint32_t sign;
std::uint32_t fraction;
std::uint32_t exp;
};
// Adapted from:
// https://developer.arm.com/documentation/ddi0596/2021-03/Shared-Pseudocode/Shared-Functions?lang=en#impl-shared.FPRecipEstimate.2
// RecipEstimate()
// ===============
// Compute estimate of reciprocal of 9-bit fixed-point number.
//
// a is in range 256 .. 511 representing a number in
// the range 0.5 <= x < 1.0.
// result is in the range 256 .. 511 representing a
// number in the range 1.0 to 511/256
std::uint32_t RecipEstimate(std::uint32_t a) {
a = a*2+1;
std::uint32_t b = (1 << 19) / a;
return ( b + 1) / 2;
}
// FPRecipEstimate()
// =================
float FPRecipEstimate(float operand) {
// ([...],sign,[...]) = FPUnpack(operand, [...], [...]);
// fraction = operand<22:0> : Zeros(29);
// exp = UInt(operand<30:23>);
float_parts parts{operand};    
// scaled = UInt('1':fraction<51:44>);
std::uint32_t scaled = 0x100 | ((parts.fraction >> 15) & 0xFF) ;
// when 32 result_exp =  253 - exp; // In range 253-254 = -1 to 253+1 = 254
parts.exp = 253 - parts.exp;
// // Scaled is in range 256 .. 511 representing a
// // fixed-point number in range [0.5 .. 1.0].
// estimate = RecipEstimate(scaled, increasedprecision);
std::uint32_t estimate = RecipEstimate(scaled);
// fraction = estimate<11:0> : Zeros(40);
parts.fraction = (estimate & 0xff ) << 15;
return float(parts);
}
int main() {
std::cout << std::setprecision(5) 
<< FPRecipEstimate(1.0f) << "n"
<< FPRecipEstimate(2.0f) << "n"
<< FPRecipEstimate(3.0f) << "n"
<< FPRecipEstimate(4.0f);
}
float_parts::float_parts(float v) {
std::uint32_t v_bits;
std::memcpy(&v_bits, &v, sizeof(float));
sign = (v_bits >> 31) & 0x1;
fraction = v_bits & ((1 << 23) - 1);
exp = (v_bits >> 23) & 0xff;
}
float_parts::operator float() const {
std::uint32_t v_bits = 
((sign & 0x1) << 31) |
(fraction & ((1 << 23) - 1)) |
((exp & 0xff) << 23);
float result;
std::memcpy(&result, &v_bits, sizeof(float));
return result;
}

生成期望值:

0.99805
0.49902
0.33301
0.24951