🚀 go-pugleaf

Some while ago, I wrote a 96 bit by 64 bit signed divide.  I probably can't
post that code fragment - and it wouldn't do much good because it was in
68K assembly.  However, it just occurred to me that I could post a C version
of the algorithm which the original requestor could adapt and translate as
needed.  It is at the end of this message.

The algorithm is a non-restoring divide using just subtracts, adds, and
shifts.  It operates directly on signed numbers without first doing an
absolute value of the inputs.  It is much slower than the methods using
the FPU or integer multiplies already posted, but it can easily be
adapted to give exact answers for large multi-precision integers.  The
implementation uses carry out or sign bits to reduce access to other
parts of the inputs.  It codes very nicely on a CPU with a C flag, but
doesn't require such a flag (as this C version demonstrates.)  (The
loops are also designed to use the 68K's DBcc instruction.  They might
be adaptable to the PPC's COUNT register.)  To simplify the C code, this
version only works on integers as big as what the compiler supports.
(To handle bigger numbers than the / operator, you need to rewrite it...)

Most places that denomShift appears could simply start a multiprecision
add (subtract) offset by one "word".  Checking for divide by zero and
overflow is also needed.  (The original use for this algorithm didn't
need these checks - intersecting two line segments could never overflow.)

Ludis Langens
ludis@netcom.com


#include <StdIO.h>

/*
 *  For more about this class of algorithm, see chapter 5 of:
 *  Introduction to Arithmetic for Digital Systems Designers
 *  by Shlomo Waser and Michael J. Flynn
 *  ISBN 0-03-0605710-7
 */

#define QTYPE short
#define QSIZE 16
#define NTYPE long
#define NMIN (NTYPE)0xc0000000
#define NMAX (NTYPE)0x3fffffff
#define DTYPE short
#define DMIN (DTYPE)0xc000
#define DMAX (DTYPE)0x3fff

QTYPE MultiPrecisionDivide(NTYPE numer, DTYPE denom)
    {
    QTYPE quotient;
    NTYPE denomShift;
    int count;

    quotient = 0;
    count = QSIZE;
    while ((DMIN <= denom) && (denom <= DMAX) && count++)
        denom <<= 1;
    while ((NMIN <= numer) && (numer <= NMAX) && count--)
        numer <<= 1;
    if (count < 0)
        goto returnQuotient;

    denomShift = denom, denomShift <<= QSIZE;
    if (numer >= 0)
        {
        if (denom >= 0)
            {
            goto pp;
            do
                {
                numer <<= 1;
              pp:
                if ((numer -= denomShift) < 0) goto ppCC;
              ppCS:
                /* !Borrow == Carry == true */
                quotient <<= 1, quotient += 1;
                }
            while (count--);
            goto returnQuotient;
            do
                {
                numer <<= 1;
                if ((numer += denomShift) >= 0) goto ppCS;
              ppCC:
                /* !Borrow == Carry == false */
                quotient <<= 1;
                }
            while (count--);
            goto returnQuotient;
            }
        else
            {
            goto pm;
            do
                {
                numer <<= 1;
              pm:
                if ((numer += denomShift) < 0) goto pmCC;
              pmCS:
                /* !Borrow == Carry == true */
                quotient <<= 1, quotient -= 1;
                }
            while (count--);
            goto returnQuotient;
            do
                {
                numer <<= 1;
                if ((numer -= denomShift) >= 0) goto pmCS;
              pmCC:
                /* !Borrow == Carry == false */
                quotient <<= 1;
                }
            while (count--);
            goto returnQuotient;
            }
        }
    else
        {
        if (denom >= 0)
            {
            goto mp;
            do
                {
                numer <<= 1;
              mp:
                if ((numer += denomShift) >= 0) goto mpCS;
              mpCC:
                /* !Borrow == Carry == false */
                quotient <<= 1, quotient -= 1;
                }
            while (count--);
            goto returnQuotient;
            do
                {
                numer <<= 1;
                if ((numer -= denomShift) < 0) goto mpCC;
              mpCS:
                if (numer == 0) goto mpEQ;
                /* !Borrow == Carry == true */
                quotient <<= 1;
                }
            while (count--);
            goto returnQuotient;
          mpEQ:
            quotient <<= 1, quotient -= 1, quotient <<= count;
            }
        else
            {
            goto mm;
            do
                {
                numer <<= 1;
              mm:
                if ((numer -= denomShift) >= 0) goto mmCS;
              mmCC:
                /* !Borrow == Carry == false */
                quotient <<= 1, quotient += 1;
                }
            while (count--);
            goto returnQuotient;
            do
                {
                numer <<= 1;
                if ((numer += denomShift) < 0) goto mmCC;
              mmCS:
                if (numer == 0) goto mmEQ;
                /* !Borrow == Carry == true */
                quotient <<= 1;
                }
            while (count--);
            goto returnQuotient;
          mmEQ:
            quotient <<= 1, quotient += 1, quotient <<= count;
            }
        }

  returnQuotient:
    return quotient;
    }

void TestMultiPrecisionDivide(NTYPE numer, DTYPE denom)
    {
    printf("%d / %d = %d (should be %d)\n", (int)numer, (int)denom,
            (int)MultiPrecisionDivide(numer, denom), (int)(numer / denom));
    }

int main(void)
    {
    TestMultiPrecisionDivide(255, 15);
    TestMultiPrecisionDivide(255, -15);
    TestMultiPrecisionDivide(-255, 15);
    TestMultiPrecisionDivide(-255, -15);
    TestMultiPrecisionDivide(240, 15);
    TestMultiPrecisionDivide(240, -15);
    TestMultiPrecisionDivide(-240, 15);
    TestMultiPrecisionDivide(-240, -15);
    TestMultiPrecisionDivide(31000, 32000);
    TestMultiPrecisionDivide(-3, 2);
    TestMultiPrecisionDivide(0x3456789a, 0x789a);
    return 0;
    }