Fixed point multiply/divide for 15.16 numbers

I am looking for an algorithm to multiply and divide fixed point 15.16 numbers.

I already have addition and subtraction. Those were easy - simple 32-bit add and subtract. With multiply and divide, I can also add many trigonometric and exponential/log functions. And I think I can deal with just multiply, as my library has a reciprocal function and I can use that to implement division: a * (1/b) = a / b. But a 32-bit multiply does not work as it ignores the radix point.

I am working on a 16-bit microcontroller, so I would like to avoid anything more than 32-bit multiply, which takes about 4 cycles on my processor. It's not crucial thou开发者_高级运维gh, I'm just trying to replace floating point math.

I have heard I need to shift or rotate the result, but I am not sure how this would help or specifically how to shift it. Any suggestions or help appreciated!

Think of is this way: your number a.b is represented as (a.b * 65536)

If you multiply a.b * c.d the value you get is (a.b * 65536) * (c.d * 65536), so to put this back in the right representation you need to divide by 65536.

When you divide a.b / c.d the value you get is (a.b * 65536) / (c.d * 65536), so to put this back in the right representation you need to multiply by 65536. You should multiply by 65536 before the divide to preserve as many bits as possible in the result.

Of course you can substitute (<< 16) for (* 65536) if that is faster on your processor. Similarly you can substitute (>> 16) for (/ 65536).

Here's a.b * c.d:

uint32_t result_low = (b * d);
uint32_t result_mid = (a * d) + (b * c);
uint32_t result_high = (a * c); 
uint32_t result = (result_high << 16) + result_mid + (result_low >> 16)

First the theory: assuming signed numbers, multiplying a Q15.16 by another Q15.16 will give you a Q(15+15+1).(16+16) = Q31.32 number. Thus, you need a 64-bit integer variable to hold the result.

If your compiler has a 64-bit integer type, just use it and let the compiler figure out how to do the 32-bit x 32-bit multiplication on the 16-bit CPU (that's what compilers are for):

int32_t a_15q16, b_15q16;
int64_t res_31q32 = (int64_t)a_15q16 * (int64_t)b_15q16;

What you do afterwards with the Q31.32 result really depends on your application.

You may be wondering why the result needs 31 integer instead of 30. In fact, the additional bit is needed only for the case where you multiply -2^15 by -2^15. If your operands are guaranteed never to equal -2^15 simultaneously, you can assume a Q30.32 result.

To find out if your compiler supports 64-bit integers, you may need to look at the compiler manual. If this is a C99 compiler, see if there is a int64_t type in your stdint.h header.

The multiplication is easily done with 64-bit multiply: (a * b) >> 16. Division, similarly, is easily done with 64 bits: (a << 16) / b. Depending on your rounding / error requirements, you might want to complicate that slightly in order to get the last bit of the output correct.