Libfixmath 中的乘法是如何执行的？

Question

我正在 CMD Batch 中移植 Wikipedia https://en.wikipedia.org/wiki/Libfixmath 提到的库。

该库使用 Q16.16 格式。

我正在尝试将其移植为 Q16.16 或 Q15.15 格式。

我看不懂的是fix16.c中的代码部分。

支持32bit的硬件情况下的乘法。

fix16_t fix16_mul(fix16_t inArg0, fix16_t inArg1)
{
    // Each argument is divided to 16-bit parts.
    //                  AB
    //          *    CD
    // -----------
    //                  BD  16 * 16 -> 32 bit products
    //               CB
    //               AD
    //              AC
    //           |----| 64 bit product
    int32_t A = (inArg0 >> 16), C = (inArg1 >> 16);
    uint32_t B = (inArg0 & 0xFFFF), D = (inArg1 & 0xFFFF);
    
    int32_t AC = A*C;
    int32_t AD_CB = A*D + C*B;
    uint32_t BD = B*D;
    
    int32_t product_hi = AC + (AD_CB >> 16);
    
    // Handle carry from lower 32 bits to upper part of result.
    uint32_t ad_cb_temp = AD_CB << 16;
    uint32_t product_lo = BD + ad_cb_temp;
    if (product_lo < BD)
        product_hi++;
...

A、B、C、D 是 16 位的数字。 我的问题是在计算进位时我不明白它是如何做的。 它取ad_cb并进行移位，即取AD_CB的低16位，将其上移并添加到BD中。 第三个操作我看不懂。

if (product_lo < BD)
        product_hi++;
...

你能更好地解释一下吗？和前面的有什么关系？

部分移植：

:fix15_mul F1 F2 Ret

setlocal

    set "F1=%1"
    set "F2=%2"

    set /A "A=F1>>15, C=F2>>15"

    set /A "B=F1 & 0x7FFF, D=F2 & 0x7FFF"
    
    set /A "AC=A*C, AD_CB=A*D + C*B,BD =B*D"
    
    set /A "product_hi = AC + (AD_CB >> 15)"
    
    rem Handle carry from lower 30 bits to upper part of result.
    set /A "ad_cb_temp=AD_CB << 15, product_lo = BD + ad_cb_temp"

    set /A "product_lo2 =(AD_CB & 0x7FFF)<<15"

    if %product_lo% lss %BD% set /A "product_hi+=1"

    set /A "return=(product_hi << 15) | (product_lo >> 15)"

( endlocal & set "%3=%Return%" )

goto :eof

编辑：仅一个焦点

Answer 1

我将尝试解释并在下面提供一个带注释的修改代码示例（未经测试），我认为它可能有效。

fix16_t fix16_mul(fix16_t inArg0, fix16_t inArg1)
{
    // Each argument is divided to 16-bit parts.
    //                  AB   *    CD
    // -----------
    //                BD  16 * 16 -> 32 bit products
    //               CB
    //               AD
    //              AC
    //             |----| 64 bit product
    int32_t  A = (inArg0 >> 16),    C = (inArg1 >> 16);
    uint32_t B = (inArg0 & 0xFFFF), D = (inArg1 & 0xFFFF);
    
    int32_t AC = A*C; // worst case bounds -2^15*(2^15-1), 2^30
    int32_t AD = A*D; // worst case bounds -2^15*(2^16-1),(2^15-1)*(2^16-1)
    int32_t BC = B*C; // ditto. The sum of these two terms could overflow 
    uint32_t BD = B*D; // bounds 0 .. (2^16-1)^2
    
    int32_t product_hi = AC + (AD >> 16) + (BC >> 16);
    
    // Handle carry from lower 32 bits to upper part of result.
    uint32_t ad_bc_temp = (AD << 16) + (BC << 16); // as unsigned
    uint32_t product_lo = BD + ad_bc_temp;
    // these are both now unsigned ints so that when adding together
    // if an overflow occurs then product_lo < BD
    if (product_lo < BD)
        product_hi++;   // so we increment the product_hi
    ...

“进位我不明白它是怎么做的。它取ad_cb并进行移位，即取AD_CB的低16位并将它们移到楼上并添加到BD中。第三个操作我无法理解。 ”

它依赖于 BD 的 unsigned int 值在发生进位时总是变小（并且这里只可能发生一次进位）。我建议您针对一些尴尬的边缘情况验证此代码，并与完整 int64 算术中完成的计算进行比较。值 0x8000ffff 和 0x7FFFFFFF 是明显的压力测试。在做这些事情时，很容易出现粗心的失误，只影响少数边缘情况。