#include <stdint.h>
#include <stdio.h>
typedef int32_t fix_t;
#define fix_ONE(w) ((fix_t)(0x00000001 << w))
#define MASK_M(w) (((fix_t)1 << w) - 1) // mask of mantissa, all LSB set, all MSB clear
#define fix_rconst(x, w) ((fix_t)((x) * fix_ONE(w) + ((x) >= 0 ? 0.5 : -0.5))) // convert a const to rounded fix_t
static const fix_t fix_overflow = 0x80000000; // 1000 0000 0000 0000 0000 0000 0000 0000
static inline float fix_to_float(fix_t a, int w) { return ((float)a / fix_ONE(w)); }
fix_t fix_mul(fix_t inArg0, fix_t inArg1, int w){
// w is fractional bitwidth
// separate interger part and fraciton part
int32_t A = (inArg0 >> w), C = (inArg1 >> w);
uint32_t B = (inArg0 & MASK_M(w)), D = (inArg1 & MASK_M(w));
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 >> w); // product_hi stands for the interger part of final result
uint32_t ad_cb_temp = AD_CB << (32-w); // get the fraction part of AD_CB
uint32_t product_lo = BD + ad_cb_temp; //product_lo stands for the fraction part of final result
if (product_lo < BD || product_lo < ad_cb_temp){ // check if product_lo overflow
product_hi++;
}
// The upper part bits should all be the same (including the sign).
if (product_hi >> 31 != product_hi >> (31-w)){
printf("Overflow in fix_mul(), please use other bitwidth \n");
return fix_overflow;
}
// combine interger part and fraction part
return (product_hi << (w)) | (product_lo >> (32-w));
}
int test_mul(void)
{
// test cases
float a = 0.267f; //0.50f;//1.267f; //-1.267f;//-1.267f;
float b = 0.849f; //0.25f;//1.849f; //1.849f; //-1.849f;
for (int w = 1; w < 28; w++) {
fix_t aa = fix_rconst(a, w);
fix_t bb = fix_rconst(b, w);
fix_t result = fix_mul(aa, bb, w);
float fresult = fix_to_float(result, w);
printf("fix_rconst(%f, %i) = %i, fix_rconst(%f, %i) = %i, result = %i, fresult=%f \n", a, w, aa, b, w, bb, result, fresult);
}
return 0;
}
int main()
{
test_mul();
// system("pause");
return 0;
}
