working fp32 multiplier using 8x8 multipliers

[suarezvictor]

This code calculates FP32 mult using 6 instances of 8x8 multipliers

//Calculates FP32 mult using 6 instances of 8x8 multipliers
//(C) 2021 Victor Suarez Rovere <[email protected]>
//LICENSE: GPLv3
//
//Compile & run:
//g++ -O3 -o fptest fp32mul.c && ./fptest

#include <stdio.h>
#include "bitregs.h"

union fp32_t
{
  uint32_t u;
  float f;
  struct
  {
    uint32_t mantissa:23;
    uint32_t exp:8;
    uint32_t sign:1;
  };
};

union uint32_bytes
{
  uint8_t a[4];
  uint32_t u;
};

//multiplier primitive
uint16_t mul8x8(uint8_t a, uint8_t b) { return uint16_t(a)*uint16_t(b); }


uint24_t mul24x24_truncated(uint24_t a, uint24_t b)
{
   //(a3+a2+a1)*(b3+b2+b1) TRUNCATED: (a3*b3+a3*b2+a3*b1 + a2*b3+a2*b2 + a1*b3)
   uint32_bytes x = {.u = a };
   uint32_bytes y = {.u = b };

   uint16_t a3b3 = mul8x8(x.a[2], y.a[2]);
   uint16_t a3b2 = mul8x8(x.a[2], y.a[1]);
   uint16_t a3b1 = mul8x8(x.a[2], y.a[0]);
   uint16_t a2b3 = mul8x8(x.a[1], y.a[2]);
   uint16_t a2b2 = mul8x8(x.a[1], y.a[1]);
   uint16_t a1b3 = mul8x8(x.a[0], y.a[2]);

   const uint8_t rounding = 0x80;
   uint10_t r0 = (a3b1 + a2b2 + a1b3 + rounding) >> 8;
   uint17_t r1 = a3b2 + a2b3;
   uint24_t r2 = a3b3 << 8;

   return r2+r1+r0;
}

//trivial implementation:
//uint24_t mul24x24_truncated(uint24_t a, uint24_t b) { return ((uint48_t(a) * uint48_t(b)) >> 24); }

uint25_t fixpoint_mul23x23(uint23_t a, uint23_t b)
{
  //(1+a)+(1+b) = 1+a+b+a*b, 0 <= inputs < 1
  uint24_t mul = mul24x24_truncated(uint24_t(a),uint24_t(b));
  return (1<<23) + a + b + (mul<<1);
}

fp32_t fp32_mul(fp32_t a, fp32_t b)
{
  fp32_t c;

  uint25_t m = fixpoint_mul23x23(a.mantissa, b.mantissa);
  uint1_t expadj = m >> 24;
  if(expadj)
    m = m >> 1;

  c.mantissa = m;
  c.exp = a.exp + b.exp + expadj - 127;
  c.sign = a.sign ^ b.sign;

  return c;
};



//////////////////////////
/////// TEST CODE ////////
//////////////////////////

void disp(const char *name, fp32_t val)
{
  printf("%s: 0x%08X %+f, unbiased exp=%d, mantissa=0x%06X %+f\n",
    name, val.u, val.f, val.exp-127, val.mantissa, 1+val.mantissa/float(1<<23));
};

fp32_t a = {.f = 1.2345 };
fp32_t b = {.f = 6.7890 };
//fp32_t a = {.f = 1.1 };
//fp32_t b = {.f = 1.1 };
//fp32_t a = {.f = 1.2345 };
//fp32_t b = {.f = 6.00789 };
//fp32_t a = {.f = 1.99 };
//fp32_t b = {.f = 1.88 };

int main()
{
  fp32_t c = fp32_mul(a, b);
  fp32_t r = {.f = a.f * b.f };

  disp("a  ", a);
  disp("b  ", b);
  disp("c  ", c);
  disp("ref", r);
  printf("DIFFERENCE %f, (LSB): %d\n", r.f-c.f, r.u-c.u);

  return 0;
}

[suarezvictor]

with other tests and corrections, precision within 1 LSB

//Calculates FP32 mult using 6 instances of 8x8 multipliers
//(C) 2021 Victor Suarez Rovere <[email protected]>
//LICENSE: GPLv3
//
//Compile & run:
//g++ -O3 -o fptest fp32mul.c && ./fptest

#include <stdio.h>
#include "bitregs.h"

#define mul24x24_truncated mul24x24_truncated8 //function to test, should be correct within +/-1 LSB

union fp32_t
{
  uint32_t u;
  float f;
  struct
  {
    uint32_t mantissa:23;
    uint32_t exp:8;
    uint32_t sign:1;
  };
  struct {
    uint32_t a11_00:12;
    uint32_t a23_12:12;
  };
  struct {
    uint32_t a08_00:9;
    uint32_t a17_09:9;
    uint32_t a26_18:9;
  };
  struct {
    uint32_t a07_00:8;
    uint32_t a15_08:8;
    uint32_t a23_16:8;
  };
};



//multiplier primitive
uint16_t mul8x8(uint8_t a, uint8_t b) { return a*b; }
//multiplier-add primitive
uint18_t mul8x8add(uint8_t a, uint8_t b, uint18_t c) { return mul8x8(a,b)+c; }

uint24_t mul24x24_truncated8(uint24_t a, uint24_t b)
{
   //(a2+a1+a0)*(b2+b1+b0) TRUNCATED: (a2*b2 + a2*b1+a1*b2 + a2*b0+a1*b1+a0*b2)
   fp32_t x = {.u = a };
   fp32_t y = {.u = b };

   const uint8_t rounding = 1<<7;
   uint8_t r0a = (x.a15_08 >> 4)*(y.a07_00 >> 4) & 0xF8; //mult4x4 to 5 msb
   uint8_t r0b = (x.a07_00 >> 4)*(y.a15_08 >> 4) & 0xF8; //mult4x4 to 5 msb

   uint18_t r11 = mul8x8add(x.a15_08, y.a15_08, r0a + r0b + rounding); //achieves 1 LSB precision
   uint18_t r02 = mul8x8add(x.a07_00, y.a23_16, r11);
   uint18_t r20 = mul8x8add(x.a23_16, y.a07_00, r02);
   uint16_t r21 = mul8x8add(x.a23_16, y.a15_08, r20 >> 8);
   uint18_t r12 = mul8x8add(x.a15_08, y.a23_16, r21);
   uint16_t r2 = mul8x8(x.a23_16, y.a23_16);

   return (r2 << 8) + r12;
}

uint24_t mul24x24_truncated24(uint24_t a, uint24_t b)
{ 
  return ((uint48_t(a) * uint48_t(b)+(1<<23)) >> 24);
}

uint25_t fixpoint_mul23x23(uint23_t a, uint23_t b)
{
  //(1+a)+(1+b) = 1+a+b+a*b, 0 <= inputs < 1
  uint24_t mul = mul24x24_truncated(uint24_t(a),uint24_t(b));
  return (1<<23) + a + b + (mul<<1);
}

fp32_t fp32_mul(fp32_t a, fp32_t b)
{
  fp32_t c;

  uint25_t m = fixpoint_mul23x23(a.mantissa, b.mantissa);
  uint16_t expadj = uint32_t(m) >> 24;
  if(expadj)
    m = m >> 1;

  uint10_t exp = a.exp + b.exp + expadj;
  c.exp = exp - 127;
  c.sign = a.sign ^ b.sign;
  c.mantissa = m;

  return c;
};



//////////////////////////
/////// TEST CODE ////////
//////////////////////////

void disp(const char *name, fp32_t val)
{
  printf("%s: 0x%08X %+g,\tunbiased exp=%+d,\tmantissa=0x%06X %+g\n",
    name, val.u, val.f, val.exp-127, val.mantissa, 1+val.mantissa/float(1<<23));
};

#include <stdlib.h> //rand(), srand()


int main()
{


  //srand(0);

  int maxerror = 0, minerror = 0;
  for(int i=0; i<10*1000*1000; ++i)
  {
    fp32_t a, b, r;

    a.u = rand();
    b.u = rand();
    fp32_t c = fp32_mul(a, b);
    r.f = a.f * b.f; //C compiler function

    bool normal = (c.exp < 127*2 && r.exp < 127*2 && a.exp != 0 && b.exp != 0 && r.exp != 0 && c.exp != 0);
    if(!normal)
      continue;

    int lsbdiff = r.u-c.u;

    const int tol = 1;
    {
      if(c.sign != r.sign)
       exit(1);

      if(lsbdiff < minerror) minerror = lsbdiff;
      if(lsbdiff > maxerror) maxerror = lsbdiff;
    }
  }
  printf("Error: min %d, max %d\n", minerror, maxerror);

  return 0;
}

[suarezvictor]

Same for 12x12

//multiplier primitive
uint24_t mul12x12(uint12_t a, uint12_t b) { return (uint32_t)a*(uint32_t)b; }
//multiplier-add primitive
uint25_t mul12x12add(uint12_t a, uint12_t b, uint25_t c) { return (uint32_t)mul12x12(a,b)+(uint32_t)c; }

uint24_t mul24x24_truncated12(uint24_t a, uint24_t b)
{
   //(a1+a0)*(b1+b0) TRUNCATED: (a1*b1 + a1*b0+a0*b1)
   fp32_t x = {.u = a };
   fp32_t y = {.u = b };

   const uint12_t rounding = 1<<11;
   uint25_t a0b0 = (mul12x12add(x.a11_00>>9, y.a11_00>>9, 0)<<6) & 0xF00; //3x3 mult, 4 bit MSB result
   uint25_t a0b1 = mul12x12add(x.a11_00, y.a23_12, a0b0 + (uint25_t) rounding);
   uint25_t a1b0 = mul12x12add(x.a23_12, y.a11_00, a0b1);
   uint25_t a1b1 = mul12x12add(x.a23_12, y.a23_12, a1b0>>12);
   return (uint24_t)a1b1;
}

[suarezvictor]

same for 9x9 mults

//multiplier primitive
uint18_t mul9x9(uint9_t a, uint9_t b) { return (uint16_t)a*(uint16_t)b; }
//multiplier-add primitive
uint20_t mul9x9add(uint9_t a, uint9_t b, uint20_t c) { return mul9x9(a,b)+c; }

uint24_t mul24x24_truncated9(uint24_t a, uint24_t b)
{
   //(a2+a1+a0)*(b2+b1+b0) TRUNCATED: (a2*b2 + a2*b1+a1*b2 + a2*b0+a1*b1+a0*b2)
   fp32_t x = {.u = a << 3 };
   fp32_t y = {.u = b << 3 };

   const uint12_t rounding = (1<<8)<<3;

   uint20_t r11 = mul9x9add(x.a17_09, y.a17_09, (uint20_t)rounding);
   uint20_t r02 = mul9x9add(x.a08_00, y.a26_18, r11);
   uint20_t r20 = mul9x9add(x.a26_18, y.a08_00, r02);
   uint20_t r21 = mul9x9add(x.a26_18, y.a17_09, r20 >> 9);
   uint20_t r12 = mul9x9add(x.a17_09, y.a26_18, r21);
   uint18_t r2 = mul9x9(x.a26_18, y.a26_18);

   return ((r2 << 9) + r12) >> 3;
}