Fixed-point unsigned division in C

Question

I need an algorithm to do unsigned fixed-point division in C. I can use at most 32-bit words.

I want to minimize the number of bits needed to represent the integer part while being able to use numbers in the range [0..15]. So apparently the minimum number of bits is 4. Problem is the algorithm which I came up only works using 5 bits. Because it compares the remainder with the divisor and then shifts the remainder until it is bigger than the divisor, if the divisor has the most significant bit 1, then the algorithm will do nothing but shift the remainder (it will never be bigger). Here's the code:

int divu(int a, int b){
   int pt_int, r, pt_frac=0;
   int i;

   pt_int = ((unsigned) a/b) << BITS_FRAC;
   r = (unsigned) a%b;

   for (i=BITS_FRAC; i>=0; i--){
      if ((unsigned) r < b)
         r <<= 1;
      else{
         r -= b;
        pt_frac += 01 << i;
        r <<= 1;
      }
   }
   return pt_int + pt_frac;
}

If you do have a solution but don't want to understand the code, please, just post it. :)

Example:

We want to divide 1.5 by 2, which results 0.75. Suppose we are using 4 bits for the integer part and 28 bits for the fraction. So our numbers is hex are:

1.5:    0x18000000
2:      0x20000000
result: 0x0c000000 

Solution 1

You have a 4.28 fixed point number, and you want to divide by a 4.28 number. You find the precision after division by subtracting the precision of the numerator from the denominator, so a straight divide would give 4.28 - 4.28 = 0 -- no significant bits. Obviously this won't work.

     1.5  [4.28] 0x18000000 
   / 2.0  [4.28] 0x20000000 
   =  0?  [0.00] 0x00000000

The ideal way to do it would be to promote the numerator to 8.56 (by multiplying by 2^28), and then doing a 64 bit divide:

                     .   .   .   .
     1.5  [8.56] 0x180000000000000
   / 2.0  [4.28]        0x20000000 
   = 0.75 [4.28]        0x0c000000

If you truly can't use 64 bit numbers then your only option is to reduce the denominator. For example you can use half the precision by dividing by 2^14

     1.5  [4.28] 0x18000000 
   / 2.0  [2.14]     0x8000
   = 0.75 [2.14]     0x3000

You can then multiply the result by the same factor to get back to a 4.28 number: 0x3000 *(1<<14) = 0x0c000000

You do lose some precision this way, but it's unavoidable without using larger numerators. For example 5.0/3.0 = 1.66667 = 0x1AAAAAA [4.28], but
((5.0<<28)/(3<<14))<<14 = 0x1AAA8000 [4.28] = 1.66662

Solution 2

As pointed out here (Example: Multiplying an Integer by the Reciprocal of a Constant), you could re-implement your division by multiplication with the reciprocal. Afterwards you should be able to represent the integer part with 4 Bits.

Comments

[Kei Nivky]

I guess it's better to do my way. The objective of not using 5 bits for the integer part is to avoid precision loss, the loss with your algorithm is greater than the 1 bit of mine.

[Sunilsingh]

Yours is more accurate, it's also quite slow. The usual motivation for fixed point is speed on a processor with slow float math. But if you want accuracy, your method is best. Losing at least one bit is unavoidable, I think.

[EML]

Does your link carry out the inverse and multiplication using full precision? In that case, this won't work - your method rounds the 'correct' answer, instead of carrying out fixed-point division. The two methods will give different results in general.