fixed point integer division ("fractional division") algorithm

Question

The Honeywell DPS8 computer (and others) have/had a "divide fractional" instruction:

"This instruction divides a 71-bit fractional dividend (including sign) by a 36-bit fractional divisor (including sign) to form a 36-bit fractional quotient (including sign) and a 36-bit fractional remainder (including sign). Bit 35 of the remainder corresponds to bit 70 of the dividend. The remainder sign is equal to the dividend sign unless the remainder is zero."

So, as I understand it, this is integer division with the decimal point way over on the left.

  .qqqqq / .ddddd

(I did scaled integer math in FORTH back in the day, but my memories of the techniques are lost in fog of time.)

To implement this instruction in a DPS8 emulator, I believe I need to start by creating two 70 bit numbers: the 71 bit dividend less it's sign bit, and the the 36 bit divisor less its sign bit and shifted 35 bits to the left so that the decimal points line up.

I think I can then form the remainder and quotient (in C) with '%' and '/', but I am unsure if those results need to be normalized (i.e. shifted).

I found an example of a "shift and subtract" algorithm "Computer Arithmetic", slide 10), but I would prefer a more straight forward implementation.

Am I on the right track, or is the solution more nuanced (fixing up the signs and detection of errors have been elided from here; those stages are well documented. The actual division is the issue.). Any pointers to C implementations of this kind of hardware emulation would be particularly helpful.

Answer 1

I do not have the definitive answer, but as a division is a division, you might find it helpful to look at some basic division routines.

Imagine that you have a 32-bit variable and you want an 8-bit fractional part. You then have an integer part between 0 and 16777215, and a fractional part which is between 0 and 255. 0xiiiiiiff (where i is the integer part, f is the fractional part).

Imagine you have a 24-bit dividend (numerator), say the value 3, and a 24-bit divisor (denominator), say the value 13. As we quickly will see, 3/13 is greater than zero and less than one. That means our fractional part is nonzero, but our integer part is filled completely with zeros.

So to do the above division using a standard divide function, we'll just bit-shift the dividend by N, thus we will get N bits of precision in our fractional part.

quotient_fp = (dividend_ip << 8) / divisor_ip

So far, so good.

But what if we want the divisor to have a fractional part, then ?

If we just shift the divisor up by 8, then we'll have a problem: (dividend_ip << 8) / (divisor_ip << 8) - because we'll obviously lose our fractional part of the quotient (result).

Instead, we'll need to shift the dividend up by as many bits as we shift the fractional part up...

((dividend_ip << 8) << 8) / (divisor_ip << 8)

...That makes it... (dividend_ip << (dividend_precision + divisor_precision) / (divisor_ip << divisor_precision)

Now, let's put our fractional part math into the picture...

(((dividend_ip << dividend_precision) | dividend_fp) << divisor_precision) / ((divisor_ip << divisor_precision) | divisor_fp)

Our quotient's precision will be the same as dividend_precision, which is 8 bits.

Unfortunately, this eats a lot of bits.

Fortunately, in your case, the integer part is not important, so you'll have a lot of room for the fractional part. Let's increase the precision to 15 bits; this can be tested using normal 32-bit integers...

(((dividend_ip << 15) | dividend_fp) << 15) / ((divisor_ip << 15) | divisor_fp)

Our quotient will now have a 15-bit precision.

OK, but since you're supplying only the fractional parts and the integer part is always zero anyway, you should be able to just toss the integer part. That makes it.... (((dividend_ip << 16) | dividend_fp) << 16) / ((divisor_ip << 16) | divisor_fp) ... reduced to ... (dividend_fp << 16) / divisor_fp

... now let's use a 64-bit integer instead, we can get 32 bits of precision in the quotient... (dividend_fp << 32) / divisor_fp

... some compilers have support for a int128_t (it can be enabled on some platforms for GCC), so you might be able to use that type, in order to get 128 bits easily. I have not tried it, but I've come across info on the Web earlier; search for int128_t, and you might find out how.

If you get the int128_t to work, you could make the dividend 128 bit, the divisor 64 bit and the quotient 64 bit... quotient_fp = ((dividend_fp << 36) / divisor) >> (64 - 36) ... in order to get 36 bits precision. Notice that since the result is in the top 36 bits of the quotient, the quotient needs to be shifted down (64 - 36) = 28 bits. You could even go as high as (128 - 36) = 92 bits precision: (dividend_fp << 92) / divisor

Now, that you probably (hopefully) have a solution, I would like to recommend that you get familiar with low-level binary divide (again; since you've been there a while ago). The best sources seem to be how hardware divides binary numbers; such as microcontrollers, CPUs and the like. Assembly language dividers are also good for getting to know the inner workings. Often 32-bit divide routines that use bit-shifting are very good sources.

Through the time, I've come across a very clever implementation for ARM in ARM assembly language. Normally I wouldn't post references or assembly language examples, but considering that the code is very small, I think it would be alright.

Taken from A Fast Hi Precision Fixed Point Divide

r0 is the numerator (dividend) r2 is the denominator (divisor)

    mov     r1,#0
    adds    r0,r0,r0
    .rept   32
    adcs    r1,r2,r1,lsl#1
    subcc   r1,r1,r2
    adcs    r0,r0,r0
    .endr

r0 is the quotient (result) r1 is the remainder (rest, modulo result)

The above routine contains the basics for an unsigned divide.

I hope this information will be useful. It may contain errors, as I have not tested any code or example mentioned. I'm confident, though, that it's not all wrong. ;)