reduce-historical/r36/cslbase/arith01.c

/*  arith01.c                         Copyright (C) 1990 Codemist Ltd */

/*
 * Arithmetic functions.
 *        Addition of generic numbers
 *        and in particular a lot of bignum support.
 *
 * Version 1.4   November 1990.
 *
 */

/* Signature: 07d5049f 23-Jan-1999 */

#include <stdarg.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

#include "machine.h"
#include "tags.h"
#include "cslerror.h"
#include "externs.h"
#include "arith.h"
#ifdef TIMEOUT
#include "timeout.h"
#endif

/*
 * I start off with a collection of utility functions that create
 * Lisp structures to represent various sorts of numbers
 * and which extract values from same.
 * The typedefs that explain the layout of these structures are in "tags.h"
 */

Lisp_Object make_one_word_bignum(int32 n)
/*
 * n is an integer - create a bignum representation of it.  This is
 * done when n is outside the range 0xf8000000 to 0x07ffffff.
 */
{   Lisp_Object w = getvector(TAG_NUMBERS, TYPE_BIGNUM, 8);
    Lisp_Object nil;
    errexit();
    bignum_digits(w)[0] = n;
    return w;
}

Lisp_Object make_two_word_bignum(int32 a1, unsigned32 a0)
/*
 * This make a 2-word bignum from the 2-word value (a1,a0), where it
 * must have been arranged already that a1 and a0 are correctly
 * normalized to put in the two words as indicated.
 */
{
    Lisp_Object w = getvector(TAG_NUMBERS, TYPE_BIGNUM, 12);
    Lisp_Object nil;
    errexit();
    bignum_digits(w)[0] = a0;
    bignum_digits(w)[1] = a1;
    bignum_digits(w)[2] = 0;
    return w;
}

#ifdef COMMON
Lisp_Object make_sfloat(double d)
/*
 * Turn a regular floating point value into a Lisp "short float", which
 * is an immediate object obtained by using the bottom 4 bits of a 32-bit
 * word as tag, and the rest as just whatever would stand for a regular
 * single precision value.  In doing the conversion here I ignore
 * rounding etc - short floats are to save heap turn-over, but will
 * not give robust numeric results.
 */
{
    Float_union w;
    w.f = (float)d;
    return (w.i & ~(int32)0xf) + TAG_SFLOAT;
}
#endif

Lisp_Object make_boxfloat(double a, int32 type)
/*
 * Make a boxed float (single, double or long according to the type specifier)
 * if type==0 this makes a short float
 */
{
    Lisp_Object r, nil;
#ifndef COMMON
    CSL_IGNORE(type);
#endif
#ifdef COMMON
    switch (type)
    {
case 0:
        {   Float_union aa;
            aa.f = (float)a;
            return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
        }
case TYPE_SINGLE_FLOAT:
        r = getvector(TAG_BOXFLOAT, TYPE_SINGLE_FLOAT, sizeof(Single_Float));
        errexit();
        single_float_val(r) = (float)a;
        return r;
default: /* TYPE_DOUBLE_FLOAT I hope */
#endif
        r = getvector(TAG_BOXFLOAT, TYPE_DOUBLE_FLOAT, sizeof(Double_Float));
        errexit();
        double_float_val(r) = a;
        return r;
#ifdef COMMON
case TYPE_LONG_FLOAT:
        r = getvector(TAG_BOXFLOAT, TYPE_LONG_FLOAT, sizeof(Long_Float));
        errexit();
        long_float_val(r) = a;
        return r;
    }
#endif
}

static double bignum_to_float(Lisp_Object v, int32 h, int *xp)
/*
 * Convert a Lisp bignum to get a floating point value.  This uses at most the
 * top 3 digits of the bignum's representation since that is enough to achieve
 * full double precision accuracy.
 * This can not overflow, because it leaves an exponent-adjustment value
 * in *xp. You need ldexp(r, *xp) afterwards.
 */
{
    int32 n = (h>>2) - 2;  /* Last index into the data */
    int x = 31*(int)n;
    int32 msd = (int32)bignum_digits(v)[n];
/* NB signed conversion on next line */
    double r = (double)msd;
    switch (n)
    {
default:        /* for very big numbers combine in 3 digits */
        r = TWO_31*r + (double)bignum_digits(v)[--n];
        x -= 31;
        /* drop through */
case 1: r = TWO_31*r + (double)bignum_digits(v)[--n];
        x -= 31;
        /* drop through */
case 0: break;  /* do no more */
    }
    *xp = x;
    return r;
}

double float_of_number(Lisp_Object a)
/*
 * Return a (double precision) floating point value for the given Lisp
 * number, or 0.0 in case of trouble.  This is often called in circumstances
 * where I already know the type of its argument and so its type-dispatch
 * is unnecessary - in doing so I am trading off performance against
 * code repetition.
 */
{
    if (is_fixnum(a)) return (double)int_of_fixnum(a);
#ifdef COMMON
    else if (is_sfloat(a))
    {   Float_union w;
        w.i = a - TAG_SFLOAT;
        return (double)w.f;
    }
#endif
    else if (is_bfloat(a))
    {   int32 h = type_of_header(flthdr(a));
        switch (h)
        {
#ifdef COMMON
    case TYPE_SINGLE_FLOAT:
            return (double)single_float_val(a);
#endif
    case TYPE_DOUBLE_FLOAT:
            return double_float_val(a);
#ifdef COMMON
    case TYPE_LONG_FLOAT:
            return (double)long_float_val(a);
#endif
    default:
            return 0.0;
        }
    }
    else
    {   Header h = numhdr(a);
        int x1;
        double r1;
        switch (type_of_header(h))
        {
    case TYPE_BIGNUM:
            r1 = bignum_to_float(a, length_of_header(h), &x1);
            return ldexp(r1, x1);
#ifdef COMMON
    case TYPE_RATNUM:
            {   int x2;
                Lisp_Object na = numerator(a);
                a = denominator(a);
                if (is_fixnum(na)) r1 = float_of_number(na), x1 = 0;
                else r1 = bignum_to_float(na, 
                              length_of_header(numhdr(na)), &x1);
                if (is_fixnum(a)) r1 = r1 / float_of_number(a), x2 = 0;
                else r1 = r1 / bignum_to_float(a,
                                   length_of_header(numhdr(a)), &x2);
/* Floating point overflow can only arise in this ldexp() */
                return ldexp(r1, x1 - x2);
            }
#endif
    default:
/*
 * If the value was non-numeric or a complex number I hand back 0.0,
 * and since I am supposed to have checked the object type already
 * this OUGHT not to arrive - bit raising an exception seems over-keen.
 */
            return 0.0;
        }
    }
}

int32 thirty_two_bits(Lisp_Object a)
/*
 * return a 32 bit integer value for the Lisp integer (fixnum or bignum)
 * passed down - ignore any higher order bits and return 0 if the arg was
 * floating, rational etc or not a number at all.  Only really wanted where
 * links between C-specific code (that might really want 32-bit values)
 * and Lisp are being coded.
 */
{
    switch ((int)a & TAG_BITS)
    {
case TAG_FIXNUM:
        return int_of_fixnum(a);
case TAG_NUMBERS:
        if (is_bignum(a))
        {   int len = bignum_length(a);
/*
 * Note that I keep 31 bits per word and use a 2s complement representation.
 * thus if I have a one-word bignum I just want its contents but in all
 * other cases I need just one bit from the next word up.
 */
            if (len == 8) return bignum_digits(a)[0]; /* One word bignum */
            return bignum_digits(a)[0] | (bignum_digits(a)[1] << 31);
        }
        /* else drop through */
case TAG_BOXFLOAT:
default:
/*
 * return 0 for all non-fixnums
 */
        return 0;
    }
}

#ifdef COMMON
Lisp_Object make_complex(Lisp_Object r, Lisp_Object i)
{
    Lisp_Object v, nil = C_nil;
/*
 * Here r and i are expected to be either both rational (which in this
 * context includes the case of integer values) or both of the same
 * floating point type.  It is assumed that this has already been
 * arranged by here.
 */
    if (i == fixnum_of_int(0)) return r;
    stackcheck2(0, r, i);
    push2(r, i);
    v = getvector(TAG_NUMBERS, TYPE_COMPLEX_NUM, sizeof(Complex_Number));
/*
 * The vector r has uninitialized contents here - dodgy.  If the call
 * to getvector succeeded then I fill it in, otherwise I will not
 * refer to it again, and I think that unreferenced vectors containing junk
 * are OK.
 */
    pop2(i, r);
    errexit();
    real_part(v) = r;
    imag_part(v) = i;
    return v;
}

Lisp_Object make_ratio(Lisp_Object p, Lisp_Object q)
/*
 * By the time this is called (p/q) must be in its lowest terms, q>0
 */
{
    Lisp_Object v, nil = C_nil;
    if (q == fixnum_of_int(1)) return p;
    stackcheck2(0, p, q);
    push2(p, q);
    v = getvector(TAG_NUMBERS, TYPE_RATNUM, sizeof(Rational_Number));
    pop2(q, p);
    errexit();
    numerator(v) = p;
    denominator(v) = q;
    return v;
}

#endif

/*
 * The next bit of code seems pretty dreadful, but I think that is just
 * what generic arithmetic is all about.  The code for plus2 is written
 * as a dispatch function into over 30 separate possible type-specific
 * versions of the code.  In a very few simple (and performance-critical)
 * cases the code is written in-line in plus2 - in particular arithmetic
 * on fixnums is done that way.  Similarly for other cases.
 * I Use one-character suffices to remind me about types:
 *          i       fixnum
 *          b       bignum
 *          r       ratio
 *          s       short float
 *          f       boxed float (single/double/long)
 *          c       complex
 *
 * Throughout this code I am going to IGNORE floating point exceptions,
 * at least for a first attempt.  Decent detection of and recovery after
 * floating point overflow seems an extra unpleasant distraction!  Note
 * that C allows me to trap the SIGFPE exception, but returning from
 * the exception handler gives undefined behaviour - one is expected
 * to longjmp out, which means accepting the cost of using setjmp.
 *
 * It would perhaps be reasonable to write the dispatch code as a big
 * macro so that the versions for plus, times etc could all be kept
 * in step - I have not done that (a) because the macro would have been
 * bigger than I like macros to be (b) it would have involved token-
 * splicing (or VERY many parameters) to generate the names of the
 * separate type-specific procedures and (c) doing it by hand allows me
 * total flexibility about coding various cases in-line.
 */

#ifdef COMMON

static Lisp_Object plusis(Lisp_Object a, Lisp_Object b)
{
    Float_union bb;
    bb.i = b - TAG_SFLOAT;
    bb.f = (float)((float)int_of_fixnum(a) + bb.f);
    return (bb.i & ~(int32)0xf) + TAG_SFLOAT;
}
#endif

/*
 * Bignums are represented as vectors where the most significant 32-bit
 * digit is treated as signed, and the remaining ones are unsigned.
 */

static Lisp_Object plusib(Lisp_Object a, Lisp_Object b)
/*
 * Add a fixnum to a bignum, returning a result as a fixnum or bignum
 * depending on its size.  This seems much nastier than one would have
 * hoped.
 */
{
    int32 len = bignum_length(b), i, sign = int_of_fixnum(a), s;
    Lisp_Object c, nil;
    len = (len >> 2) - 1;
    if (len == 1)
    {   int32 t;
/*
 * Partly because it will be a common case and partly because it has
 * various special cases I have special purpose code to cope with
 * adding a fixnum to a one-word bignum.
 */
        s = (int32)bignum_digits(b)[0] + sign;
        t = s + s;
        if (top_bit_set(s ^ t))  /* needs to turn into two-word bignum */
        {   c = getvector(TAG_NUMBERS, TYPE_BIGNUM, 12);
            errexit();
            if (top_bit_set(s))
            {   bignum_digits(c)[0] = clear_top_bit(s);
                bignum_digits(c)[1] = (unsigned32)(-1);
            }
            else
            {   bignum_digits(c)[0] = s;
                bignum_digits(c)[1] = 0;
            }
            bignum_digits(c)[2] = 0;  /* set padding word to zero to be tidy */
            return c;
        }
        t = s & fix_mask;    /* Will it fit as a fixnum? */
        if (t == 0 || t == fix_mask) return fixnum_of_int(s);
        /* here the result is a one-word bignum */
        c = getvector(TAG_NUMBERS, TYPE_BIGNUM, 8);
        errexit();
        bignum_digits(c)[0] = s;
        return c;
    }
/*
 * Now, after all the silly cases have been handled, I have a calculation
 * which seems set to give a multi-word result.  The result here can at
 * least never shrink to a fixnum since subtracting a fixnum can at
 * most shrink the length of a number by one word.  I omit the stack-
 * check here in the hope that code here never nests enough for trouble.
 */
    push(b);
    c = getvector(TAG_NUMBERS, TYPE_BIGNUM, 4+(len<<2));
    pop(b);
    errexit();
    s = bignum_digits(b)[0] + clear_top_bit(sign);
    bignum_digits(c)[0] = clear_top_bit(s);
    if (sign >= 0) sign = 0; else sign = 0x7fffffff; /* extend the sign */

    len--;
    for (i=1; i<len; i++)
    {   s = bignum_digits(b)[i] + sign + top_bit(s);
        bignum_digits(c)[i] = clear_top_bit(s);
    }
    /* Now just the most significant digit remains to be processed  */
    if (sign != 0) sign = -1;
    {   s = bignum_digits(b)[i] + sign + top_bit(s);
        if (!signed_overflow(s))         /* did it overflow? */
        {
/*
 * Here the most significant digit did not produce an overflow, but maybe
 * what we actually had was some cancellation and the MSD is now zero
 * or -1, so that the number should shrink...
 */
            if ((s == 0 && (bignum_digits(c)[i-1] & 0x40000000) == 0) ||
                (s == -1 && (bignum_digits(c)[i-1] & 0x40000000) != 0))
            {   /* shrink the number */
                numhdr(c) -= pack_hdrlength(1L);
                if (s == -1) bignum_digits(c)[i-1] |= ~0x7fffffff;
/*
 * Now sometimes the shrinkage will leave a padding word, sometimes it
 * will really allow me to save space.
 */
                if ((i & 1) == 0)
                {   bignum_digits(c)[i] = 0;   /* leave the unused word tidy */
                    return c;
                }
/*
 * Having shrunk the number I am leaving a doubleword of unallocated space
 * in the heap.  Dump a header word into it to make it look like an
 * 8-byte bignum since that will allow the garbage collector to handle it.
 * It I left it containing arbitrary junk I could wreck myself.
 */
                bignum_digits(c)[i] = make_bighdr(2L);
                return c;
            }
            bignum_digits(c)[i] = s;  /* length unchanged */
            return c;
        }
/*
 * Here the result is one word longer than the input-bignum.
 * Once again SOMTIMES this will not involve allocating more store,
 * but just encroaching into the previously unused word that was padding
 * things out to a multiple of 8 bytes.
 */
        if ((i & 1) == 1)
        {   bignum_digits(c)[i++] = clear_top_bit(s);
            bignum_digits(c)[i] = top_bit_set(s) ? -1 : 0;
            numhdr(c) += pack_hdrlength(1L);
            return c;
        }
        push(c);
        b = getvector(TAG_NUMBERS, TYPE_BIGNUM, 12+(len<<2));
        pop(c);
        errexit();
        for (i=0; i<=len; i++)
            bignum_digits(b)[i] = bignum_digits(c)[i];
        bignum_digits(b)[i++] = clear_top_bit(s);
        bignum_digits(b)[i] = top_bit_set(s) ? -1 : 0;
        return b;
    }
}

#ifdef COMMON
static Lisp_Object plusir(Lisp_Object a, Lisp_Object b)
/*
 * fixnum and ratio, but also valid for bignum and ratio.
 * Note that if the inputs were in lowest terms there is no need for
 * and GCD calculations here.
 */
{
    Lisp_Object nil;
    push(b);
    a = times2(a, denominator(b));
    nil = C_nil;
    if (!exception_pending()) a = plus2(a, numerator(stack[0]));
    pop(b);
    errexit();
    return make_ratio(a, denominator(b));
}

static Lisp_Object plusic(Lisp_Object a, Lisp_Object b)
/*
 * real of any sort plus complex.
 */
{
    Lisp_Object nil;
    push(b);
    a = plus2(a, real_part(b));
    pop(b);
    errexit();
/*
 * make_complex() takes responsibility for mapping #C(n 0) onto n
 */
    return make_complex(a, imag_part(b));
}
#endif

static Lisp_Object plusif(Lisp_Object a, Lisp_Object b)
/*
 * Fixnum plus boxed-float.
 */
{
    double d = (double)int_of_fixnum(a) + float_of_number(b);
    return make_boxfloat(d, type_of_header(flthdr(b)));
}

#ifdef COMMON
#define plussi(a, b) plusis(b, a)

#define plussb(a, b) plusbs(b, a)

#define plussr(a, b) plusrs(b, a)

#define plussc(a, b) plusic(a, b)

#endif

static Lisp_Object plussf(Lisp_Object a, Lisp_Object b)
/*
 * This can be used for any rational value plus a boxed-float. plusif()
 * is separated just for (minor) efficiency reasons.
 */
{
    double d = float_of_number(a) + float_of_number(b);
    return make_boxfloat(d, type_of_header(flthdr(b)));
}

#define plusbi(a, b) plusib(b, a)

#ifdef COMMON

static Lisp_Object plusbs(Lisp_Object a, Lisp_Object b)
{
    double d = float_of_number(a) + float_of_number(b);
    return make_sfloat(d);
}

#endif

Lisp_Object lengthen_by_one_bit(Lisp_Object a, int32 msd)
/*
 * (a) is a bignum, and arithmetic on it has (just) caused overflow
 * in its top word - I just need to stick on another word. (msd) is the
 * current top word, and its sign will be used to decide on the value
 * that must be appended.
 */
{
    int32 len = bignum_length(a);
/*
 * Sometimes I need to allocate a new vector and copy data across into it
 */
    if ((len & 4) == 0)
    {   Lisp_Object b, nil;
        int32 i;
        push(a);
        b = getvector(TAG_NUMBERS, TYPE_BIGNUM, len+4);
        pop(a);
        errexit();
        len = (len >> 2) - 1;
        for (i=0; i<len; i++)
            bignum_digits(b)[i] = clear_top_bit(bignum_digits(a)[i]);
        bignum_digits(b)[i] = top_bit_set(msd) ? -1 : 0;
        bignum_digits(b)[i+1] = 0;
        return b;
    }
    else
/*
 * .. whereas sometimes I have a spare word already available.
 */
    {   numhdr(a) += pack_hdrlength(1L);
        len = (len >> 2) - 1;
        bignum_digits(a)[len-1] = clear_top_bit(bignum_digits(a)[len-1]);
        bignum_digits(a)[len] = top_bit_set(msd) ? -1 : 0;
        return a;
    }
}

static Lisp_Object plusbb(Lisp_Object a, Lisp_Object b)
/*
 * add two bignums.
 */
{
    int32 la = bignum_length(a),
          lb = bignum_length(b),
          i, s, carry;
    Lisp_Object c, nil;
    if (la < lb)    /* maybe swap order of args */
    {   Lisp_Object t = a;
        int32 t1;
        a = b; b = t;
        t1 = la; la = lb; lb = t1;
    }
/*
 * now (a) is AT LEAST as long as b.  I have special case code for
 * when both args are single-word bignums, since I expect that to be
 * an especially common case.
 */
    if (la == 8)    /* and hence b also has only 1 digit */
    {   int32 va = bignum_digits(a)[0],
              vb = bignum_digits(b)[0],
              vc = va + vb;
        if (signed_overflow(vc)) /* we have a 2-word bignum result */
        {   Lisp_Object w = getvector(TAG_NUMBERS, TYPE_BIGNUM, 12);
            errexit();
            bignum_digits(w)[0] = clear_top_bit(vc);
            bignum_digits(w)[1] = top_bit_set(vc) ? -1 : 0;
            bignum_digits(w)[2] = 0;
            return w;
        }
/*
 * here the result fits into one word - maybe it will squash down into
 * a fixnum?
 */
        else
        {   vb = vc & fix_mask;
            if (vb == 0 || vb == fix_mask) return fixnum_of_int(vc);
            else return make_one_word_bignum(vc);
        }
    }
    push2(a, b);
    c = getvector(TAG_NUMBERS, TYPE_BIGNUM, la);
    pop2(b, a);
    errexit();
    la = (la >> 2) - 2;
    lb = (lb >> 2) - 2;
    carry = 0;
/*
 * Add all but the top digit of b
 */
    for (i=0; i<lb; i++)
    {   carry = bignum_digits(a)[i] + bignum_digits(b)[i] + top_bit(carry);
        bignum_digits(c)[i] = clear_top_bit(carry);
    }
    if (la == lb) s = bignum_digits(b)[i];
    else
/*
 * If a is strictly longer than b I sign extend b here and add in as many
 * copies of 0 or -1 as needbe to get up to the length of a.
 */
    {   s = bignum_digits(b)[i];
        carry =  bignum_digits(a)[i] + clear_top_bit(s) + top_bit(carry);
        bignum_digits(c)[i] = clear_top_bit(carry);
        if (s < 0) s = -1; else s = 0;
        for (i++; i<la; i++)
        {   carry = bignum_digits(a)[i] + clear_top_bit(s) + top_bit(carry);
            bignum_digits(c)[i] = clear_top_bit(carry);
        }
    }
/*
 * the most significant digit is added using signed arithmetic so that I
 * can tell if it overflowed.
 */
    carry = bignum_digits(a)[i] + s + top_bit(carry);
    if (!signed_overflow(carry))
    {
/*
 * Here the number has not expanded - but it may be shrinking, and it can
 * shrink by any number of words, all the way down to a fixnum maybe.  Note
 * that I started with at least a 2-word bignum here.
 */
        int32 msd;
        bignum_digits(c)[i] = carry;
        if (carry == 0)
        {   int32 j = i-1;
            while ((msd = bignum_digits(c)[j]) == 0 && j > 0) j--;
/*
 * ... but I may need a zero word on the front if the next word down
 * has its top bit set... (top of 31 bits, that is)
 */
            if ((msd & 0x40000000) != 0)
            {   j++;
                if (i == j) return c;
            }
            if (j == 0)
            {   int32 s = bignum_digits(c)[0];
                if ((s & fix_mask) == 0) return fixnum_of_int(s);
            }
/*
 * If I am shrinking by one word and had an even length to start with
 * I do not have to mess about so much.
 */
            if (j == i-1 && (i & 1) == 0)
            {   numhdr(c) -= pack_hdrlength(1L);
                return c;
            }
            numhdr(c) -= pack_hdrlength(i - j);
            i = (i+1) | 1;
            j = (j+1) | 1;     /* Round up to odd index */
/*
 * I forge a header word to allow the garbage collector to skip over
 * (and in due course reclaim) the space that turned out not to be needed.
 */
            bignum_digits(c)[j] = make_bighdr(i - j);
            return c;
        }
/*
 * Now do all the same sorts of things but this time for negative numbers.
 */
        else if (carry == -1)
        {   int32 j = i-1;
            msd = carry;    /* in case j = 0 */
            while ((msd = bignum_digits(c)[j]) == 0x7fffffff && j > 0) j--;
            if ((msd & 0x40000000) == 0)
            {   j++;
                if (i == j) return c;
            }
            if (j == 0)
            {   int32 s = bignum_digits(c)[0] | ~0x7fffffff;
                if ((s & fix_mask) == fix_mask) return fixnum_of_int(s);
            }
            if (j == i-1 && (i & 1) == 0)
            {   bignum_digits(c)[i] = 0;
                bignum_digits(c)[i-1] |= ~0x7fffffff;
                numhdr(c) -= pack_hdrlength(1);
                return c;
            }
            numhdr(c) -= pack_hdrlength(i - j);
            bignum_digits(c)[j+1] = 0;
            bignum_digits(c)[j] |= ~0x7fffffff;
            i = (i+1) | 1;
            j = (j+1) | 1;     /* Round up to odd index */
            bignum_digits(c)[j] = make_bighdr(i - j);
            return c;
        }
        return c;
    }
    else
    {   bignum_digits(c)[i] = carry;
        return lengthen_by_one_bit(c, carry);
    }
}

#ifdef COMMON
#define plusbr(a, b) plusir(a, b)

#define plusbc(a, b) plusic(a, b)
#endif

#define plusbf(a, b) plussf(a, b)

#ifdef COMMON
#define plusri(a, b) plusir(b, a)

#define plusrs(a, b) plusbs(a, b)

#define plusrb(a, b) plusri(a, b)

static Lisp_Object plusrr(Lisp_Object a, Lisp_Object b)
/*
 * Adding two ratios involves some effort to keep the result in
 * lowest terms.
 */
{
    Lisp_Object nil = C_nil;
    Lisp_Object na = numerator(a), nb = numerator(b);
    Lisp_Object da = denominator(a), db = denominator(b);
    Lisp_Object w = nil;
    push5(na, nb, da, db, nil);
#define g   stack[0]
#define db  stack[-1]
#define da  stack[-2]
#define nb  stack[-3]
#define na  stack[-4]
    g = gcd(da, db);
    nil = C_nil;
    if (exception_pending()) goto fail;
/*
 * all the calls to quot2() in this procedure are expected - nay required -
 * to give exact integer quotients.
 */
    db = quot2(db, g);
    nil = C_nil;
    if (exception_pending()) goto fail;
    g = quot2(da, g);
    nil = C_nil;
    if (exception_pending()) goto fail;
    na = times2(na, db);
    nil = C_nil;
    if (exception_pending()) goto fail;
    nb = times2(nb, g);
    nil = C_nil;
    if (exception_pending()) goto fail;
    na = plus2(na, nb);
    nil = C_nil;
    if (exception_pending()) goto fail;
    da = times2(da, db);
    nil = C_nil;
    if (exception_pending()) goto fail;
    g = gcd(na, da);
    nil = C_nil;
    if (exception_pending()) goto fail;
    na = quot2(na, g);
    nil = C_nil;
    if (exception_pending()) goto fail;
    da = quot2(da, g);
    nil = C_nil;
    if (exception_pending()) goto fail;
    w = make_ratio(na, da);
/*
 * All the goto statements and the label seem a fair way of expressing
 * the common action that has to be taken if an error or interrupt is
 * detected during any of the intermediate steps here.  Anyone who
 * objects can change it if they really want...
 */
fail:
    popv(5);
    return w;
#undef na
#undef nb
#undef da
#undef db
#undef g
}

#define plusrc(a, b) plusic(a, b)

#define plusrf(a, b) plussf(a, b)

#define plusci(a, b) plusic(b, a)

#define pluscs(a, b) plussc(b, a)

#define pluscb(a, b) plusbc(b, a)

#define pluscr(a, b) plusrc(b, a)

static Lisp_Object pluscc(Lisp_Object a, Lisp_Object b)
/*
 * Add complex values.
 */
{
    Lisp_Object c, nil;
    push2(a, b);
    c = plus2(imag_part(a), imag_part(b));
    pop2(b, a);
    errexit();
    a = plus2(real_part(a), real_part(b));
    errexit();
    return make_complex(a, c);
}

#define pluscf(a, b) plusfc(b, a)
#endif

#define plusfi(a, b) plusif(b, a)

#ifdef COMMON
#define plusfs(a, b) plussf(b, a)
#endif

#define plusfb(a, b) plusbf(b, a)

#ifdef COMMON
#define plusfr(a, b) plusrf(b, a)

#define plusfc(a, b) plusic(a, b)
#endif

static Lisp_Object plusff(Lisp_Object a, Lisp_Object b)
/*
 * Add two boxed floats - the type of the result must match the
 * longer of the types of the arguments, hence the extra
 * messing about.
 */
{
#ifdef COMMON
    int32 ha = type_of_header(flthdr(a)), hb = type_of_header(flthdr(b));
#endif
    double d;
/*
 * This is written as a declaration followed by a separate assignment to
 * d because I hit a compiler bug on a VAX once otherwise.
 */
    d = float_of_number(a) + float_of_number(b);
#ifdef COMMON
    if (ha == TYPE_LONG_FLOAT || hb == TYPE_LONG_FLOAT)
        ha = TYPE_LONG_FLOAT;
    else if (ha == TYPE_DOUBLE_FLOAT || hb == TYPE_DOUBLE_FLOAT)
        ha = TYPE_DOUBLE_FLOAT;
    else ha = TYPE_SINGLE_FLOAT;
    return make_boxfloat(d, ha);
#else
    return make_boxfloat(d, TYPE_DOUBLE_FLOAT);
#endif
}

/*
 * and now for the dispatch code...
 */

/*
 * The following verifies that a number is properly formatted - a
 * fixnum if small enough or a decently normalised bignum.  For use when
 * there is suspicion of a bug wrt such matters. Call is
 *   validate_number("msg", numberToCheck, nX, nY)
 * where nX and nY must be numbers and are shown in any
 * diagnostic.
 */

void validate_number(char *s, Lisp_Object a, Lisp_Object b, Lisp_Object c)
{
    int32 la, w, msd;
    if (!is_numbers(a)) return;
    la = (length_of_header(numhdr(a))>>2)-2;
    if (la < 0)
        {   trace_printf("%s: number with no digits (%.8x)\n", s, numhdr(a));
            if (is_number(b)) prin_to_trace(b), trace_printf("\n");
            if (is_number(c)) prin_to_trace(c), trace_printf("\n");
            my_exit(EXIT_FAILURE);
        }
    if (la == 0)
    {   msd = bignum_digits(a)[0];
        w = msd & fix_mask;
        if (w == 0 || w == fix_mask)
        {   trace_printf("%s: %.8x should be fixnum\n", s, msd);
            if (is_number(b)) prin_to_trace(b), trace_printf("\n");
            if (is_number(c)) prin_to_trace(c), trace_printf("\n");
            my_exit(EXIT_FAILURE);
        }
        if (signed_overflow(msd))
        {   trace_printf("%s: %.8x should be two-word\n", s, msd);
            if (is_number(b)) prin_to_trace(b), trace_printf("\n");
            if (is_number(c)) prin_to_trace(c), trace_printf("\n");
            my_exit(EXIT_FAILURE);
        }
        return;
    }
    msd = bignum_digits(a)[la];
    if (signed_overflow(msd))
    {   trace_printf("%s: %.8x should be longer\n", s, msd);
        if (is_number(b)) prin_to_trace(b), trace_printf("\n");
        if (is_number(c)) prin_to_trace(c), trace_printf("\n");
        my_exit(EXIT_FAILURE);
    }
    if (msd == 0 && ((msd = bignum_digits(a)[la-1]) & 0x40000000) == 0)
    {   trace_printf("%s: 0: %.8x should be shorter\n", s, msd);
        if (is_number(b)) prin_to_trace(b), trace_printf("\n");
        if (is_number(c)) prin_to_trace(c), trace_printf("\n");
        my_exit(EXIT_FAILURE);
    }
    if (msd == -1 && ((msd = bignum_digits(a)[la-1]) & 0x40000000) != 0)
    {   trace_printf("%s: -1: %.8x should be shorter\n", s, msd);
        if (is_number(b)) prin_to_trace(b), trace_printf("\n");
        if (is_number(c)) prin_to_trace(c), trace_printf("\n");
        my_exit(EXIT_FAILURE);
    }
}


Lisp_Object plus2(Lisp_Object a, Lisp_Object b)
/*
 * I probably want to change the specification of plus2 so that the fixnum +
 * fixnum case is always expected to be done before the main body of the code
 * is entered.  Well maybe even if I do that it then costs very little to
 * include the fixnum code here as well, so I will not delete it.
 */
{
    switch ((int)a & TAG_BITS)
    {
case TAG_FIXNUM:
        switch ((int)b & TAG_BITS)
        {
    case TAG_FIXNUM:
/*
 * This is where fixnum + fixnum arithmetic happens - the case I most want to
 * make efficient.
 */
            {   int32 r = int_of_fixnum(a) + int_of_fixnum(b);
                int32 t = r & fix_mask;
                if (t == 0 || t == fix_mask) return fixnum_of_int(r);
                else return make_one_word_bignum(r);
            }
#ifdef COMMON
    case TAG_SFLOAT:
            return plusis(a, b);
#endif
    case TAG_NUMBERS:
            {   int32 hb = type_of_header(numhdr(b));
                switch (hb)
                {
        case TYPE_BIGNUM:
                return plusib(a, b);
#ifdef COMMON
        case TYPE_RATNUM:
                return plusir(a, b);
        case TYPE_COMPLEX_NUM:
                return plusic(a, b);
#endif
        default:
                return aerror1("bad arg for plus", b);
                }
            }
    case TAG_BOXFLOAT:
            return plusif(a, b);
    default:
            return aerror1("bad arg for plus",  b);
        }
#ifdef COMMON
case TAG_SFLOAT:
        switch (b & TAG_BITS)
        {
    case TAG_FIXNUM:
            return plussi(a, b);
    case TAG_SFLOAT:
            {   Float_union aa, bb;
                aa.i = a - TAG_SFLOAT;
                bb.i = b - TAG_SFLOAT;
                aa.f = (float)(aa.f + bb.f);
                return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
            }
    case TAG_NUMBERS:
            {   int32 hb = type_of_header(numhdr(b));
                switch (hb)
                {
        case TYPE_BIGNUM:
                return plussb(a, b);
        case TYPE_RATNUM:
                return plussr(a, b);
        case TYPE_COMPLEX_NUM:
                return plussc(a, b);
        default:
                return aerror1("bad arg for plus",  b);
                }
            }
    case TAG_BOXFLOAT:
            return plussf(a, b);
    default:
            return aerror1("bad arg for plus",  b);
        }
#endif
case TAG_NUMBERS:
        {   int32 ha = type_of_header(numhdr(a));
            switch (ha)
            {
    case TYPE_BIGNUM:
                switch ((int)b & TAG_BITS)
                {
            case TAG_FIXNUM:
                    return plusbi(a, b);
#ifdef COMMON
            case TAG_SFLOAT:
                    return plusbs(a, b);
#endif
            case TAG_NUMBERS:
                    {   int32 hb = type_of_header(numhdr(b));
                        switch (hb)
                        {
                case TYPE_BIGNUM:
                        return plusbb(a, b);
#ifdef COMMON
                case TYPE_RATNUM:
                        return plusbr(a, b);
                case TYPE_COMPLEX_NUM:
                        return plusbc(a, b);
#endif
                default:
                        return aerror1("bad arg for plus",  b);
                        }
                    }
            case TAG_BOXFLOAT:
                    return plusbf(a, b);
            default:
                    return aerror1("bad arg for plus",  b);
                }
#ifdef COMMON
    case TYPE_RATNUM:
                switch (b & TAG_BITS)
                {
            case TAG_FIXNUM:
                    return plusri(a, b);
            case TAG_SFLOAT:
                    return plusrs(a, b);
            case TAG_NUMBERS:
                    {   int32 hb = type_of_header(numhdr(b));
                        switch (hb)
                        {
                case TYPE_BIGNUM:
                        return plusrb(a, b);
                case TYPE_RATNUM:
                        return plusrr(a, b);
                case TYPE_COMPLEX_NUM:
                        return plusrc(a, b);
                default:
                        return aerror1("bad arg for plus",  b);
                        }
                    }
            case TAG_BOXFLOAT:
                    return plusrf(a, b);
            default:
                    return aerror1("bad arg for plus",  b);
                }
    case TYPE_COMPLEX_NUM:
                switch (b & TAG_BITS)
                {
            case TAG_FIXNUM:
                    return plusci(a, b);
            case TAG_SFLOAT:
                    return pluscs(a, b);
            case TAG_NUMBERS:
                    {   int32 hb = type_of_header(numhdr(b));
                        switch (hb)
                        {
                case TYPE_BIGNUM:
                        return pluscb(a, b);
                case TYPE_RATNUM:
                        return pluscr(a, b);
                case TYPE_COMPLEX_NUM:
                        return pluscc(a, b);
                default:
                        return aerror1("bad arg for plus",  b);
                        }
                    }
            case TAG_BOXFLOAT:
                    return pluscf(a, b);
            default:
                    return aerror1("bad arg for plus",  b);
                }
#endif
    default:    return aerror1("bad arg for plus",  a);
            }
        }
case TAG_BOXFLOAT:
        switch ((int)b & TAG_BITS)
        {
    case TAG_FIXNUM:
            return plusfi(a, b);
#ifdef COMMON
    case TAG_SFLOAT:
            return plusfs(a, b);
#endif
    case TAG_NUMBERS:
            {   int32 hb = type_of_header(numhdr(b));
                switch (hb)
                {
        case TYPE_BIGNUM:
                return plusfb(a, b);
#ifdef COMMON
        case TYPE_RATNUM:
                return plusfr(a, b);
        case TYPE_COMPLEX_NUM:
                return plusfc(a, b);
#endif
        default:
                return aerror1("bad arg for plus",  b);
                }
            }
    case TAG_BOXFLOAT:
            return plusff(a, b);
    default:
            return aerror1("bad arg for plus",  b);
        }
default:
        return aerror1("bad arg for plus",  a);
    }
}

Lisp_Object difference2(Lisp_Object a, Lisp_Object b)
{
    Lisp_Object nil;
    switch ((int)b & TAG_BITS)
    {
case TAG_FIXNUM:
        if (is_fixnum(a))
        {
            int32 r = int_of_fixnum(a) - int_of_fixnum(b);
            int32 t = r & fix_mask;
            if (t == 0 || t == fix_mask) return fixnum_of_int(r);
            else return make_one_word_bignum(r);
        }
        else if (b != ~0x7ffffffe) return plus2(a, 2*TAG_FIXNUM-b);
        else
        {   push(a);
            b = make_one_word_bignum(-int_of_fixnum(b));
            break;
        }
case TAG_NUMBERS:
        push(a);
        if (type_of_header(numhdr(b)) == TYPE_BIGNUM) b = negateb(b);
        else b = negate(b);
        break;
case TAG_BOXFLOAT:
default:
        push(a);
        b = negate(b);
        break;
    }
    pop(a);
    errexit();
    return plus2(a, b);
}

Lisp_Object add1(Lisp_Object p)
/*
 * Increment a number.  Short cut when the number is a fixnum, otherwise
 * just hand over to the general addition code.
 */
{
    if (is_fixnum(p))
    {   unsigned32 r = (int32)p + 0x10;
                                  /* fixnums have data shifted left 4 bits */
        if (r == ~0x7ffffffe)     /* The ONLY possible overflow case here  */
            return make_one_word_bignum(1 + int_of_fixnum(p));
        else return (Lisp_Object)r;
    }
    else return plus2(p, fixnum_of_int(1));
}

Lisp_Object sub1(Lisp_Object p)
/*
 * Decrement a number.  Short cut when the number is a fixnum, otherwise
 * just hand over to the general addition code.
 */
{
    if (is_fixnum(p))
    {   if (p == ~0x7ffffffe)     /* The ONLY possible overflow case here  */
            return make_one_word_bignum(int_of_fixnum(p) - 1);
        else return (Lisp_Object)(p - 0x10);
    }
    else return plus2(p, fixnum_of_int(-1));
}

/* end of arith01.c */