3

I am wondering how to implement IEEE-754 32-bit single precision floating point division in binary with no division hardware and no floating point hardware?

I have shifting hardware, add, subtract, and multiply.

I have already implemented floating point multiplication, addition, and subtraction using 16-bit words.

I am implementing these instructions on a proprietary multicore processor and writing my code in assembly. Beforehand, I am using matlab to verify my algorithm.

I know I need to subtract the exponents, but how do i perform unsigned division on the mantissas?

11
  • 1
    You haven't included any example code that represents what you mean and the language you're using. They are needed to answer your query.
    – Lion
    Jan 25, 2012 at 22:25
  • 1
    It is only too much fluff if you're dumping a thousand-line stack trace or an entire program.
    – BoltClock
    Jan 25, 2012 at 22:29
  • @BoltClock Less is more sometimes.
    – Veridian
    Jan 25, 2012 at 22:32
  • Not pleased but vote it up as you say to make you pleased.
    – Lion
    Jan 25, 2012 at 22:36
  • 2
    There is nothing wrong with what you are doing. No reason why you cant get the precision you are looking for without floating point hardware. A good floating point hardware verification platform written against no floating point hardware is testfloat which uses softfloat jhauser.us/arithmetic/SoftFloat.html as the reference. Other resources are divms.uiowa.edu/~jones/bcd/divide.html if you know ahead of time what the divisor is, and hackersdelight.org and others like it for generic division routines.
    – old_timer
    Jan 26, 2012 at 4:36

2 Answers 2

3

Depends on how complicated you want to make it. Keeping it reasonably simple, you could try division by reciprocal approximation.

Rather than calculating: (n / d) you'd work out: n * (1 / d).

To do this you'd need to work out the reciprocal using some method, for example, Newton-Raphson which uses Newton's method to calculate successively more accurate estimates of the reciprocal of the divisor until it's "adequately" accurate for your purpose before doing the final multiplication step.

EDIT

Just seen your update. This may, or may not, be useful for you after all!

0
3

If your hardware has a sufficiently fast integer multiplier (say 4-5 cycles), using an iterative approach to compute recip = 1 / divisor is likely to be the fastest approach. You would then compute the quotient as dividend * recip. It is very helpful for the necessary fixed-point computations if the hardware offers either integer multiplication with double-width result (i.e. the full product) or a "mulhi" instruction that delivers the upper 32 bits of the product of two 32 bit integers.

You will need to re-scale the fixed-point computations on the fly to retain close to 32 bits in intermediate results to derive a result that is accurate to the 24 bits that are required for the final result after rounding.

The fastest approach is likely generating a 9-bit starting approximation "approx" to 1 / divisor followed by a cubically convergent iteration for the reciprocal:

e = 1 - divisor * approx
e = e * e + e
recip = e * approx + approx

It is easiest to precompute the starting approximation and sstore it in an array of 256 bytes, indexed by bits 23:16 of the divisor (i.e. the 8 most significant fractional bits of the mantissa). As all approximation values are numbers in the range 0x100 ... 0x1FF (corresponding to 0.5 to 0.998046875) there is no need to store the most significant bit of each value as it will be '1'. Simply add 0x100 to the table element retrieved to get the final value of the initial approximation.

If you cannot afford 256 bytes of table storage, an alternative way to generate a starting approximation accurate to 9 bits would be a polynomial of degree 3 that approximates 1 / (1+f) where f is the fractional part of the divisor mantissa, m. Since with IEEE-754, m is known to be in [1.0,2.0), f is in [0.0,1.0).

Correct rounding to nearest-or-even (if required) could be implemented by the back-multiplication of the preliminary quotient by the divisor to establish the remainder, and selecting the final quotient such that the remainder is minimized.

The following code demonstrates the approximation principles discussed above, using the simpler case of the reciprocal, with proper rounding according to IEEE-754 nearest-or-even mode, and with appropriate handling of special cases (zero, denormals, infinity, NaNs). It assumes that a 32-bit IEEE-754 single-precision float can be transferred bit-wise from and to a 32-bit unsigned int. All operations are then performed on 32-bit integers.

unsigned int frcp_rn_core (unsigned int z)
{
  unsigned int x, y;
  int sign;
  int expo;

  sign = z & 0x80000000;
  expo = (z >> 23) & 0xff;
  x = expo - 1;
  if (x > 0xfd) {
    /* handle special cases */
    x = z << 1;
    /* zero or small denormal returns infinity of like sign */
    if (x <= 0x00400000) {
      return sign | 0x7f800000;
    }
    /* infinity returns zero of like sign */
    else if (x == 0xff000000) {
      return sign;
    }
    /* convert SNaNs to QNaNs */
    else if (x > 0xff000000) {
      return z | 0x00400000;
    } 
    /* large denormal, normalize it */      
    else {
      y = x < 0x00800000;
      z = x << y;
      expo = expo - y;
    }
  }
  y = z & 0x007fffff;
#if USE_TABLE
  z = 256 + rcp_tab[y >> 15];  /* approx */
#else
  x = y << 3;
  z =                    0xe39ad7a0;
  z = mul_32_hi (x, z) + 0x0154bde4;
  z = mul_32_hi (x, z) + 0xfff87521;
  z = mul_32_hi (x, z) + 0x00001ffa;
  z = z >> 4;
#endif /* USE_TABLE */
  y = y | 0x00800000;
  /* cubically convergent approximation to reciprocal */
  x = (unsigned int)-(int)(y * z); /* x = 1 - arg * approx */
  x = mul_32_hi (x, x) + x;        /* x = x * x + x */
  z = z << 15;
  z = mul_32_hi (x, z) + z;        /* approx = x * approx + approx */
  /* compute result exponent */
  expo = 252 - expo;
  if (expo >= 0) {
    /* result is a normal */
    x = y * z + y;
    z = (expo << 23) + z;
    z = z | sign;
    /* round result */
    if ((int)x <= (int)(y >> 1)) {
      z++;
    }
    return z;
  } 
  /* result is a denormal */
  expo = -expo;
  z = z >> expo;
  x = y * z + y;
  z = z | sign;
  /* round result */
  if ((int)x <= (int)(y >> 1)) {
    z++;
  }
  return z;
}

The function mul_32_high() is a placeholder for a machine-specific operation that returns the upper 32 bits of a signed multiplication of two 32-bit bit integers. A semi-portable implementation in lieu of a machine-specific version is

/* 32-bit int, 64-bit long long int */
int mul_32_hi (int a, int b)
{
  return (int)(unsigned int)(((unsigned long long)(((long long)a)*b)) >> 32);
}

The 256-entry reciprocal table used by the table-based variant can be constructed as follows:

static unsigned char rcp_tab[256];  
for (int i = 0; i < 256; i++) {
  rcp_tab[i] = (unsigned char)(((1./(1.+((i+.5)/256.)))*512.+.5)-256.);
}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.