c - 32-bit signed integer multiplication without using 64-bit data type

Question

I want to do 32-bit signed integer multiplication without using a 64-bit data type. My inputs are in Q1.31 (both) format.

input1 = A32 (Ah Al) - higher, lower half's of A32
input2 = B32 (Bh Bl) - higher, lower half's of B32

Result should be in Q1.31 format, leave the overflow case.

I need C code. Please provide the explanation with formats also.

Answer 1

Signed Q1.31 format is a fully fractional format capable of representing operands between -1 and almost +1. The scale factor is 2³¹. This means that when each Q1.31 operand is stored in a 32-bit signed integer, we can generate the Q1.31 product by computing the full double-width product of the signed integers, then right shifting the result by 31 bits. The right shift is necessary because the product includes the scale factor twice, and the shift acts as a division that removes one instance of the scale factor.

We can compute the double-width product of two 32-bit integers by separately computing the upper and lower 32 bits of the full product. The lower 32 bits are computed trivially as the ordinary product of the two inputs. To compute the upper 32 bits, we need to write a function mul32hi(). In order to avoid using a wider type (i.e. one using more than 32 bits) in intermediate computations, we need to split the original operands into halves, compute their partial products, and then sum these partial products appropriately.

Note that various processors provide a hardware instruction that implements the functionality of mul32hi(). In this case one would want to use an appropriate intrinsic, or a bit of inline assembly code if no intrinsic exists, rather than use the emulation code presented here.

It helps to reduce the problem first to the corresponding unsigned multiplication, umul32hi(), then derive the signed result from that via the definition of 2's complement representation (which is assumed in the following C code):

#include <stdint.h>

/* compute the upper 32 bits of the product of two unsigned 32-bit integers */
uint32_t umul32hi (uint32_t a, uint32_t b)
{
    /* split operands into halves */
    uint32_t al = (uint16_t)a;
    uint32_t ah = a >> 16;
    uint32_t bl = (uint16_t)b;
    uint32_t bh = b >> 16;
    /* compute partial products */
    uint32_t p0 = al * bl;
    uint32_t p1 = al * bh;
    uint32_t p2 = ah * bl;
    uint32_t p3 = ah * bh;
    /* sum partial products */
    uint32_t cy = ((p0 >> 16) + (uint16_t)p1 + (uint16_t)p2) >> 16;
    return p3 + (p2 >> 16) + (p1 >> 16) + cy;
}

/* compute the upper 32 bits of the product of two signed 32-bit integers */
int32_t mul32hi (int32_t a, int32_t b)
{
    return umul32hi (a, b) - ((a < 0) ? b : 0) - ((b < 0) ? a : 0);
}

/* compute the full 64-bit product of two signed 32-bit integers */
void mul32wide (int32_t a, int32_t b, int32_t *rhi, int32_t *rlo)
{
    *rlo = a * b;           /* bits <31:0> of the product a * b */
    *rhi = mul32hi (a, b);  /* bits <63:32> of the product a * b */
}

/* compute the product of two signed Q1.31 fixed-point numbers */    
int32_t mul_q_1_31 (int32_t a, int32_t b)
{
    int32_t hi, lo;
    mul32wide (a, b, &hi, &lo);
    /* Q1.31 is scaled by 2**31, trim out scale factor */
    return (int32_t)(((uint32_t)hi << 1) | ((uint32_t)lo >> 31));
}

I interpreted the request to "leave the overflow case" to mean to ignore overflow. As a consequence, multiplying -1 (0x80000000) by -1 (0x80000000) with mul_q_1_31() will return -1 (0x80000000).