计算两个字的有符号双字乘积，已知低位字乘积。

Question

计算两个字的有符号双字乘积，已知低位字乘积。

5

在《黑客的乐趣》一书中，有一个算法可以计算两个（有符号）字的双字乘积。

函数muldws1使用四次乘法和五次加法来计算两个字的双字。

在该代码的末尾有一行被注释掉的代码。

/* w[1] = u*v;                  // Alternative. */

这种替代方法使用了五次乘法和四次加法，即将一次加法换成了一次乘法。

但我认为这种替代方法可以改进。我还没有提到硬件方面的内容。假设有一个假想的CPU，它可以计算两个字的乘积的低位，但不能计算高位（例如对于32位字，32x32的低32位）。在这种情况下，我觉得这个算法可以改进。以下是我的构想，假设是32位字（相同的概念也适用于64位字）。

void muldws1_improved(int w[], int32_t x, int32_t y) {
    uint16_t xl = x; int16_t xh = x >> 16;
    uint16_t yl = y; int16_t yh = y >> 16;

    uint32 lo = x*y;
    int32_t t = xl*yh + xh*yl;

    uint16_t tl = t; int16_t th = t >>16;
    uint16_t loh = lo >> 16;

    int32_t cy = loh<tl; //carry
    int32_t hi = xh*yh + th + cy;
    w[0] = hi; w[1] = lo;
}

这个方法使用了四次乘法，三次加法和一次比较。这个改进比我希望的要小。

这能被改进吗？有更好的方法来确定进位标志吗？ 我应该指出，我还假设硬件没有进位标志（例如没有ADDC指令），但可以比较字（例如word1<word）。

编辑：正如Sander De Dycker所指出的那样，我的函数未通过单元测试。这里有一个通过单元测试但效率较低的版本。我认为它可以改进。

void muldws1_improved_v2(int w[], int32_t x, int32_t y) {
    uint16_t xl = x; int16_t xh = x >> 16;
    uint16_t yl = y; int16_t yh = y >> 16;

    uint32_t lo = x*y;
    int32_t  t2 = xl*yh;
    int32_t  t3 = xh*yl;
    int32_t  t4 = xh*yh;

    uint16_t t2l = t2; int16_t t2h = t2 >>16;
    uint16_t t3l = t3; int16_t t3h = t3 >>16;
    uint16_t loh = lo >> 16;

    uint16_t t = t2l + t3l;
    int32_t carry = (t<t2l) + (loh<t);
    int32_t hi = t4 + t2h + t3h + carry;
    w[0] = hi; w[1] = lo;
}

这个函数使用了四次乘法、五次加法和两次比较，比原来的函数更劣。

- Z boson

如果您的代码是正确的，那么它也会胜过使用更少移位和算术&操作的黑客之乐版本。当然，这可能取决于 CPU 架构，16 位移位可能是免费或非常便宜的。 - John Bollinger

你的代码未通过原始代码中的 {0x7fffffff, 0x7eeeeeee, 0x3f777776,0x81111112} 单元测试。 - Sander De Dycker

顺便提一下：该函数被定义为返回 int64_t，但实际上它并不返回任何值。 - joop

@joop，是的，我知道，让我来修复它。在我的代码中，我将其转换为int64_t并返回以进行比较，但我修改了它，使其看起来像来自Hacker's Delight的函数。 - Z boson

1

小问题：int w[] 应该改为 int32_t w[] 或者更好的方式是：int32_t *whi, uint32_t *wlo。 - chux - Reinstate Monica

显示剩余3条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Z boson · Accepted Answer

我的问题中的muldws1_improved函数有两个问题。其中一个是当我执行xl*yh + xh*yl时，它错过了进位。这就是为什么它未通过单元测试的原因。~~但另一个问题是有一些带符号*无符号乘积需要比C代码中看到的机器逻辑更多。~~（见下面的编辑）。我找到了一个更好的解决方案, 就是先优化无符号乘积函数muldwu1，然后再进行操作。

muldwu1(w,x,y);
w[0] -= ((x<0) ? y : 0)  + ((y<0) ? x : 0);

我尝试使用低位 lo = x*y 来改进 muldwu1（是的，这个函数通过了 Hacker's delight 的单元测试），以纠正符号。

void muldwu1_improved(uint32_t w[], uint32_t x, uint32_t y) {
    uint16_t xl = x; uint16_t xh = x >> 16;
    uint16_t yl = y; uint16_t yh = y >> 16;

    uint32_t lo   = x*y;    //32x32 to 32
    uint32_t t1   = xl*yh;  //16x16 to 32
    uint32_t t2   = xh*yl;  //16x16 to 32
    uint32_t t3   = xh*yh;  //16x16 to 32

    uint32_t t    = t1 + t2;
    uint32_t tl   = 0xFFFF & t;
    uint32_t th   = t >> 16;
    uint32_t loh  = lo >> 16;

    uint32_t cy   = ((t<t1) << 16) + (loh<tl); //carry
             w[1] = lo;
             w[0] = t3 + th + cy;
}

这个函数比《黑客秘笈》中的原始函数少了一个加法，但需要进行两次比较。

 1 mul32x32 to 32
 3 mul16x16 to 32
 4 add32
 5 shift logical (or shuffles)
 1 and
 2 compare32
***********
16 operations

编辑：

我对《黑客的趣味》（第二版）中关于mulhs和mulhu算法的一句话感到困扰。

该算法需要16个基本的RISC指令，无论是有符号还是无符号版本，其中有四个乘法指令。

我在仅使用16个SSE指令实现了无符号算法，但我的有符号版本需要更多的指令。我找到了原因，现在可以回答自己的问题了。

我没有找到比《黑客的趣味》中更好的版本的原因是他们假设的RISC处理器有一个计算两个字的乘积的低位字的指令。换句话说，他们的算法已经针对这种情况进行了优化，因此很难找到比他们已经有的更好的版本。

他们列出备选方案的原因是他们认为乘法（和除法）可能比其他指令更昂贵，因此他们将替代方案留作优化的情况。

因此，C代码并没有隐藏重要的机器逻辑。它假设机器可以将字*字计算为低位字。

为什么这很重要？在他们的算法中，他们首先执行

u0 = u >> 16;

以后

t = u0*v1 + k;

如果u = 0x80000000，则u0 = 0xffff8000。但是，如果您的CPU只能使用半字乘积得到完整字，则u0的高半字将被忽略，您将得到错误的有符号结果。

在这种情况下，您应该计算无符号的上半字，然后使用hi -= ((x<0) ? y : 0) + ((y<0) ? x : 0);进行校正，就像我之前所述的那样。

我感兴趣的原因是因为英特尔的SIMD指令（从SSE2到AVX2）没有一个可以执行64x64到64的指令，它们只有32x32到64。这就是为什么我的有符号版本需要更多的指令的原因。

但是，AVX512有一个64x64到64的指令。因此，在AVX512中，有符号版本应该需要与无符号版本相同数量的指令。然而，由于64x64到64的指令可能比32x32到64的指令慢得多，因此执行无符号版本并进行校正可能更有意义。