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 

Answer 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