得到整数的高半部分和低半乘

Getting the high half and low half of a full integer multiply

本文关键字：高半部 整数      更新时间：2023-10-16

我以三个值 A，B，C(无符号 32 位整数(开头。我必须获得两个值 D，E(也是无符号的 32 位整数(。哪里

D = high(A*C);
E = low(A*C) + high(B*C);

我希望两个 32 位 uint 的乘以会产生 64 位结果。"高"和"低"只是我用 64 位结果标记前 32 位和后 32 位的 64 位。

我尝试获得一些 allready 功能代码的优化代码。我在巨大的循环中有一小段代码，它只有几个命令行，但它消耗了几乎所有的计算时间(几个小时计算的物理模拟(。这就是为什么我尝试优化这一小部分的原因，其余代码可以保持更"用户良好安排"。

有一些 SSE 指令适合上述计算例程。gcc编译器可能会做优化的工作。但是，如果有必要，我不拒绝直接在 SSE 入侵中编写一些代码的选项。

请耐心等待我对SSE的低经验。我将尝试以符号方式为 SSE 编写算法。订购口罩或理解结构可能会有一些错误。

  1. 将四个 32 位整数按顺序存储到一个 128 位寄存器中：A、B、C、C。
  2. 将指令(可能是 pmuludq(应用到上述 128 位寄存器中，该寄存器将 32 位整数对相乘并返回 64 位整数对作为结果。因此，它应该同时计算A*C乘法和B*C乘法并返回两个 64 位值。
  3. 我希望我有新的 128 位寄存器值 P，Q，R，S(四个 32 位块(，其中 P，Q 是 A*C 的 64 位结果，R，S 是 B*C 的 64 位结果。然后我继续将寄存器处的值重新排列为顺序 P，Q，0，R
    4. 取第一个 64 位
  4. P，Q 并添加第二个 64 位 0，R。结果是一个新的 64 位值。
  5. 将结果的前 32 位读取为 D，将结果的最后 32 位读取为 E。

此算法应返回正确的 E 和 D 值。

我的问题：

c++ 中是否有静态代码可以生成与上述 1-5 SSE 算法类似的 SSE 例程？我选择具有更高性能的解决方案。如果算法对于标准 c++ 命令有问题，有没有办法在 SSE 中编写算法？

我使用 TDM-GCC 4.9.2 64 位编译器。

(注：问题经建议后修改(

(注2：我从这个 http://sci.tuomastonteri.fi/programming/sse 中得到了使用SSE以获得更好性能的灵感(

除非有多个输入要并行处理，否则不需要向量。 Clang和GCC在优化编写代码的"正常"方式方面已经做得很好：投射到两倍的大小，乘以，然后移位以获得高半。 编译器可识别此模式。

他们注意到操作数开始时为 32bit，因此在转换为 64b 后，上半部分全部为零。 因此，他们可以使用 x86 的 mul insn 进行 32b*32b->64b 乘法，而不是执行完全扩展精度的 64b 乘法。 在 64 位模式下，它们对__uint128_t版本的代码执行相同的操作。

这两个函数都可以编译成相当好的代码(每个乘法一个mul或imul(。 gcc -m32不支持 128b 类型，但我不会讨论它，因为 1.您只询问了 32 位值的全乘法和 2。当您希望某些内容快速运行时，您应该始终使用 64 位代码。 如果你正在做全乘法，结果不适合寄存器，clang会避免很多额外的mov指令，因为gcc对此很愚蠢。 这个小测试函数为 gcc 错误报告提供了一个很好的测试用例。

该 godbolt 链接包含一个在循环中调用 this 的函数，将结果存储在数组中。 它通过一堆随机播放自动矢量化，但如果有多个输入要并行处理，它看起来仍然像加速。 不同的输出格式在乘法后可能会减少洗牌，例如可能存储单独的数组用于D和E。

我包括 128b 版本是为了表明编译器可以处理这个问题，即使它不是微不足道的(例如，只需执行 64 位imul指令即可在 32 位输入上执行 64*64->64b 乘法，在将函数输入时可能位于输入寄存器中的任何高位归零后。

当针对 Haswell CPU 及更新版本时，gcc 和 clang 可以使用 mulx BMI2 指令。 (我在 godbolt 链接中使用了-mno-bmi2 -mno-avx2来保持 asm 更简单。 如果你有一个哈斯韦尔CPU，只需使用-O3 -march=haswell。 mulx dest1, dest2, src1 dest1:dest2 = rdx * src1，而mul src1 rdx:rax = rax * src1。 所以mulx有两个只读输入(一个隐式：edx/rdx (，和两个只写输出。 这使编译器可以使用更少的mov指令执行全乘法，以将数据传入和传出mul的隐式寄存器。 这只是一个小的加速，特别是因为 64 位mulx在 Haswell 上有 4 个周期延迟而不是 3 个周期延迟。 (奇怪的是，64位mul和mulx比32位mul和mulx略便宜。

// compiles to good code: you can and should do this sort of thing:
#include <stdint.h>
struct DE { uint32_t D,E; };
struct DE f_structret(uint32_t A, uint32_t B, uint32_t C) {
  uint64_t AC = A * (uint64_t)C;
  uint64_t BC = B * (uint64_t)C;
  uint32_t D = AC >> 32;         // high half
  uint32_t E = AC + (BC >> 32);  // We could cast to uint32_t before adding, but don't need to
  struct DE retval = { D, E };
  return retval;
}

#ifdef __SIZEOF_INT128__  // IDK the "correct" way to detect __int128_t support
struct DE64 { uint64_t D,E; };
struct DE64 f64_structret(uint64_t A, uint64_t B, uint64_t C) {
  __uint128_t AC = A * (__uint128_t)C;
  __uint128_t BC = B * (__uint128_t)C;
  uint64_t D = AC >> 64;         // high half
  uint64_t E = AC + (BC >> 64);
  struct DE64 retval = { D, E };
  return retval;
}
#endif

如果我理解正确，您希望计算 A*B 中潜在溢出的数量。如果是，那么您有 2 个不错的选择 - "使用两倍大的变量"(为 uint64 编写 128 位数学函数 - 这并不难(或等待我明天发布((和"使用浮点类型"：(浮点数(A(*浮点数(B((/浮点数(C(由于精度损失最小(假设浮点数为 4 字节，双精度 8 字节，长双精度 16 字节长(，并且浮点数和 uint32 都需要 4 字节的内存(使用双精度表示uint64_t因为它应该是 8 字节长(：

#include <iostream>
#include <conio.h>
#include <stdint.h>
using namespace std;
int main()
{
    uint32_t a(-1), b(-1);
    uint64_t result1;
    float result2;
    result1 = uint64_t(a)*uint64_t(b)/4294967296ull;    // >>32 would be faster and less memory consuming
    result2 = float(a)*float(b)/4294967296.0f;
    cout.precision(20);
    cout<<result1<<'n'<<result2;
    getch();
    return 0;
}

生产：

4294967294
4294967296

但是如果你想要真正精确和正确的答案，我建议使用两倍的大类型进行计算。

现在我想到了 - 你可以对 uint64 使用长双精度，对

uint32 使用双精度，而不是为 uint64 编写函数，但我认为不能保证长双精度是 128 位，你必须检查它。我会选择更通用的选择。

编辑：

You can write function to calculate that without using anything more
than A, B and result variable which would be of the same type as A.
Just add rightmost bit of (where Z equals B*(A>>pass_number&1)) Z<<0,
Z<<1, Z<<2 (...) Z<<X in first pass, Z<<-1, Z<<0, Z<<1 (...) Z<<(X-1)
for second (there should be X passes), while right shifting the result
by 1 (the just computed bit becomes irrelevant to us after it's
computed as it won't participate in calculation anymore, and it would
be erased anyway after dividing by 2^X (doing >>X)

(不得不放入"代码"，因为我是新来的，找不到另一种方法来防止格式化脚本吃掉一半(

这只是一个快速的想法。您必须检查它的正确性(对不起，但我现在真的很累 - 但结果不应该在任何计算点溢出，因为如果我是正确的，最大进位的值将是 2X，并且算法本身似乎很好(。

如果您仍然需要帮助，我明天会为此编写代码。

相关文章：