From: Tom Marshall
       To: All
  Subject: Fixed point math

I have been doing a fair amount of programming with fixed point math
recently.  Fixed addition and subtraction are fairly obvious operations.
Fixed multiplication isn't too bad.  But fixed division always makes me
stop and think for a minute or scribble on a piece of scratch paper.
So, I decided to make a few macros for doing fixed multiplication and
division.  Maybe someone reading this will find it useful.  Here's the
concept:

First off, this discussion will assume that we are using at least a 386
and all numbers are stored in 32-bit registers in "16.16" format.  That
is, the high 16 bits represent the integer portion and the low 16 bits
represent the fractional portion.  Here's a diagram:

           EAX (hi)         EAX (lo)
      -----------------------------------
Bit : | 31 30 .. 17 16 | 15 14 ..  1  0 |
      -----------------------------------
        |        |  |    |  |
        |        |  |    |   \-- 2^-2 (ie. quarters) place
        |        |  |     \----- 2^-1 (ie. halves) place
        |        |   \---------- 2^0  (ie. ones) place
        |         \------------- 2^1  (ie. twos) place
         \---------------------- 2^15 place or sign bit

If we call our number N, this representation can be thought of as
N*2^16.  The "decimal point" is between bit 16 and 15.  We can represent
numbers in the range of +/- 65536 for unsigned or +/- 32768 for signed.

Okay, now how do we manipulate these numbers?

Adding and subtracting is really simple.  Just use the normal ADD and
SUB instructions.

When multiplying, you must remember that multiplying two 32- bit
registers leaves a 64-bit result.  Using a little 8th grade math gives
us this:

(A*2^16) * (B*2^16) = A*B*2^32

The result is a 64-bit number with 32 bits in the integer portion and 32
bits in the fractional portion.  Notice that we can't use those nifty
new forms of IMUL for doing this!  The integer portion is in the high
register (EDX) and the fractional portion is in the low register (EAX).
To convert back to our original form, we shift EDX:EAX right by 16 bits
and disregard EDX.  The result must be in the range of +/- 65536 or we
will lose integer precision converting it. Fraction precision will
always be lost b ecause we only have room for 16 fraction bits.  But
that's the price for super fast fixed math.

Okay, now for the fun part:  division.  Using our mathematical
representation again, we get this:

(A*2^16) / (B*2^16) = A/B

That's no good because we're left with an integer!  We need to multiply
the numerator by 2^16 before we divide to ensure that the answer has 16
bits of fractional precision. So let's try again:

(A*2^16*2^16) / (B*2^16) = A/B*2^16

Ok, that's better.  The result is in exactly the form we want!  Now what
about the remainder?  It's not good for much at all because we're not
dividing integers.  Don't try to use it to round your result .. it won't
work!

Now for the useful part... macros!  These are all in TASM format, so
switch the macro declarations to use them under MASM (ie. MACRO foo =>
foo MACRO).  That's it for now... someone please let me know if you
found this stuff useful. I can explain a sine table or whatever next.

MACRO           FIXMUL                  ;EAX * EBX => EAX (EDX used)
                xor     edx,edx         ;  All numbers unsigned 16.16
                mul     ebx             ;  Overflow not checked
                shrd    eax,edx,16
ENDM

MACRO           IFIXMUL                 ;EAX * EBX => EAX (EDX used)
                cdq                     ;  All numbers signed 16.16
                imul    ebx             ;  Overflow not checked
                shrd    eax,edx,16
ENDM

MACRO           FIXDIV                  ;EAX / EBX => EAX (EDX = Rem)
                xor     edx,edx         ;  All numbers unsigned 16.16
                ror     eax,16          ;  Overflow / fault not checked
                xchg    ax,dx
                div     ebx
ENDM

MACRO           IFIXDIV                 ;EAX / EBX => EAX (EDX = Rem)
                cdq                     ;  All numbers signed 16.16
                ror     eax,16          ;  Overflow / fault not checked
                xchg    ax,dx
                idiv    ebx
ENDM