使用 pre-SSE4 计算 vector2 double 的 floor 和 ceil答案

【问题标题】：Calculating floor & ceil of vector2 double using pre-SSE4使用 pre-SSE4 计算 vector2 double 的 floor 和 ceil
【发布时间】：2020-07-28 16:28:22
【问题描述】：

这可以通过 sse4.1 内在函数 _mm_floor_pd 和 _mm_ceil_pd 来完成翻译成roundpd xmm,xmm,1和roundpd xmm,xmm,2

使用 SSE/SSE2/SSE3 计算的最佳方法是什么？

【问题讨论】：

如果您的值不会溢出 int32_t，您可以使用打包转换到 int 并返回。对于非负值，trunc = floor，因此您可以保存/恢复 cvttpd2dq / cvtdq2pd 周围的符号位。或者添加和减去一个巨大的正数（或负数）可能会更快，但我认为这会让你离你最近，而不是地板/天花板，没有其他技巧。
这对你有帮助吗Fast Rounding of Floating Point Numbers...？还可以查看现代编译器的汇编输出floor/ceil at godbolt。
相关：Efficient SSE FP `floor()` / `ceil()` / `round()` Rounding Functions Without SSE4.1? - 这是单精度的，它转换为整数并返回（因此对于全范围的浮点数甚至不安全）。您可以利用任何有限范围的因素，例如幅度 How to efficiently perform double/int64 conversions with SSE/AVX? 已使用 bithacks 打包转换为 int64 / uint64。
我想问的正是双打。标题已编辑
floor/ceil for double vector2 at godbolt

标签： c++ assembly sse simd intrinsics

【解决方案1】：

这是在 SSE4.1 之前的 CPU 上执行 floor/ceil 的代码。请注意，使用 '-ffast-math' 会破坏它！

#include <cmath>
#include <emmintrin.h>
#include <cstdio> // for printf

#ifdef _MSC_VER
#define __attribute__(P)
#endif

struct vec2d {
    double x;
    double y;
};

static __m128d mm_blendv_pd(const __m128d& a, const __m128d& b, const __m128d& mask) noexcept {
  return _mm_or_pd(_mm_andnot_pd(mask, a), _mm_and_pd(b, mask));
}

__attribute__((optimize("-fno-associative-math")))
vec2d _floor(vec2d v) noexcept {
  __m128d src = _mm_set_pd(v.x, v.y);
  __m128d maxn = _mm_set_pd(4503599627370496.0, 4503599627370496.0);  // pow(2, 52)
  __m128d magic = _mm_set_pd(6755399441055744.0, 6755399441055744.0); // pow(2, 52) + pow(2, 51)
  __m128d msk = _mm_cmpnlt_pd(src, maxn);
  __m128d rounded = _mm_sub_pd(_mm_add_pd(src, magic), magic); //! -ffast-math will break this!
  __m128d maybeone = _mm_and_pd(_mm_cmplt_pd(src, rounded), _mm_set_pd(1.0, 1.0));
  __m128d res = mm_blendv_pd(_mm_sub_pd(rounded, maybeone), src, msk);
  return {_mm_cvtsd_f64(_mm_unpackhi_pd(res, res)), _mm_cvtsd_f64(res)};
}

__attribute__((optimize("-fno-associative-math")))
vec2d _ceil(vec2d v) noexcept {
  __m128d src = _mm_set_pd(v.x, v.y);
  __m128d maxn = _mm_set_pd(4503599627370496.0, 4503599627370496.0);  // pow(2, 52)
  __m128d magic = _mm_set_pd(6755399441055744.0, 6755399441055744.0); // pow(2, 52) + pow(2, 51)
  __m128d msk = _mm_cmpnlt_pd(src, maxn);
  __m128d rounded = _mm_sub_pd(_mm_add_pd(src, magic), magic); //! -ffast-math will break this!
  __m128d maybeone = _mm_and_pd(_mm_cmpnle_pd(src, rounded), _mm_set_pd(1.0, 1.0));
  __m128d res = mm_blendv_pd(_mm_add_pd(rounded, maybeone), src, msk);
  return {_mm_cvtsd_f64(_mm_unpackhi_pd(res, res)), _mm_cvtsd_f64(res)};
}


int main() {
    vec2d v{3.1,4.6};

    vec2d v2 = _floor(v);
    vec2d v3 = _ceil(v);

    printf(" v2: %f %f\n",v2.x,v2.y);
    printf(" v3: %f %f\n",v3.x,v3.y);

    return 0;
}

Live Code

它基于这篇博文C++ Compilers and FP Rounding on X86 ，但是那里的代码有错误。

【讨论】：

-4503599627370495.5 失败
@chtz：感谢您报告此事。这意味着这段代码不能普遍使用:(