Import upstream version 5.3.

[mup] / mup / mup / rational.c

/* Copyright (c) 1995, 1996, 2001 by Arkkra Enterprises */
/* All rights reserved */
/*
 *      rational.c      functions to do operations on rational numbers
 *
 *      Contents:
 *
 *              radd(), rsub(), rmul(), rdiv(), rneg(), rinv(), rrai(), rred(),
 *              ator(), rtoa(); also gtrat(), called by macros GT,GE,LT,LE
 *
 *              ratmsg(), add64_64(), mul32_64(), divmod64(), red64_64()
 *
 *              The first group of functions are for the user.  The second
 *              are for internal use only.
 *
 *      Description:
 *              The functions in this file do operations on rational numbers.
 *              The rational arguments to functions that the users can call
 *              must be in standard form (lowest terms, with positive
 *              denominator) except for rred().  There are checks for division
 *              by zero and for overflow of numerators and denominators.
 *              (The absolute values of each are limited to MAXLONG, defined
 *              in rational.h.)  If there is an error, the external int
 *              raterrno is set to RATDIV0 or RATOVER, as the case may be,
 *              and *raterrfuncp is checked.  If nonzero, it is assumed to
 *              point at the user's error handler, and it is called with a
 *              parameter equal to raterrno.  Otherwise, a message is printed
 *              to stderr.  In any case, the answer returned to the user is
 *              0/1.  If there was no error, raterrno is set to RATNOERR.
 *
 *              In general, the functions assume they are being called with
 *              valid parameters.  If they are not, results are not guaranteed
 *              to be correct.  However, they are defensive enough so that
 *              invalid parameters will not cause a crash in these routines.
 *              They will not always detect invalid parameters, but if they
 *              do, they will use the raterrno/raterrfuncp mechanism described
 *              above, with the value RATPARM.
 *
 *              These routines depend on a INT32B being a 32-bit number,
 *              stored in two's complement form, and UINT32B being the same
 *              for unsigned.  See rational.h.  Numerators and denominators
 *              are assumed to be INT32B.  Furthermore, the number 0x80000000
 *              is not allowed.  The routines should work on any machine and
 *              compiler where these requirements are met.
 *
 *              Internally, when 64-bit numbers are used, they are represented
 *              by an array of two INT32B.  The 0 subscript contains the low
 *              order bits and the 1 subscript contains the high order bits.
 *              The numbers are usually used as two's complement signed
 *              integers, so the high bit of the 1 subscript is a sign bit.
 */

#ifndef stderr
#       include <stdio.h>
#endif

#ifndef isspace
#       include <ctype.h>
#endif

#ifndef _RATIONAL
#       include "rational.h"
#endif



/*
 * Define SMALL to be a number that uses less than half as many bits as
 * MAXLONG (15 as compared to 31).  Define a SMALLRAT as a rational number
 * whose numerator (absolute value) and denominator are in that range.
 * The denominator is assumed to be positive.
 */
#define SMALL           0x7fff
#define SMALLRAT(q)     ((q).n <= SMALL && (q).n >= -SMALL && (q).d <= SMALL)


/*
 * This macro checks whether a 64-bit integer is actually less than or
 * equal to MAXLONG in absolute value, 0x7fffffff.
 */
#define INT32(n)        ((n)[1] ==  0 && (n)[0] >= 0 ||         \
                         (n)[1] == -1 && (n)[0] < 0 && (n)[0] != 0x80000000)


/*
 * Return whether the first 64-bit number equals the second.
 * To be equal, both words must be equal.
 */
#define EQ64(x, y)      ( (x)[1] == (y)[1] && (x)[0] == (y)[0] )


/*
 * Return whether the first 64-bit number is greater than the second.
 * If the high order words are equal, use the low order words, unsigned;
 * otherwise, just use the high order words.
 */
#define GT64(x, y) (                                            \
        (x)[1] == (y)[1] ?                                      \
                (UINT32B)(x)[0] > (UINT32B)(y)[0]               \
        :                                                       \
                (x)[1] > (y)[1]                                 \
)


/*
 * Return whether the first 64-bit number is less than the second.
 * If the high order words are equal, use the low order words, unsigned;
 * otherwise, just use the high order words.
 */
#define LT64(x, y) (                                            \
        (x)[1] == (y)[1] ?                                      \
                (UINT32B)(x)[0] < (UINT32B)(y)[0]               \
        :                                                       \
                (x)[1] < (y)[1]                                 \
)


/*
 * Return whether the first *unsigned* 64-bit number is less than or equal to
 * the second.  If the high order words are equal, use the low order words;
 * otherwise, just use the high order words.
 */
#define LEU64(x, y) (                                           \
        (x)[1] == (y)[1] ?                                      \
                (UINT32B)(x)[0] <= (UINT32B)(y)[0]              \
        :                                                       \
                (UINT32B)(x)[1] <= (UINT32B)(y)[1]              \
)


/*
 * Negate a 64-bit number.
 */
#define NEG64(x)        {                                               \
        (x)[1] = ~(x)[1];       /* one's complement */                  \
        (x)[0] = -(x)[0];       /* two's complement */                  \
        if ((x)[0] == 0)        /* if "carry" must inc high word */     \
                (x)[1]++;                                               \
}


/*
 * Shift a 64-bit number left one bit as unsigned (not that it matters).
 */
#define SHL1U64(x)      {                                               \
        (x)[1] <<= 1;           /* shift high word */                   \
        if ((x)[0] < 0)         /* if high bit of low word is set */    \
                (x)[1]++;       /* shift it into the high word */       \
        (x)[0] <<= 1;           /* shift low word */                    \
}


/*
 * Shift a 64-bit number right one bit as unsigned.
 */
#define SHR1U64(x)      {                                               \
        (x)[0] = (UINT32B)(x)[0] >> 1;  /* shift low word */            \
        if ((x)[1] & 1)                 /* if low bit of high word set*/\
                (x)[0] |= 0x80000000;   /* shift it into low word */    \
        (x)[1] = (UINT32B)(x)[1] >> 1;  /* shift low word */            \
}



/* declare as static the functions that are only used internally */
#ifdef __STDC__
static void ratmsg(int code);
static void add64_64(INT32B a[], INT32B x[], INT32B y[]);
static void mul32_64(INT32B a[], INT32B x, INT32B y);
static void divmod64(INT32B x[], INT32B y[], INT32B q[], INT32B r[]);
static void red64_64(INT32B num[], INT32B den[]);
#else
static void ratmsg(), add64_64(), mul32_64(), divmod64(), red64_64();
#endif


int raterrno;                   /* set to error type upon return to user */
void (*raterrfuncp)();          /* error handler functions to be called */

static RATIONAL zero = {0,1};
\f
/*
 *      radd()          add two rational numbers
 *
 *      This function adds two rational numbers.  They must be in standard
 *      form.
 *
 *      Parameters:     x       the first number
 *                      y       the second number
 *
 *      Return value:   The sum (x + y), if it can be represented as a
 *                      RATIONAL, else 0/1.
 *
 *      Side effects:   If radd() succeeds, it sets raterrno to RATNOERR.
 *                      Otherwise, the numerator or denominator must have
 *                      overflowed, so it sets raterrno to RATOVER and
 *                      either prints a message or calls a user-supplied
 *                      error handler.
 */

RATIONAL
radd(x, y)

RATIONAL x, y;

{
        RATIONAL a;                     /* the answer */
        INT32B bign[2];                 /* 64-bit numerator */
        INT32B bigd[2];                 /* 64-bit denominator */
        INT32B bigt[2];                 /* temp storage */


        raterrno = RATNOERR;            /* no error yet */

        /*
         * If the numbers are small enough, do it the easy way, since there is
         * then no danger of overflow.
         */
        if (SMALLRAT(x) && SMALLRAT(y)) {
                a.n = x.n * y.d + x.d * y.n;
                a.d = x.d * y.d;
                rred(&a);               /* reduce to standard form */
                return(a);
        }

        /*
         * To avoid overflow during the calculations, use two INT32B to
         * hold numbers.
         */
        mul32_64(bign, x.n, y.d);       /* get first part of numerator */
        mul32_64(bigt, x.d, y.n);       /* get second part of numerator */
        add64_64(bign, bign, bigt);     /* add to get full numerator */
        mul32_64(bigd, x.d, y.d);       /* get denominator */
        red64_64(bign, bigd);           /* reduce */

        /* overflow if the result can't fit in a RATIONAL */
        if ( ! INT32(bign) || ! INT32(bigd) ) {
                ratmsg(RATOVER);        /* set raterrno, report error */
                return(zero);
        }

        a.n = bign[0];                  /* set answer */
        a.d = bigd[0];

        return(a);
}
\f
/*
 *      rsub()          subtract two rational numbers
 *
 *      This function subtracts two rational numbers.  They must be in standard
 *      form.
 *
 *      Parameters:     x       the first number
 *                      y       the second number
 *
 *      Return value:   The difference (x - y), if it can be represented as a
 *                      RATIONAL, else 0/1.
 *
 *      Side effects:   If rsub() succeeds, it sets raterrno to RATNOERR.
 *                      Otherwise, the numerator or denominator must have
 *                      overflowed, so it sets raterrno to RATOVER and
 *                      either prints a message or calls a user-supplied
 *                      error handler.
 */

RATIONAL
rsub(x, y)

RATIONAL x, y;

{
        /*
         * Just negate the second operand and add.  We could call rneg() to
         * negate y, but why waste the time?
         */
        y.n = -y.n;
        return(radd(x, y));
}
\f
/*
 *      rmul()          multiply two rational numbers
 *
 *      This function multiplies two rational numbers.  They must be in standard
 *      form.
 *
 *      Parameters:     x       the first number
 *                      y       the second number
 *
 *      Return value:   The product (x * y), if it can be represented as a
 *                      RATIONAL, else 0/1.
 *
 *      Side effects:   If rsub() succeeds, it sets raterrno to RATNOERR.
 *                      Otherwise, the numerator or denominator must have
 *                      overflowed, so it sets raterrno to RATOVER and
 *                      either prints a message or calls a user-supplied
 *                      error handler.
 */

RATIONAL
rmul(x, y)

RATIONAL x, y;

{
        RATIONAL a;                     /* the answer */
        INT32B bign[2];                 /* 64-bit numerator */
        INT32B bigd[2];                 /* 64-bit denominator */


        raterrno = RATNOERR;            /* no error yet */

        /*
         * If the numbers are small enough, do it the easy way, since there is
         * then no danger of overflow.
         */
        if (SMALLRAT(x) && SMALLRAT(y)) {
                a.n = x.n * y.n;
                a.d = x.d * y.d;
                rred(&a);               /* reduce to standard form */
                return(a);
        }

        /*
         * To avoid overflow during the calculations, use two INT32B to
         * hold numbers.
         */
        mul32_64(bign, x.n, y.n);       /* get numerator */
        mul32_64(bigd, x.d, y.d);       /* get denominator */
        red64_64(bign, bigd);           /* reduce */

        /* overflow if the result can't fit in a RATIONAL */
        if ( ! INT32(bign) || ! INT32(bigd) ) {
                ratmsg(RATOVER);        /* set raterrno, report error */
                return(zero);
        }

        a.n = bign[0];                  /* set answer */
        a.d = bigd[0];

        return(a);
}
\f
/*
 *      rdiv()          divide two rational numbers
 *
 *      This function divides two rational numbers.  They must be in standard
 *      form.
 *
 *      Parameters:     x       the first number
 *                      y       the second number
 *
 *      Return value:   The quotient (x / y), if it is defined and can be
 *                      represented as a RATIONAL, else 0/1.
 *
 *      Side effects:   If rdiv() succeeds, it sets raterrno to RATNOERR.
 *                      Otherwise, either the second number was zero or the
 *                      numerator or denominator overflowed.  In this case,
 *                      it sets raterrno to RATDIV0 or RATOVER, respectively,
 *                      and either prints a message or calls a user-supplied
 *                      error handler.
 */

RATIONAL
rdiv(x, y)

RATIONAL x, y;

{
        RATIONAL r;                     /* reciprocal of y */


        r = rinv(y);                    /* first find 1/y */

        if (raterrno != RATNOERR)       /* if y was 0, return failure now */
                return(zero);

        /*
         * Return  x * r.   Whether rmul() succeeds or fails, we still just want
         * to leave raterrno the same and return what rmul() returns.
         */
        return(rmul(x, r));
}
\f
/*
 *      rneg()          negate a rational number
 *
 *      This function negates a rational number.  It must be in standard form.
 *
 *      Parameters:     x       the number
 *
 *      Return value:   The negative (-x).
 *
 *      Side effects:   It sets raterrno to RATNOERR.
 */

RATIONAL
rneg(x)

RATIONAL x;

{
        raterrno = RATNOERR;            /* no errors are possible */

        x.n = -x.n;

        /* answer is already in standard form since x was */
        return(x);
}
\f
/*
 *      rinv()          invert a rational number
 *
 *      This function inverts a rational number.  It must be in standard form.
 *
 *      Parameters:     x       the number
 *
 *      Return value:   The reciprocal (1 / x), if it is defined, else 0/1.
 *
 *      Side effects:   If rinv() succeeds, it sets raterrno to RATNOERR.
 *                      Otherwise, the second number must have been zero,
 *                      so it sets raterrno to RATDIV0 and either prints a
 *                      message or calls a user-supplied error handler.
 */

RATIONAL
rinv(x)

RATIONAL x;

{
        RATIONAL a;                     /* the answer */


        /* check for division by 0 */
        if (ZE(x)) {
                ratmsg(RATDIV0);        /* set raterrno, report error */
                return(zero);
        }

        raterrno = RATNOERR;            /* no errors from here on */

        a.n = x.d;                      /* flip numerator and denominator */
        a.d = x.n;

        if (a.d < 0) {                  /* if x was negative, reverse signs */
                a.n = -a.n;
                a.d = -a.d;
        }

        return(a);
}
\f
/*
 *      rrai()          raise a rational number to an integral power
 *
 *      This function raises a rational number to an integral power.  The
 *      rational number must be in standard form.
 *
 *      Parameters:     x       the rational number
 *                      n       the power, an integer
 *
 *      Return value:   The result (x to the nth power), if it is defined and
 *                      can be represented as a RATIONAL, else 0/1.
 *
 *      Side effects:   If rrai() succeeds, it sets raterrno to RATNOERR.
 *                      Otherwise, either zero is being raised to a non-
 *                      positive power, or the numerator or denominator
 *                      overflowed.  In this case, it sets raterrno to
 *                      RATDIV0 or RATOVER, respectively, and either prints
 *                      a message or calls a user-supplied error handler.
 */

RATIONAL
rrai(x, n)

RATIONAL x;
register int n;

{
        static RATIONAL one = {1,1};

        RATIONAL a;                     /* the answer */
        register int i;                 /* loop counter */


        /* it is undefined to raise zero to a nonpositive power */
        if (ZE(x) && n <= 0) {
                ratmsg(RATDIV0);        /* set raterrno, report error */
                return(zero);
        }

        raterrno = RATNOERR;            /* no error yet */

        a = one;                        /* init to 1 */
        if (n >= 0) {
                for (i = 0; i < n; i++) {
                        a = rmul(a, x);         /* mul again by x */
                        if (raterrno != RATNOERR)
                                return(zero);
                }
        } else {
                for (i = 0; i > n; i--) {
                        a = rdiv(a, x);         /* div again by x */
                        if (raterrno != RATNOERR)
                                return(zero);
                }
        }

        return(a);
}
\f
/*
 *      rred()          reduce a rational number to standard form
 *
 *      This function puts a rational number into standard form; that is,
 *      numerator and denominator will be relatively prime and the denominator
 *      will be positive.  On input, they may be any integers whose absolute
 *      values do not exceed MAXLONG.
 *
 *      Parameters:     ap      pointer to the rational number
 *
 *      Return value:   None.
 *
 *      Side effects:   If ap->d is 0, the function sets raterrno to RATDIV0,
 *                      either prints a message or calls a user-supplied
 *                      error handler, and sets *ap to 0/1.  Otherwise, it
 *                      sets raterrno to RATNOERR and puts *ap in standard form.
 */

void
rred(ap)

register RATIONAL *ap;

{
        register INT32B b, c, r; /* temp variables for Euclidean algorithm */
        register int sign;      /* answer is pos (1) or neg (-1) */


        /*
         * Since the numerator and denominator can be anything <= MAXLONG,
         * we must guard against division by 0.
         */
        if (ap->d == 0) {
                ratmsg(RATDIV0);        /* set raterrno, report error */
                *ap = zero;
                return;
        }

        raterrno = RATNOERR;            /* no errors possible from here on */

        if (ap->n == 0) {               /* if so, answer is "0/1" */
                ap->d = 1;
                return;
        }

        /* now figure out sign of answer, and make n & d positive */
        sign = 1;                       /* init to positive */
        if (ap->n < 0) {                /* reverse if numerator neg */
                sign = -sign;
                ap->n = -(ap->n);
        }
        if (ap->d < 0) {                /* reverse if denominator neg */
                sign = -sign;
                ap->d = -(ap->d);
        }

        /* now check whether numerator or denominator are equal */
        if (ap->n == ap->d) {           /* if so, answer is +1 or -1 */
                ap->n = sign;
                ap->d = 1;
                return;
        }

        if (ap->n < ap->d) {            /* set so that c > b */
                c = ap->d;
                b = ap->n;
        } else {
                c = ap->n;
                b = ap->d;
        }

        /* use Euclidean Algorithm to find greatest common divisor of c & b */
        do {
                r = c % b;
                c = b;
                b = r;
        } while (r != 0);

        /* now c is the greatest common divisor */

        ap->n /= c;             /* divide out greatest common divisor */
        ap->d /= c;             /* divide out greatest common divisor */

        if (sign < 0)           /* put sign in if result should be negative */
                ap->n = -(ap->n);

        return;
}
\f
/*
 *      ator()          convert an ascii string to a rational number
 *
 *      This function takes an ascii string as input and interprets it as
 *      a rational number.  White space may precede the number, but the
 *      number may not contain white space.  The numerator may be preceded
 *      by a minus sign.  The denomintor is optional, but if present, must
 *      not contain a sign.  In short, the number must match one of the
 *      following lex regular expressions, which starts where s points and
 *      ends before the first character not matching the pattern:
 *              [ \t\n]*-?[0-9]+
 *              [ \t\n]*-?[0-9]+\/[0-9]+
 *      Further restrictions are that the absolute values of numerator and
 *      denominator cannot exceed MAXLONG, and the denominator cannot be 0.
 *      If neither pattern is matched, or the further restrictions are
 *      violated, the function sets *rp to 0/1 and returns NULL.  Otherwise,
 *      it sets *rp to the result in standard form, and returns a pointer to
 *      the first char after the number found.
 *
 *      Parameters:     rp      pointer to where the answer goes
 *                      s       string containing ascii rational number
 *
 *      Return value:   If a valid rational number is found, the function
 *                      returns a pointer to the next char in the string
 *                      following the number.  Otherwise it returns NULL.
 *
 *      Side effects:   If ator() succeeds, it sets *rp to the result.
 *                      Otherwise, it sets it to 0/1.
 */

char *
ator(rp, s)

register RATIONAL *rp;
register char s[];

{
        register char *p;       /* point somewhere in s[] */
        int sign;               /* 1 means positive, -1 negative */


        /* skip by white space */
        for (p = s; isspace(*p); p++)
                ;

        /* init sign to positive; then reverse it if a dash is found */
        sign = 1;
        if (p[0] == '-') {
                sign = -1;
                p++;
        }

        /* fail if there are no digits */
        if ( ! isdigit(*p) ) {
                *rp = zero;
                return(NULL);
        }

        /*
         * Collect the numerator digits, and defend against overflow.
         */
        rp->n = 0;
        while ( isdigit(*p) ) {
                if (rp->n > MAXLONG / 10) {
                        *rp = zero;
                        return(NULL);
                }
                rp->n *= 10;
                if (rp->n > MAXLONG - (*p - '0')) {
                        *rp = zero;
                        return(NULL);
                }
                rp->n += *p++ - '0';
        }

        if (sign < 0)                   /* make negative if necessary */
                rp->n = -(rp->n);


        /*
         * If there is to be a denominator, collect its digits.  Otherwise,
         * set it to 1.  Defend against overflow.
         */
        if (*p == '/') {
                p++;
                if ( ! isdigit(*p) ) {  /* must be digit (no '-' allowed) */
                        *rp = zero;
                        return(NULL);
                }
                rp->d = 0;
                while ( isdigit(*p) ) {
                        if (rp->d > MAXLONG / 10) {
                                *rp = zero;
                                return(NULL);
                        }
                        rp->d *= 10;
                        if (rp->d > MAXLONG - (*p - '0')) {
                                *rp = zero;
                                return(NULL);
                        }
                        rp->d += *p++ - '0';
                }
                if (rp->d == 0) {       /* zero denominator is a failure */
                        *rp = zero;
                        return(NULL);
                }
        } else {
                rp->d = 1;              /* no denominator; assume 1 */
        }

        rred(rp);                       /* reduce the fraction */

        return(p);                      /* first char after the number */
}
\f
/*
 *      rtoa()          convert a rational number to an ascii string
 *
 *      This function takes a rational number as input converts it into
 *      an ascii string.  If the denominator is 1, it will not be printed.
 *      The number must be in standard form.
 *
 *      Parameters:     s       pointer to where the answer goes
 *                      rp      pointer to the rational number
 *
 *      Return value:   The function returns a pointer to the next char in
 *                      the string following the number.
 *
 *      Side effects:   The function sets s[] to the result.
 */

char *
rtoa(s, rp)

register char s[];
RATIONAL *rp;

{
        register INT32B num, den;       /* copy of num and den from *rp */
        register int i;                 /* index into t[] */
        char t[12];                     /* temp answer string */


        num = rp->n;                    /* copy num and den for efficiency */
        den = rp->d;

        if (num < 0) {                  /* if num is negative */
                *s++ = '-';             /* output minus sign */
                num = -num;             /* and make num positive */
        }

        i = 0;
        do {                            /* calc digits in reverse order */
                t[i++] = num % 10 + '0';
                num /= 10;
        } while (num > 0);              /* always loop at least once so 0="0"*/

        while (--i >= 0)                /* copy digits to answer string */
                *s++ = t[i];

        if (den != 1) {                 /* if a denominator is needed */
                *s++ = '/';             /* fraction bar */
                i = 0;
                do {                    /* calc digits in reverse order */
                        t[i++] = den % 10 + '0';
                        den /= 10;
                } while (den > 0);

                while (--i >= 0)        /* copy digits to answer string */
                        *s++ = t[i];
        }

        return(s);
}
\f
/*
 *      gtrat()         decide whether one rational is greater than another
 *
 *      This function decides whether its first parameter is greater than
 *      its second.  It is used by the macros GT, GE, LT, and LE.
 *      The numbers must be in standard form.  (Actually, all that is
 *      matters is that the denominators be positive.)
 *
 *      Parameters:     x       the first rational
 *                      y       the second rational
 *
 *      Return value:   1 if x > y, otherwise 0
 *
 *      Side effects:   none
 */

int
gtrat(x, y)

RATIONAL x, y;

{
        INT32B a[2];            /* temp holding areas for 64-bit numbers */
        INT32B b[2];


        /* if no overflow possible, cross-multiply and return truth value */
        /* note:  this depends on positive denominators */
        if (SMALLRAT(x) && SMALLRAT(y))
                return(x.n * y.d > x.d * y.n);

        /*
         * The numbers are too big; we have to do it the hard way to avoid
         * overflow.  Cross-multiply.  Note: this depends on positive
         * denominators.
         */
        mul32_64(a, x.n, y.d);
        mul32_64(b, x.d, y.n);

        return(GT64(a, b));
}
\f
/*
 *      ratmsg()        handle rational error of type "code"
 *
 *      This function sets raterrno.  Then calls the user's error handler,
 *      if there is one, or else prints a message to standard error.
 *
 *      Parameters:     code    the error code
 *
 *      Return value:   None.
 *
 *      Side effects:   raterrno is set; then either a message is printed
 *                      to standard error or the user's error handler is
 *                      called.
 */

static void
ratmsg(code)

int code;

{
        raterrno = code;                /* set global error flag */

        if (raterrfuncp == 0) {
                /* no user trap exists, so print message from here */
                switch (code) {
                case RATOVER:
                        (void)fputs("rational overflow\n", stderr);
                        break;

                case RATDIV0:
                        (void)fputs("rational division by zero\n", stderr);
                        break;

                case RATPARM:
                        (void)fputs("invalid number passed to rational number routine\n", stderr);
                        break;

                default:
                        (void)fputs("error in rational routines\n", stderr);
                        break;
                }
        } else {
                /* call user trap function to handle the error */
                (*raterrfuncp)(code);
        }
}
\f
/*
 *      add64_64()      add 64-bit numbers to get a 64-bit number
 *
 *      This function adds two 64-bit signed numbers to get a 64-bit
 *      signed number.  It is assumed that the result will not overflow.
 *      Any of the inputs may be the same arrays.
 *
 *      Parameters:     a       answer goes here
 *                      x       the first input
 *                      y       the second input
 *
 *      Return value:   none
 *
 *      Side effects:   a is set to the result.
 */

static void
add64_64(a, x, y)

INT32B a[2];
INT32B x[2];
INT32B y[2];

{
        INT32B t[2];                    /* temp storage */


        /* first add low and high parts separately */
        /* use temp storage in case a[] is the same array as x[] or y[] */
        t[0] = x[0] + y[0];
        t[1] = x[1] + y[1];

        /* figure out if the low part carries into the high part */
        if (x[0] < 0 && y[0] < 0) {             /* both high order bits set */
                t[1]++;                         /* must be a carry */
        } else if (x[0] < 0 || y[0] < 0) {      /* exactly one high bit set */
                if (t[0] >= 0)                  /* if result high bit clear */
                        t[1]++;                 /* must be a carry */
        }

        a[0] = t[0];                    /* copy results */
        a[1] = t[1];
}
\f
/*
 *      mul32_64()      multiply 32-bit numbers to get a 64-bit number
 *
 *      This function multiplies two 32-bit signed numbers to get a 64-bit
 *      signed number.  The numbers must not equal 0x80000000.  Overflow
 *      cannot occur.
 *
 *      Parameters:     a       answer goes here
 *                      x       the first 32-bit number
 *                      y       the second 32-bit number
 *
 *      Return value:   none
 *
 *      Side effects:   a is set to the result.
 */

static void
mul32_64(a, x, y)

INT32B a[2];
INT32B x;
INT32B y;

{
        INT32B t[2];                    /* temp storage for inner terms */
        INT32B xl, xh;                  /* low and high 16 bits of x */
        INT32B yl, yh;                  /* low and high 16 bits of y */
        int sign;                       /* sign of the result */


        /* make both numbers positive and determine the sign of the result */
        sign = 1;                               /* start at positive */
        if (x < 0) {
                x = -x;
                sign = -sign;
        }
        if (y < 0) {
                y = -y;
                sign = -sign;
        }

        /* break x and y into high and low pieces */
        xl = x & 0xffff;                /* 0 <= xl <= 0xffff */
        xh = x >> 16;                   /* 0 <= xh <= 0x7fff */
        yl = y & 0xffff;                /* 0 <= yl <= 0xffff */
        yh = y >> 16;                   /* 0 <= yh <= 0x7fff */

        /* multiply the outer parts */
        a[0] = xl * yl;                 /* 0 <= a[0] <= 0xfffe0001 */
        a[1] = xh * yh;                 /* 0 <= a[1] <= 0x3fff0001 */

        /* multiply the inner parts and break the result in two pieces */
        t[0] = xl * yh + xh * yl;       /* 0 <= t[0] <= 0xfffd0002 */
        t[1] = (UINT32B)t[0] >> 16;             /* 0 <= t[1] <= 0x0000fffd */
        t[0] <<= 16;                            /* 0 <= t[0] <= 0xffff0000 */

        /* add the two partial products */
        add64_64(a, a, t);              /* 0 <= a <= 0x3fffffff00000001 */

        /* if the answer is supposed to be negative, negate it */
        if (sign < 0)
                NEG64(a);
}
\f
/*
 *      divmod64()      find quotient and remainder of two 64-bit numbers
 *
 *      This function takes two 64-bit numbers and divides the first by the
 *      second, to get a quotient and remainder, both 64 bits.  It is assumed
 *      that the first number is nonnegative and the second number is positive.
 *      q and r must be different arrays.
 *
 *      Parameters:     x       first number (dividend)
 *                      y       second number (divisor)
 *                      q       quotient
 *                      r       remainder
 *
 *      Return value:   none
 *
 *      Side effects:   q and r are altered to be the results
 *                      x and y are not altered
 */

static void
divmod64(x, y, q, r)

INT32B x[2];
INT32B y[2];
INT32B q[2];
INT32B r[2];

{
        INT32B s[2];            /* temp storage for divisor */
        INT32B t[2];            /* temp storage for scratch */
        register int shift;     /* how far has y been shifted? */


        r[0] = x[0];            /* copy dividend to remainder place */
        r[1] = x[1];
        s[0] = y[0];            /* copy divisor to temp storage */
        s[1] = y[1];

        /* shift divisor left until greater than dividend */
        /* compare as unsigned so no problem if it gets shifted into sign bit */
        for (shift = 0; LEU64(s, r); shift++)
                SHL1U64(s);

        SHR1U64(s);             /* shift it back right one, so <= dividend */
        shift--;

        q[0] = 0;               /* start quotient at 0 */
        q[1] = 0;

        /*
         * Loop once for each bit shifted.
         */
        for ( ; shift >= 0; shift--) {
                /*
                 * If the current divisor does not exceed what's left of the
                 * dividend, subtract it from it, and record that by setting
                 * the low order bit of the current quotient.
                 */
                if ( ! GT64(s, r) ) {
                        t[0] = s[0];
                        t[1] = s[1];
                        NEG64(t);
                        add64_64(r, r, t);
                        q[0] |= 1;
                }

                /* shift quotient left and divisor right */
                SHL1U64(q);
                SHR1U64(s);
        }

        SHR1U64(q);             /* shift quotient right */
}
\f
/*
 *      red64_64()      reduce a 64 bit over 64 bit rational to lowest terms
 *
 *      This function takes two 64-bit numbers as numerator and denominator
 *      of a rational number, and reduces them to lowest terms, with the
 *      denominator positive.  If the user called this package correctly, the
 *      denominator cannot be zero, but to be defensive we check for that, and
 *      if it happens, set raterrno to RATPARM and either print a message or
 *      call a user-supplied error handler.
 *
 *      Parameters:     num     numerator
 *                      den     denominator
 *
 *      Return value:   none
 *
 *      Side effects:   num and den are altered to be the result
 */

static void
red64_64(num, den)

INT32B num[2];
INT32B den[2];

{
        INT32B b[2], c[2], r[2]; /* temp variables for Euclidean algorithm */
        INT32B junk[2];         /* placeholder for calling divmod64 */
        int sign;               /* answer is pos (1) or neg (-1) */


        if (den[1] == 0 && den[0] == 0) { /* if den == 0 */
                /*
                 * This is an error.  The user must have called a routine with
                 * an invalid number for us to get here, in fact a number with
                 * a zero denominator, since "den" here always was created as
                 * the product of denominators.  Report the error and return
                 * zero.
                 */
                num[1] = 0;             /* set num = 0 */
                num[0] = 0;
                den[1] = 0;             /* set den = 1 */
                den[0] = 1;
                ratmsg(RATPARM);        /* set raterrno, report error */
                return;
        }

        if (num[1] == 0 && num[0] == 0) { /* if num == 0 */
                den[1] = 0;             /* set den = 1; answer is 0/1 */
                den[0] = 1;
                return;
        }

        /* now figure out sign of answer, and make num & den positive */
        sign = 1;                       /* init to positive */
        if (num[1] < 0) {               /* if numerator neg */
                sign = -sign;           /* reverse the sign */
                NEG64(num);
        }
        if (den[1] < 0) {               /* if denominator neg */
                sign = -sign;           /* reverse the sign */
                NEG64(den);
        }

        /* now check whether numerator or denominator is larger */
        if (EQ64(num, den)) {
                num[0] = sign;          /* answer is +1 or -1 */

                if (sign < 0)           /* set high order word to sign bit */
                        num[1] = -1;
                else
                        num[1] = 0;

                den[1] = 0;             /* set den to 1 */
                den[0] = 1;

                return;
        }

        /* set up c and b so that one is num, the other den, and c > b */
        if (LT64(num, den)) {                   /* if num < den */
                c[0] = den[0];                  /* c = den */
                c[1] = den[1];
                b[0] = num[0];                  /* b = num */
                b[1] = num[1];
        } else {
                c[0] = num[0];                  /* c = num */
                c[1] = num[1];
                b[0] = den[0];                  /* b = den */
                b[1] = den[1];
        }

        /* use Euclidean Algorithm to find greatest common divisor of c & b */
        do {
                divmod64(c, b, junk, r);        /* r = c % b */
                c[0] = b[0];                    /* c = b */
                c[1] = b[1];
                b[0] = r[0];                    /* b = r */
                b[1] = r[1];
        } while (r[0] != 0 || r[1] != 0);       /* while r != 0 */

        /* now c is the greatest common divisor of num and den */

        /* divide out the greatest common divisor and put the sign in */
        divmod64(num, c, num, junk);            /* num /= c */
        if (sign < 0)                           /* if should be negative */
                NEG64(num);                     /* negate the numerator */
        divmod64(den, c, den, junk);            /* den /= c */

        return;
}