Usenet 1994 October/usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso/unix/volume27/calc-2.9.0/part12

Text File | 1993-12-07 | 61KB | 2,783 lines

Newsgroups: comp.sources.unix
From: [e-mail] (David I. Bell)
Subject: v27i139: calc-2.9.0 - arbitrary precision C-like programmable calculator, Part12/19
References: <[e-mail]>
Sender: [e-mail]
Approved: [e-mail]

Submitted-By: [e-mail] (David I. Bell)
Posting-Number: Volume 27, Issue 139
Archive-Name: calc-2.9.0/part12

#!/bin/sh
# this is part 12 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/zfunc.c continued
#
CurArch=12
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
     exit 1; fi
( read Scheck
  if test "$Scheck" != $CurArch
  then echo "Please unpack part $Scheck next!"
       exit 1;
  else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file calc2.9.0/zfunc.c"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/zfunc.c
X     * is unimportant.
X     */
X    highbit = (highbit + k - 1) / k;
X    try.len = (highbit / BASEB) + 1;
X    try.v = alloc(try.len);
X    try.v[try.len-1] = ((HALF)1 << (highbit % BASEB));
X    try.sign = 0;
X    old.v = alloc(try.len);
X    old.len = 1;
X    *old.v = 0;
X    old.sign = 0;
X    /*
X     * Main divide and average loop
X     */
X    for (;;) {
X        zpowi(try, k1, &temp);
X        zquo(z1, temp, &quo);
X        zfree(temp);
X        i = zrel(try, quo);
X        if (i <= 0) {
X            /*
X             * Current try is less than or equal to the root since it is
X             * less than the quotient. If the quotient is equal to the try,
X             * we are all done.  Also, if the try is equal to the old value,
X             * we are done since no improvement occurred.
X             * If not, save the improved value and loop some more.
X             */
X            if ((i == 0) || (zcmp(old, try) == 0)) {
X                zfree(quo);
X                zfree(old);
X                try.sign = (HALF)sign;
X                zquicktrim(try);
X                *dest = try;
X                return;
X            }
X            old.len = try.len;
X            zcopyval(try, old);
X        }
X        /* average current try and quotent for the new try */
X        zmul(try, k1, &temp);
X        zfree(try);
X        zadd(quo, temp, &temp2);
X        zfree(temp);
X        zfree(quo);
X        zquo(temp2, z2, &try);
X        zfree(temp2);
X    }
X}
X
X
X/*
X * Test to see if a number is an exact square or not.
X */
XBOOL
Xzissquare(z)
X    ZVALUE z;        /* number to be tested */
X{
X    long n, i;
X    ZVALUE tmp;
X
X    if (z.sign)        /* negative */
X        return FALSE;
X    while ((z.len > 1) && (*z.v == 0)) {    /* ignore trailing zero words */
X        z.len--;
X        z.v++;
X    }
X    if (zisleone(z))    /* zero or one */
X        return TRUE;
X    n = *z.v & 0xf;        /* check mod 16 values */
X    if ((n != 0) && (n != 1) && (n != 4) && (n != 9))
X        return FALSE;
X    n = *z.v & 0xff;    /* check mod 256 values */
X    i = 0x80;
X    while (((i * i) & 0xff) != n)
X        if (--i <= 0)
X            return FALSE;
X    n = zsqrt(z, &tmp);    /* must do full square root test now */
X    zfree(tmp);
X    return n;
X}
X
X
X/*
X * Return a trivial hash value for an integer.
X */
XHASH
Xzhash(z)
X    ZVALUE z;
X{
X    HASH hash;
X    int i;
X
X    hash = z.len * 1000003;
X    if (z.sign)
X        hash += 10000019;
X    for (i = z.len - 1; i >= 0; i--)
X        hash = hash * 79372817 + z.v[i] + 10000079;
X    return hash;
X}
X
X/* END CODE */
SHAR_EOF
echo "File calc2.9.0/zfunc.c is complete"
chmod 0644 calc2.9.0/zfunc.c || echo "restore of calc2.9.0/zfunc.c fails"
set `wc -c calc2.9.0/zfunc.c`;Sum=$1
if test "$Sum" != "34638"
then echo original size 34638, current size $Sum;fi
echo "x - extracting calc2.9.0/zio.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/zio.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Scanf and printf routines for arbitrary precision integers.
X */
X
X#include "stdarg.h"
X#include "zmath.h"
X
X
X#define    OUTBUFSIZE    200        /* realloc size for output buffers */
X
X#define    PUTCHAR(ch)        math_chr(ch)
X#define    PUTSTR(str)        math_str(str)
X#define    PRINTF1(fmt, a1)    math_fmt(fmt, a1)
X#define    PRINTF2(fmt, a1, a2)    math_fmt(fmt, a1, a2)
X
X
Xlong    _outdigits_ = 20;        /* default digits for output */
Xint    _outmode_ = MODE_INITIAL;    /* default output mode */
X
X
X/*
X * Output state that has been saved when diversions are done.
X */
Xtypedef struct iostate IOSTATE;
Xstruct iostate {
X    IOSTATE *oldiostates;        /* previous saved state */
X    long outdigits;            /* digits for output */
X    int outmode;            /* output mode */
X    FILE *outfp;            /* file unit for output (if any) */
X    char *outbuf;            /* output string buffer (if any) */
X    long outbufsize;        /* current size of string buffer */
X    long outbufused;        /* space used in string buffer */
X    BOOL outputisstring;        /* TRUE if output is to string buffer */
X};
X
X
Xstatic IOSTATE    *oldiostates = NULL;    /* list of saved output states */
Xstatic FILE    *outfp = stdout;    /* file unit for output */
Xstatic char    *outbuf = NULL;        /* current diverted buffer */
Xstatic BOOL    outputisstring = FALSE;
Xstatic long    outbufsize;
Xstatic long    outbufused;
X
X
X/*
X * Routine to output a character either to a FILE
X * handle or into a string.
X */
Xvoid
Xmath_chr(ch)
X{
X    char    *cp;
X
X    if (!outputisstring) {
X        fputc(ch, outfp);
X        return;
X    }
X    if (outbufused >= outbufsize) {
X        cp = (char *)realloc(outbuf, outbufsize + OUTBUFSIZE + 1);
X        if (cp == NULL)
X            math_error("Cannot realloc output string");
X        outbuf = cp;
X        outbufsize += OUTBUFSIZE;
X    }
X    outbuf[outbufused++] = (char)ch;
X}
X
X
X/*
X * Routine to output a null-terminated string either
X * to a FILE handle or into a string.
X */
Xvoid
Xmath_str(str)
X    char    *str;
X{
X    char    *cp;
X    int    len;
X
X    if (!outputisstring) {
X        fputs(str, outfp);
X        return;
X    }
X    len = strlen(str);
X    if ((outbufused + len) > outbufsize) {
X        cp = (char *)realloc(outbuf, outbufsize + len + OUTBUFSIZE + 1);
X        if (cp == NULL)
X            math_error("Cannot realloc output string");
X        outbuf = cp;
X        outbufsize += (len + OUTBUFSIZE);
X    }
X    memcpy(&outbuf[outbufused], str, len);
X    outbufused += len;
X}
X
X
X/*
X * Output a null-terminated string either to a FILE handle or into a string,
X * padded with spaces as needed so as to fit within the specified width.
X * If width is positive, the spaces are added at the front of the string.
X * If width is negative, the spaces are added at the end of the string.
X * The complete string is always output, even if this overflows the width.
X * No characters within the string are handled specially.
X */
Xvoid
Xmath_fill(str, width)
X    char *str;
X    long width;
X{
X    if (width > 0) {
X        width -= strlen(str);
X        while (width-- > 0)
X            PUTCHAR(' ');
X        PUTSTR(str);
X    } else {
X        width += strlen(str);
X        PUTSTR(str);
X        while (width++ < 0)
X            PUTCHAR(' ');
X    }
X}
X
X
X/*
X * Routine to output a printf-style formatted string either
X * to a FILE handle or into a string.
X */
X#ifdef VARARGS
X# define VA_ALIST fmt, va_alist
X# define VA_DCL char *fmt; va_dcl
X#else
X# ifdef __STDC__
X#  define VA_ALIST char *fmt, ...
X#  define VA_DCL
X# else
X#  define VA_ALIST fmt
X#  define VA_DCL char *fmt;
X# endif
X#endif
X/*VARARGS*/
Xvoid
Xmath_fmt(VA_ALIST)
X    VA_DCL
X{
X    va_list ap;
X    char buf[200];
X
X#ifdef VARARGS
X    va_start(ap);
X#else
X    va_start(ap, fmt);
X#endif
X    vsprintf(buf, fmt, ap);
X    va_end(ap);
X    math_str(buf);
X}
X
X
X/*
X * Flush the current output stream.
X */
Xvoid
Xmath_flush()
X{
X    if (!outputisstring)
X        fflush(outfp);
X}
X
X
X/*
X * Divert further output so that it is saved into a string that will be
X * returned later when the diversion is completed.  The current state of
X * output is remembered for later restoration.  Diversions can be nested.
X * Output diversion is only intended for saving output to "stdout".
X */
Xvoid
Xmath_divertio()
X{
X    register IOSTATE *sp;
X
X    sp = (IOSTATE *) malloc(sizeof(IOSTATE));
X    if (sp == NULL)
X        math_error("No memory for diverting output");
X    sp->oldiostates = oldiostates;
X    sp->outdigits = _outdigits_;
X    sp->outmode = _outmode_;
X    sp->outfp = outfp;
X    sp->outbuf = outbuf;
X    sp->outbufsize = outbufsize;
X    sp->outbufused = outbufused;
X    sp->outputisstring = outputisstring;
X
X    outbufused = 0;
X    outbufsize = 0;
X    outbuf = (char *) malloc(OUTBUFSIZE + 1);
X    if (outbuf == NULL)
X        math_error("Cannot allocate divert string");
X    outbufsize = OUTBUFSIZE;
X    outputisstring = TRUE;
X    oldiostates = sp;
X}
X
X
X/*
X * Undivert output and return the saved output as a string.  This also
X * restores the output state to what it was before the diversion began.
X * The string needs freeing by the caller when it is no longer needed.
X */
Xchar *
Xmath_getdivertedio()
X{
X    register IOSTATE *sp;
X    char *cp;
X
X    sp = oldiostates;
X    if (sp == NULL)
X        math_error("No diverted state to restore");
X    cp = outbuf;
X    cp[outbufused] = '\0';
X    oldiostates = sp->oldiostates;
X    _outdigits_ = sp->outdigits;
X    _outmode_ = sp->outmode;
X    outfp = sp->outfp;
X    outbuf = sp->outbuf;
X    outbufsize = sp->outbufsize;
X    outbufused = sp->outbufused;
X    outbuf = sp->outbuf;
X    outputisstring = sp->outputisstring;
X    return cp;
X}
X
X
X/*
X * Clear all diversions and set output back to the original destination.
X * This is called when resetting the global state of the program.
X */
Xvoid
Xmath_cleardiversions()
X{
X    while (oldiostates)
X        free(math_getdivertedio());
X}
X
X
X/*
X * Set the output routines to output to the specified FILE stream.
X * This interacts with output diversion in the following manner.
X *    STDOUT    diversion    action
X *    ----    ---------    ------
X *    yes    yes        set output to diversion string again.
X *    yes    no        set output to stdout.
X *    no    yes        set output to specified file.
X *    no    no        set output to specified file.
X */
Xvoid
Xmath_setfp(newfp)
X    FILE *newfp;
X{
X    outfp = newfp;
X    outputisstring = (oldiostates && (newfp == stdout));
X}
X
X
X/*
X * Set the output mode for numeric output.
X * This also returns the previous mode.
X */
Xint
Xmath_setmode(newmode)
X{
X    int oldmode;
X
X    if ((newmode <= MODE_DEFAULT) || (newmode > MODE_MAX))
X        math_error("Setting illegal output mode");
X    oldmode = _outmode_;
X    _outmode_ = newmode;
X    return oldmode;
X}
X
X
X/*
X * Set the number of digits for float or exponential output.
X * This also returns the previous number of digits.
X */
Xlong
Xmath_setdigits(newdigits)
X    long newdigits;
X{
X    long olddigits;
X
X    if (newdigits < 0)
X        math_error("Setting illegal number of digits");
X    olddigits = _outdigits_;
X    _outdigits_ = newdigits;
X    return olddigits;
X}
X
X
X/*
X * Print an integer value as a hex number.
X * Width is the number of columns to print the number in, including the
X * sign if required.  If zero, no extra output is done.  If positive,
X * leading spaces are typed if necessary. If negative, trailing spaces are
X * typed if necessary.  The special characters 0x appear to indicate the
X * number is hex.
X */
X/*ARGSUSED*/
Xvoid
Xzprintx(z, width)
X    ZVALUE z;
X    long width;
X{
X    register HALF *hp;    /* current word to print */
X    int len;        /* number of halfwords to type */
X    char *str;
X
X    if (width) {
X        math_divertio();
X        zprintx(z, 0L);
X        str = math_getdivertedio();
X        math_fill(str, width);
X        free(str);
X        return;
X    }
X    len = z.len - 1;
X    if (zisneg(z))
X        PUTCHAR('-');
X    if ((len == 0) && (*z.v <= (FULL) 9)) {
X        len = '0' + *z.v;
X        PUTCHAR(len);
X        return;
X    }
X    hp = z.v + len;
X    PRINTF1("0x%x", (FULL) *hp--);
X    while (--len >= 0)
X        PRINTF1("%04x", (FULL) *hp--);
X}
X
X
X/*
X * Print an integer value as a binary number.
X * The special characters 0b appear to indicate the number is binary.
X */
X/*ARGSUSED*/
Xvoid
Xzprintb(z, width)
X    ZVALUE z;
X    long width;
X{
X    register HALF *hp;    /* current word to print */
X    int len;        /* number of halfwords to type */
X    HALF val;        /* current value */
X    HALF mask;        /* current mask */
X    int didprint;        /* nonzero if printed some digits */
X    int ch;            /* current char */
X    char *str;
X
X    if (width) {
X        math_divertio();
X        zprintb(z, 0L);
X        str = math_getdivertedio();
X        math_fill(str, width);
X        free(str);
X        return;
X    }
X    len = z.len - 1;
X    if (zisneg(z))
X        PUTCHAR('-');
X    if ((len == 0) && (*z.v <= (FULL) 1)) {
X        len = '0' + *z.v;
X        PUTCHAR(len);
X        return;
X    }
X    hp = z.v + len;
X    didprint = 0;
X    PUTSTR("0b");
X    while (len-- >= 0) {
X        val = *hp--;
X        mask = (1 << (BASEB - 1));
X        while (mask) {
X            ch = '0' + ((mask & val) != 0);
X            if (didprint || (ch != '0')) {
X                PUTCHAR(ch);
X                didprint = 1;
X            }
X            mask >>= 1;
X        }
X    }
X}
X
X
X/*
X * Print an integer value as an octal number.
X * The number begins with a leading 0 to indicate that it is octal.
X */
X/*ARGSUSED*/
Xvoid
Xzprinto(z, width)
X    ZVALUE z;
X    long width;
X{
X    register HALF *hp;    /* current word to print */
X    int len;        /* number of halfwords to type */
X    int num1, num2;        /* numbers to type */
X    int rem;        /* remainder number of halfwords */
X    char *str;
X
X    if (width) {
X        math_divertio();
X        zprinto(z, 0L);
X        str = math_getdivertedio();
X        math_fill(str, width);
X        free(str);
X        return;
X    }
X    if (zisneg(z))
X        PUTCHAR('-');
X    len = z.len;
X    if ((len == 1) && (*z.v <= (FULL) 7)) {
X        num1 = '0' + *z.v;
X        PUTCHAR(num1);
X        return;
X    }
X    hp = z.v + len - 1;
X    rem = len % 3;
X    switch (rem) {    /* handle odd amounts first */
X        case 0:
X            num1 = (((FULL) hp[0]) << 8) + (((FULL) hp[-1]) >> 8);
X            num2 = (((FULL) (hp[-1] & 0xff)) << 16) + ((FULL) hp[-2]);
X            rem = 3;
X            break;
X        case 1:
X            num1 = 0;
X            num2 = (FULL) hp[0];
X            break;
X        case 2:
X            num1 = (((FULL) hp[0]) >> 8);
X            num2 = (((FULL) (hp[0] & 0xff)) << 16) + ((FULL) hp[-1]);
X            break;
X    }
X    if (num1)
X        PRINTF2("0%o%08o", num1, num2);
X    else
X        PRINTF1("0%o", num2);
X    len -= rem;
X    hp -= rem;
X    while (len > 0) {    /* finish in groups of 3 halfwords */
X        num1 = (((FULL) hp[0]) << 8) + (((FULL) hp[-1]) >> 8);
X        num2 = (((FULL) (hp[-1] & 0xff)) << 16) + ((FULL) hp[-2]);
X        PRINTF2("%08o%08o", num1, num2);
X        hp -= 3;
X        len -= 3;
X    }
X}
X
X
X/*
X * Print a decimal integer to the terminal.
X * This works by dividing the number by 10^2^N for some N, and
X * then doing this recursively on the quotient and remainder.
X * Decimals supplies number of decimal places to print, with a decimal
X * point at the right location, with zero meaning no decimal point.
X * Width is the number of columns to print the number in, including the
X * decimal point and sign if required.  If zero, no extra output is done.
X * If positive, leading spaces are typed if necessary. If negative, trailing
X * spaces are typed if necessary.  As examples of the effects of these values,
X * (345,0,0) = "345", (345,2,0) = "3.45", (345,5,8) = "  .00345".
X */
Xvoid
Xzprintval(z, decimals, width)
X    ZVALUE z;        /* number to be printed */
X    long decimals;        /* number of decimal places */
X    long width;        /* number of columns to print in */
X{
X    int depth;        /* maximum depth */
X    int n;            /* current index into array */
X    int i;            /* number to print */
X    long leadspaces;    /* number of leading spaces to print */
X    long putpoint;        /* digits until print decimal point */
X    long digits;        /* number of digits of raw number */
X    BOOL output;        /* TRUE if have output something */
X    BOOL neg;        /* TRUE if negative */
X    ZVALUE quo, rem;    /* quotient and remainder */
X    ZVALUE leftnums[32];    /* left parts of the number */
X    ZVALUE rightnums[32];    /* right parts of the number */
X
X    if (decimals < 0)
X        decimals = 0;
X    if (width < 0)
X        width = 0;
X    neg = (z.sign != 0);
X
X    leadspaces = width - neg - (decimals > 0);
X    z.sign = 0;
X    /*
X     * Find the 2^N power of ten which is greater than or equal
X     * to the number, calculating it the first time if necessary.
X     */
X    _tenpowers_[0] = _ten_;
X    depth = 0;
X    while ((_tenpowers_[depth].len < z.len) || (zrel(_tenpowers_[depth], z) <= 0)) {
X        depth++;
X        if (_tenpowers_[depth].len == 0)
X            zsquare(_tenpowers_[depth-1], &_tenpowers_[depth]);
X    }
X    /*
X     * Divide by smaller 2^N powers of ten until the parts are small
X     * enough to output.  This algorithm walks through a binary tree
X     * where each node is a piece of the number to print, and such that
X     * we visit left nodes first.  We do the needed recursion in line.
X     */
X    digits = 1;
X    output = FALSE;
X    n = 0;
X    putpoint = 0;
X    rightnums[0].len = 0;
X    leftnums[0] = z;
X    for (;;) {
X        while (n < depth) {
X            i = depth - n - 1;
X            zdiv(leftnums[n], _tenpowers_[i], &quo, &rem);
X            if (!ziszero(quo))
X                digits += (1L << i);
X            n++;
X            leftnums[n] = quo;
X            rightnums[n] = rem;
X        }
X        i = leftnums[n].v[0];
X        if (output || i || (n == 0)) {
X            if (!output) {
X                output = TRUE;
X                if (decimals > digits)
X                    leadspaces -= decimals;
X                else
X                    leadspaces -= digits;
X                while (--leadspaces >= 0)
X                    PUTCHAR(' ');
X                if (neg)
X                    PUTCHAR('-');
X                if (decimals) {
X                    putpoint = (digits - decimals);
X                    if (putpoint <= 0) {
X                        PUTCHAR('.');
X                        while (++putpoint <= 0)
X                            PUTCHAR('0');
X                        putpoint = 0;
X                    }
X                }
X            }
X            i += '0';
X            PUTCHAR(i);
X            if (--putpoint == 0)
X                PUTCHAR('.');
X        }
X        while (rightnums[n].len == 0) {
X            if (n <= 0)
X                return;
X            if (leftnums[n].len)
X                zfree(leftnums[n]);
X            n--;
X        }
X        zfree(leftnums[n]);
X        leftnums[n] = rightnums[n];
X        rightnums[n].len = 0;
X    }
X}
X
X
X/*
X * Read an integer value in decimal, hex, octal, or binary.
X * Hex numbers are indicated by a leading "0x", binary with a leading "0b",
X * and octal by a leading "0".  Periods are skipped over, but any other
X * extraneous character stops the scan.
X */
Xvoid
Xatoz(s, res)
X    register char *s;
X    ZVALUE *res;
X{
X    ZVALUE z, ztmp, digit;
X    HALF digval;
X    BOOL minus;
X    long shift;
X
X    minus = FALSE;
X    shift = 0;
X    if (*s == '+')
X        s++;
X    else if (*s == '-') {
X        minus = TRUE;
X        s++;
X    }
X    if (*s == '0') {        /* possibly hex, octal, or binary */
X        s++;
X        if ((*s >= '0') && (*s <= '7')) {
X            shift = 3;
X        } else if ((*s == 'x') || (*s == 'X')) {
X            shift = 4;
X            s++;
X        } else if ((*s == 'b') || (*s == 'B')) {
X            shift = 1;
X            s++;
X        }
X    }
X    digit.v = &digval;
X    digit.len = 1;
X    digit.sign = 0;
X    z = _zero_;
X    while (*s) {
X        digval = *s++;
X        if ((digval >= '0') && (digval <= '9'))
X            digval -= '0';
X        else if ((digval >= 'a') && (digval <= 'f') && shift)
X            digval -= ('a' - 10);
X        else if ((digval >= 'A') && (digval <= 'F') && shift)
X            digval -= ('A' - 10);
X        else if (digval == '.')
X            continue;
X        else
X            break;
X        if (shift)
X            zshift(z, shift, &ztmp);
X        else
X            zmuli(z, 10L, &ztmp);
X        zfree(z);
X        zadd(ztmp, digit, &z);
X        zfree(ztmp);
X    }
X    ztrim(&z);
X    if (minus && !ziszero(z))
X        z.sign = 1;
X    *res = z;
X}
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/zio.c || echo "restore of calc2.9.0/zio.c fails"
set `wc -c calc2.9.0/zio.c`;Sum=$1
if test "$Sum" != "14342"
then echo original size 14342, current size $Sum;fi
echo "x - extracting calc2.9.0/zmath.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/zmath.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Extended precision integral arithmetic primitives
X */
X
X#include "zmath.h"
X
X
XHALF _twoval_[] = { 2 };
XHALF _zeroval_[] = { 0 };
XHALF _oneval_[] = { 1 };
XHALF _tenval_[] = { 10 };
X
XZVALUE _zero_ = { _zeroval_, 1, 0};
XZVALUE _one_ = { _oneval_, 1, 0 };
XZVALUE _ten_ = { _tenval_, 1, 0 };
X
X
X/*
X * mask of given bits, rotated thru all bit positions twice
X *
X * bitmask[i]      (1 << (i-1)),  for  -BASEB*4<=i<=BASEB*4
X * rotmask[j][k] (1 << ((j+k-1)%BASEB)),  for  -BASEB*2<=j,k<=BASEB*2
X */
Xstatic HALF *bmask;        /* actual rotation thru 8 cycles */
Xstatic HALF **rmask;        /* actual rotation pointers thru 2 cycles */
XHALF *bitmask;            /* bit rotation, norm 0 */
X#if 0
Xstatic HALF **rotmask;        /* pointers to bit rotations, norm index */
X#endif
X
XBOOL _math_abort_;        /* nonzero to abort calculations */
X
X
Xstatic void dadd MATH_PROTO((ZVALUE z1, ZVALUE z2, long y, long n));
Xstatic BOOL dsub MATH_PROTO((ZVALUE z1, ZVALUE z2, long y, long n));
Xstatic void dmul MATH_PROTO((ZVALUE z, FULL x, ZVALUE *dest));
X
X
X#ifdef ALLOCTEST
Xstatic long nalloc = 0;
Xstatic long nfree = 0;
X#endif
X
X
XHALF *
Xalloc(len)
X    long len;
X{
X    HALF *hp;
X
X    if (_math_abort_)
X        math_error("Calculation aborted");
X    hp = (HALF *) malloc(len * sizeof(HALF) + 2);
X    if (hp == 0)
X        math_error("Not enough memory");
X#ifdef ALLOCTEST
X    ++nalloc;
X#endif
X    return hp;
X}
X
X
X#ifdef ALLOCTEST
Xvoid
Xfreeh(h)
X    HALF *h;
X{
X    if ((h != _zeroval_) && (h != _oneval_)) {
X        free(h);
X        ++nfree;
X    }
X}
X
X
Xvoid
XallocStat()
X{
X    fprintf(stderr, "nalloc: %ld nfree: %ld kept: %ld\n",
X        nalloc, nfree, nalloc - nfree);
X}
X#endif
X
X
X/*
X * Convert a normal integer to a number.
X */
Xvoid
Xitoz(i, res)
X    long i;
X    ZVALUE *res;
X{
X    long diddle, len;
X
X    res->len = 1;
X    res->sign = 0;
X    diddle = 0;
X    if (i == 0) {
X        res->v = _zeroval_;
X        return;
X    }
X    if (i < 0) {
X        res->sign = 1;
X        i = -i;
X        if (i < 0) {    /* fix most negative number */
X            diddle = 1;
X            i--;
X        }
X    }
X    if (i == 1) {
X        res->v = _oneval_;
X        return;
X    }
X    len = 1 + (((FULL) i) >= BASE);
X    res->len = len;
X    res->v = alloc(len);
X    res->v[0] = (HALF) (i + diddle);
X    if (len == 2)
X        res->v[1] = (HALF) (i / BASE);
X}
X
X
X/*
X * Convert a number to a normal integer, as far as possible.
X * If the number is out of range, the largest number is returned.
X */
Xlong
Xztoi(z)
X    ZVALUE z;
X{
X    long i;
X
X    if (zisbig(z)) {
X        i = MAXFULL;
X        return (z.sign ? -i : i);
X    }
X    i = (zistiny(z) ? z1tol(z) : z2tol(z));
X    return (z.sign ? -i : i);
X}
X
X
X/*
X * Make a copy of an integer value
X */
Xvoid
Xzcopy(z, res)
X    ZVALUE z, *res;
X{
X    res->sign = z.sign;
X    res->len = z.len;
X    if (zisleone(z)) {    /* zero or plus or minus one are easy */
X        res->v = (z.v[0] ? _oneval_ : _zeroval_);
X        return;
X    }
X    res->v = alloc(z.len);
X    zcopyval(z, *res);
X}
X
X
X/*
X * Add together two integers.
X */
Xvoid
Xzadd(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE dest;
X    HALF *p1, *p2, *pd;
X    long len;
X    FULL carry;
X    SIUNION sival;
X
X    if (z1.sign && !z2.sign) {
X        z1.sign = 0;
X        zsub(z2, z1, res);
X        return;
X    }
X    if (z2.sign && !z1.sign) {
X        z2.sign = 0;
X        zsub(z1, z2, res);
X        return;
X    }
X    if (z2.len > z1.len) {
X        pd = z1.v; z1.v = z2.v; z2.v = pd;
X        len = z1.len; z1.len = z2.len; z2.len = len;
X    }
X    dest.len = z1.len + 1;
X    dest.v = alloc(dest.len);
X    dest.sign = z1.sign;
X    carry = 0;
X    pd = dest.v;
X    p1 = z1.v;
X    p2 = z2.v;
X    len = z2.len;
X    while (len--) {
X        sival.ivalue = ((FULL) *p1++) + ((FULL) *p2++) + carry;
X        *pd++ = sival.silow;
X        carry = sival.sihigh;
X    }
X    len = z1.len - z2.len;
X    while (len--) {
X        sival.ivalue = ((FULL) *p1++) + carry;
X        *pd++ = sival.silow;
X        carry = sival.sihigh;
X    }
X    *pd = (HALF)carry;
X    zquicktrim(dest);
X    *res = dest;
X}
X
X
X/*
X * Subtract two integers.
X */
Xvoid
Xzsub(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    register HALF *h1, *h2, *hd;
X    long len1, len2;
X    FULL carry;
X    SIUNION sival;
X    ZVALUE dest;
X
X    if (z1.sign != z2.sign) {
X        z2.sign = z1.sign;
X        zadd(z1, z2, res);
X        return;
X    }
X    len1 = z1.len;
X    len2 = z2.len;
X    if (len1 == len2) {
X        h1 = z1.v + len1 - 1;
X        h2 = z2.v + len2 - 1;
X        while ((len1 > 0) && ((FULL)*h1 == (FULL)*h2)) {
X            len1--;
X            h1--;
X            h2--;
X        }
X        if (len1 == 0) {
X            *res = _zero_;
X            return;
X        }
X        len2 = len1;
X    }
X    dest.sign = z1.sign;
X    carry = ((len1 < len2) || ((len1 == len2) && ((FULL)*h1 < (FULL)*h2)));
X    h1 = z1.v;
X    h2 = z2.v;
X    if (carry) {
X        carry = len1;
X        len1 = len2;
X        len2 = carry;
X        h1 = z2.v;
X        h2 = z1.v;
X        dest.sign = !dest.sign;
X    }
X    hd = alloc(len1);
X    dest.v = hd;
X    dest.len = len1;
X    len1 -= len2;
X    carry = 0;
X    while (--len2 >= 0) {
X        sival.ivalue = (BASE1 - ((FULL) *h1++)) + *h2++ + carry;
X        *hd++ = BASE1 - sival.silow;
X        carry = sival.sihigh;
X    }
X    while (--len1 >= 0) {
X        sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
X        *hd++ = BASE1 - sival.silow;
X        carry = sival.sihigh;
X    }
X    if (hd[-1] == 0)
X        ztrim(&dest);
X    *res = dest;
X}
X
X
X#if 0
X/*
X * Multiply two integers together.
X */
Xvoid
Xzmul(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    register HALF *s1, *s2, *sd, *h1;
X    FULL mul, carry;
X    long len1, len2;
X    SIUNION sival;
X    ZVALUE dest;
X
X    if (ziszero(z1) || ziszero(z2)) {
X        *res = _zero_;
X        return;
X    }
X    if (zisone(z1)) {
X        zcopy(z2, res);
X        return;
X    }
X    if (zisone(z2)) {
X        zcopy(z1, res);
X        return;
X    }
X    dest.len = z1.len + z2.len;
X    dest.v = alloc(dest.len);
X    dest.sign = (z1.sign != z2.sign);
X    zclearval(dest);
X    /*
X     * Swap the numbers if necessary to make the second one smaller.
X     */
X    if (z1.len < z2.len) {
X        len1 = z1.len;
X        z1.len = z2.len;
X        z2.len = len1;
X        s1 = z1.v;
X        z1.v = z2.v;
X        z2.v = s1;
X    }
X    /*
X     * Multiply the first number by each 'digit' of the second number
X     * and add the result to the total.
X     */
X    sd = dest.v;
X    s2 = z2.v;
X    for (len2 = z2.len; len2--; sd++) {
X        mul = (FULL) *s2++;
X        if (mul == 0)
X            continue;
X        h1 = sd;
X        s1 = z1.v;
X        carry = 0;
X        len1 = z1.len;
X        while (len1 >= 4) {    /* expand the loop some */
X            len1 -= 4;
X            sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + carry;
X            *h1++ = sival.silow;
X            sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + ((FULL) sival.sihigh);
X            *h1++ = sival.silow;
X            sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + ((FULL) sival.sihigh);
X            *h1++ = sival.silow;
X            sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + ((FULL) sival.sihigh);
X            *h1++ = sival.silow;
X            carry = sival.sihigh;
X        }
X        while (--len1 >= 0) {
X            sival.ivalue = (mul * ((FULL) *s1++)) + ((FULL) *h1) + carry;
X            *h1++ = sival.silow;
X            carry = sival.sihigh;
X        }
X        while (carry) {
X            sival.ivalue = ((FULL) *h1) + carry;
X            *h1++ = sival.silow;
X            carry = sival.sihigh;
X        }
X    }
X    ztrim(&dest);
X    *res = dest;
X}
X#endif
X
X
X/*
X * Multiply an integer by a small number.
X */
Xvoid
Xzmuli(z, n, res)
X    ZVALUE z;
X    long n;
X    ZVALUE *res;
X{
X    register HALF *h1, *sd;
X    FULL low;
X    FULL high;
X    FULL carry;
X    long len;
X    SIUNION sival;
X    ZVALUE dest;
X
X    if ((n == 0) || ziszero(z)) {
X        *res = _zero_;
X        return;
X    }
X    if (n < 0) {
X        n = -n;
X        z.sign = !z.sign;
X    }
X    if (n == 1) {
X        zcopy(z, res);
X        return;
X    }
X    low = ((FULL) n) & BASE1;
X    high = ((FULL) n) >> BASEB;
X    dest.len = z.len + 2;
X    dest.v = alloc(dest.len);
X    dest.sign = z.sign;
X    /*
X     * Multiply by the low digit.
X     */
X    h1 = z.v;
X    sd = dest.v;
X    len = z.len;
X    carry = 0;
X    while (len--) {
X        sival.ivalue = ((FULL) *h1++) * low + carry;
X        *sd++ = sival.silow;
X        carry = sival.sihigh;
X    }
X    *sd = (HALF)carry;
X    /*
X     * If there was only one digit, then we are all done except
X     * for trimming the number if there was no last carry.
X     */
X    if (high == 0) {
X        dest.len--;
X        if (carry == 0)
X            dest.len--;
X        *res = dest;
X        return;
X    }
X    /*
X     * Need to multiply by the high digit and add it into the
X     * previous value.  Clear the final word of rubbish first.
X     */
X    *(++sd) = 0;
X    h1 = z.v;
X    sd = dest.v + 1;
X    len = z.len;
X    carry = 0;
X    while (len--) {
X        sival.ivalue = ((FULL) *h1++) * high + ((FULL) *sd) + carry;
X        *sd++ = sival.silow;
X        carry = sival.sihigh;
X    }
X    *sd = (HALF)carry;
X    zquicktrim(dest);
X    *res = dest;
X}
X
X
X/*
X * Divide two numbers by their greatest common divisor.
X * This is useful for reducing the numerator and denominator of
X * a fraction to its lowest terms.
X */
Xvoid
Xzreduce(z1, z2, z1res, z2res)
X    ZVALUE z1, z2, *z1res, *z2res;
X{
X    ZVALUE tmp;
X
X    if (zisleone(z1) || zisleone(z2))
X        tmp = _one_;
X    else
X        zgcd(z1, z2, &tmp);
X    if (zisunit(tmp)) {
X        zcopy(z1, z1res);
X        zcopy(z2, z2res);
X    } else {
X        zquo(z1, tmp, z1res);
X        zquo(z2, tmp, z2res);
X    }
X    zfree(tmp);
X}
X
X
X/*
X * Divide two numbers to obtain a quotient and remainder.
X * This algorithm is taken from
X * Knuth, The Art of Computer Programming, vol 2: Seminumerical Algorithms.
X * Slight modifications were made to speed this mess up.
X */
Xvoid
Xzdiv(z1, z2, res, rem)
X    ZVALUE z1, z2, *res, *rem;
X{
X    long i, j, k;
X    register HALF *q, *pp;
X    SIUNION pair;        /* pair of halfword values */
X    HALF h2, v2;
X    long y;
X    FULL x;
X    ZVALUE ztmp1, ztmp2, ztmp3, quo;
X
X    if (ziszero(z2))
X        math_error("Division by zero");
X    if (ziszero(z1)) {
X        *res = _zero_;
X        *rem = _zero_;
X        return;
X    }
X    if (zisone(z2)) {
X        zcopy(z1, res);
X        *rem = _zero_;
X        return;
X    }
X    i = 32768;
X    j = 0;
X    k = (long) z2.v[z2.len - 1];
X    while (! (k & i)) {
X        j ++;
X        i >>= 1;
X    }
X    ztmp1.v = alloc(z1.len + 1);
X    ztmp1.len = z1.len + 1;
X    zcopyval(z1, ztmp1);
X    ztmp1.v[z1.len] = 0;
X    ztmp1.sign = 0;
X    ztmp2.v = alloc(z2.len);
X    ztmp2.len = z2.len;
X    ztmp2.sign = 0;
X    zcopyval(z2, ztmp2);
X    if (zrel(ztmp1, ztmp2) < 0) {
X        rem->v = ztmp1.v;
X        rem->sign = z1.sign;
X        rem->len = z1.len;
X        zfree(ztmp2);
X        *res = _zero_;
X        return;
X    }
X    quo.len = z1.len - z2.len + 1;
X    quo.v = alloc(quo.len);
X    quo.sign = z1.sign != z2.sign;
X    zclearval(quo);
X
X    ztmp3.v = zalloctemp(z2.len + 1);
X
X    /*
X     * Normalize z1 and z2
X     */
X    zshiftl(ztmp1, j);
X    zshiftl(ztmp2, j);
X
X    k = ztmp1.len - ztmp2.len;
X    q = quo.v + quo.len;
X    y = ztmp1.len - 1;
X    h2 = ztmp2.v [ztmp2.len - 1];
X    v2 = 0;
X    if (ztmp2.len >= 2)
X        v2 = ztmp2.v [ztmp2.len - 2];
X    for (;k--; --y) {
X        pp = ztmp1.v + y - 1;
X        pair.silow = pp[0];
X        pair.sihigh = pp[1];
X        if (ztmp1.v[y] == h2)
X            x = BASE1;
X        else
X            x = pair.ivalue / h2;
X        if (x) {
X            while (pair.ivalue - x * h2 < BASE &&
X                v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
X                    --x;
X            }
X            dmul(ztmp2, x, &ztmp3);
X#ifdef divblab
X            printf(" x: %ld\n", x);
X            printf("ztmp1: ");
X            printz(ztmp1);
X            printf("ztmp2: ");
X            printz(ztmp2);
X            printf("ztmp3: ");
X            printz(ztmp3);
X#endif
X            if (dsub(ztmp1, ztmp3, y, ztmp2.len)) {
X                --x;
X                /*
X                printf("adding back\n");
X                */
X                dadd(ztmp1, ztmp2, y, ztmp2.len);
X            }
X        }
X        ztrim(&ztmp1);
X        *--q = (HALF)x;
X    }
X    zshiftr(ztmp1, j);
X    *rem = ztmp1;
X    ztrim(rem);
X    zfree(ztmp2);
X    ztrim(&quo);
X    *res = quo;
X}
X
X
X/*
X * Return the quotient and remainder of an integer divided by a small
X * number.  A nonzero remainder is only meaningful when both numbers
X * are positive.
X */
Xlong
Xzdivi(z, n, res)
X    ZVALUE z, *res;
X    long n;
X{
X    register HALF *h1, *sd;
X    FULL val;
X    HALF divval[2];
X    ZVALUE div;
X    ZVALUE dest;
X    long len;
X
X    if (n == 0)
X        math_error("Division by zero");
X    if (ziszero(z)) {
X        *res = _zero_;
X        return 0;
X    }
X    if (n < 0) {
X        n = -n;
X        z.sign = !z.sign;
X    }
X    if (n == 1) {
X        zcopy(z, res);
X        return 0;
X    }
X    /*
X     * If the division is by a large number, then call the normal
X     * divide routine.
X     */
X    if (n & ~BASE1) {
X        div.sign = 0;
X        div.len = 2;
X        div.v = divval;
X        divval[0] = (HALF) n;
X        divval[1] = ((FULL) n) >> BASEB;
X        zdiv(z, div, res, &dest);
X        n = (zistiny(dest) ? z1tol(dest) : z2tol(dest));
X        zfree(dest);
X        return n;
X    }
X    /*
X     * Division is by a small number, so we can be quick about it.
X     */
X    len = z.len;
X    dest.sign = z.sign;
X    dest.len = len;
X    dest.v = alloc(len);
X    h1 = z.v + len - 1;
X    sd = dest.v + len - 1;
X    val = 0;
X    while (len--) {
X        val = ((val << BASEB) + ((FULL) *h1--));
X        *sd-- = val / n;
X        val %= n;
X    }
X    zquicktrim(dest);
X    *res = dest;
X    return val;
X}
X
X
X/*
X * Return the quotient of two numbers.
X * This works the same as zdiv, except that the remainer is not returned.
X */
Xvoid
Xzquo(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    long i, j, k;
X    register HALF *q, *pp;
X    SIUNION pair;            /* pair of halfword values */
X    HALF h2, v2;
X    long y;
X    FULL x;
X    ZVALUE ztmp1, ztmp2, ztmp3, quo;
X
X    if (ziszero(z2))
X        math_error("Division by zero");
X    if (ziszero(z1)) {
X        *res = _zero_;
X        return;
X    }
X    if (zisone(z2)) {
X        zcopy(z1, res);
X        return;
X    }
X    i = 32768;
X    j = 0;
X    k = (long) z2.v[z2.len - 1];
X    while (! (k & i)) {
X        j ++;
X        i >>= 1;
X    }
X    ztmp1.v = alloc(z1.len + 1);
X    ztmp1.len = z1.len + 1;
X    zcopyval(z1, ztmp1);
X    ztmp1.v[z1.len] = 0;
X    ztmp1.sign = 0;
X    ztmp2.v = alloc(z2.len);
X    ztmp2.len = z2.len;
X    ztmp2.sign = 0;
X    zcopyval(z2, ztmp2);
X    if (zrel(ztmp1, ztmp2) < 0) {
X        zfree(ztmp1);
X        zfree(ztmp2);
X        *res = _zero_;
X        return;
X    }
X    quo.len = z1.len - z2.len + 1;
X    quo.v = alloc(quo.len);
X    quo.sign = z1.sign != z2.sign;
X    zclearval(quo);
X
X    ztmp3.v = zalloctemp(z2.len + 1);
X
X    /*
X     * Normalize z1 and z2
X     */
X    zshiftl(ztmp1, j);
X    zshiftl(ztmp2, j);
X
X    k = ztmp1.len - ztmp2.len;
X    q = quo.v + quo.len;
X    y = ztmp1.len - 1;
X    h2 = ztmp2.v [ztmp2.len - 1];
X    v2 = 0;
X    if (ztmp2.len >= 2)
X        v2 = ztmp2.v [ztmp2.len - 2];
X    for (;k--; --y) {
X        pp = ztmp1.v + y - 1;
X        pair.silow = pp[0];
X        pair.sihigh = pp[1];
X        if (ztmp1.v[y] == h2)
X            x = BASE1;
X        else
X            x = pair.ivalue / h2;
X        if (x) {
X            while (pair.ivalue - x * h2 < BASE &&
X                v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
X                    --x;
X            }
X            dmul(ztmp2, x, &ztmp3);
X            if (dsub(ztmp1, ztmp3, y, ztmp2.len)) {
X                --x;
X                dadd(ztmp1, ztmp2, y, ztmp2.len);
X            }
X        }
X        ztrim(&ztmp1);
X        *--q = (HALF)x;
X    }
X    zfree(ztmp1);
X    zfree(ztmp2);
X    ztrim(&quo);
X    *res = quo;
X}
X
X
X/*
X * Compute the remainder after dividing one number by another.
X * This is only defined for positive z2 values.
X * The result is normalized to lie in the range 0 to z2-1.
X */
Xvoid
Xzmod(z1, z2, rem)
X    ZVALUE z1, z2, *rem;
X{
X    long i, j, k, neg;
X    register HALF *pp;
X    SIUNION pair;            /* pair of halfword values */
X    HALF h2, v2;
X    long y;
X    FULL x;
X    ZVALUE ztmp1, ztmp2, ztmp3;
X
X    if (ziszero(z2))
X        math_error("Division by zero");
X    if (zisneg(z2))
X        math_error("Non-positive modulus");
X    if (ziszero(z1) || zisunit(z2)) {
X        *rem = _zero_;
X        return;
X    }
X    if (zistwo(z2)) {
X        if (zisodd(z1))
X            *rem = _one_;
X        else
X            *rem = _zero_;
X        return;
X    }
X    neg = z1.sign;
X    z1.sign = 0;
X
X    /*
X     * Do a quick check to see if the absolute value of the number
X     * is less than the modulus.  If so, then the result is just a
X     * subtract or a copy.
X     */
X    h2 = z1.v[z1.len - 1];
X    v2 = z2.v[z2.len - 1];
X    if ((z1.len < z2.len) || ((z1.len == z2.len) && (h2 < v2))) {
X        if (neg)
X            zsub(z2, z1, rem);
X        else
X            zcopy(z1, rem);
X        return;
X    }
X
X    /*
X     * Do another quick check to see if the number is positive and
X     * between the size of the modulus and twice the modulus.
X     * If so, then the answer is just another subtract.
X     */
X    if (!neg && (z1.len == z2.len) && (h2 > v2) &&
X        (((FULL) h2) < 2 * ((FULL) v2)))
X    {
X        zsub(z1, z2, rem);
X        return;
X    }
X
X    /*
X     * If the modulus is an exact power of two, then the result
X     * can be obtained by ignoring the high bits of the number.
X     * This truncation assumes that the number of words for the
X     * number is at least as large as the number of words in the
X     * modulus, which is true at this point.
X     */
X    if (((v2 & -v2) == v2) && zisonebit(z2)) {    /* ASSUMES 2'S COMP */
X        i = zhighbit(z2);
X        z1.len = (i + BASEB - 1) / BASEB;
X        zcopy(z1, &ztmp1);
X        i %= BASEB;
X        if (i)
X            ztmp1.v[ztmp1.len - 1] &= ((((HALF) 1) << i) - 1);
X        ztmp2.len = 0;
X        goto gotanswer;
X    }
X
X    /*
X     * If the modulus is one less than an exact power of two, then
X     * the result can be simplified similarly to "casting out 9's".
X     * Only do this simplification for large enough modulos.
X     */
X    if ((z2.len > 1) && (z2.v[0] == BASE1) && zisallbits(z2)) {
X        i = -(zhighbit(z2) + 1);
X        zcopy(z1, &ztmp1);
X        z1 = ztmp1;
X        while ((k = zrel(z1, z2)) > 0) {
X            ztmp1 = _zero_;
X            while (!ziszero(z1)) {
X                zand(z1, z2, &ztmp2);
X                zadd(ztmp2, ztmp1, &ztmp3);
X                zfree(ztmp1);
X                zfree(ztmp2);
X                ztmp1 = ztmp3;
X                zshift(z1, i, &ztmp2);
X                zfree(z1);
X                z1 = ztmp2;
X            }
X            zfree(z1);
X            z1 = ztmp1;
X        }
X        if (k == 0) {
X            zfree(ztmp1);
X            *rem = _zero_;
X            return;
X        }
X        ztmp2.len = 0;
X        goto gotanswer;
X    }
X
X    /*
X     * Must actually do the divide.
X     */
X    i = 32768;
X    j = 0;
X    k = (long) z2.v[z2.len - 1];
X    while (! (k & i)) {
X        j ++;
X        i >>= 1;
X    }
X    ztmp1.v = alloc(z1.len + 1);
X    ztmp1.len = z1.len + 1;
X    zcopyval(z1, ztmp1);
X    ztmp1.v[z1.len] = 0;
X    ztmp1.sign = 0;
X    ztmp2.v = alloc(z2.len);
X    ztmp2.len = z2.len;
X    ztmp2.sign = 0;
X    zcopyval(z2, ztmp2);
X    if (zrel(ztmp1, ztmp2) < 0)
X        goto gotanswer;
X
X    ztmp3.v = zalloctemp(z2.len + 1);
X
X    /*
X     * Normalize z1 and z2
X     */
X    zshiftl(ztmp1, j);
X    zshiftl(ztmp2, j);
X
X    k = ztmp1.len - ztmp2.len;
X    y = ztmp1.len - 1;
X    h2 = ztmp2.v [ztmp2.len - 1];
X    v2 = 0;
X    if (ztmp2.len >= 2)
X        v2 = ztmp2.v [ztmp2.len - 2];
X    for (;k--; --y) {
X        pp = ztmp1.v + y - 1;
X        pair.silow = pp[0];
X        pair.sihigh = pp[1];
X        if (ztmp1.v[y] == h2)
X            x = BASE1;
X        else
X            x = pair.ivalue / h2;
X        if (x) {
X            while (pair.ivalue - x * h2 < BASE &&
X                v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
X                    --x;
X            }
X            dmul(ztmp2, x, &ztmp3);
X            if (dsub(ztmp1, ztmp3, y, ztmp2.len))
X                dadd(ztmp1, ztmp2, y, ztmp2.len);
X        }
X        ztrim(&ztmp1);
X    }
X    zshiftr(ztmp1, j);
X
Xgotanswer:
X    ztrim(&ztmp1);
X    if (ztmp2.len)
X        zfree(ztmp2);
X    if (neg && !ziszero(ztmp1)) {
X        zsub(z2, ztmp1, rem);
X        zfree(ztmp1);
X    } else
X        *rem = ztmp1;
X}
X
X
X/*
X * Calculate the mod of an integer by a small number.
X * This is only defined for positive moduli.
X */
Xlong
Xzmodi(z, n)
X    ZVALUE z;
X    long n;
X{
X    register HALF *h1;
X    FULL val;
X    HALF divval[2];
X    ZVALUE div;
X    ZVALUE temp;
X    long len;
X
X    if (n == 0)
X        math_error("Division by zero");
X    if (n < 0)
X        math_error("Non-positive modulus");
X    if (ziszero(z) || (n == 1))
X        return 0;
X    if (zisone(z))
X        return 1;
X    /*
X     * If the modulus is by a large number, then call the normal
X     * modulo routine.
X     */
X    if (n & ~BASE1) {
X        div.sign = 0;
X        div.len = 2;
X        div.v = divval;
X        divval[0] = (HALF) n;
X        divval[1] = ((FULL) n) >> BASEB;
X        zmod(z, div, &temp);
X        n = (zistiny(temp) ? z1tol(temp) : z2tol(temp));
X        zfree(temp);
X        return n;
X    }
X    /*
X     * The modulus is by a small number, so we can do this quickly.
X     */
X    len = z.len;
X    h1 = z.v + len - 1;
X    val = 0;
X    while (len--)
X        val = ((val << BASEB) + ((FULL) *h1--)) % n;
X    if (z.sign)
X        val = n - val;
X    return val;
X}
X
X
X/*
X * Return whether or not one number exactly divides another one.
X * Returns TRUE if division occurs with no remainder.
X * z1 is the number to be divided by z2.
X */
XBOOL
Xzdivides(z1, z2)
X    ZVALUE z1, z2;        /* numbers to test division into and by */
X{
X    ZVALUE temp;
X    long cv;
X
X    z1.sign = 0;
X    z2.sign = 0;
X    /*
X     * Take care of obvious cases first
X     */
X    if (zisleone(z2)) {    /* division by zero or one */
X        if (*z2.v == 0)
X            math_error("Division by zero");
X        return TRUE;
X    }
X    if (ziszero(z1))    /* everything divides zero */
X        return TRUE;
X    if (z1.len < z2.len)    /* quick size comparison */
X        return FALSE;
X    if ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))    /* more */
X        return FALSE;
X    if (zisodd(z1) && ziseven(z2))    /* can't divide odd by even */
X        return FALSE;
X    if (zlowbit(z1) < zlowbit(z2))    /* can't have smaller power of two */
X        return FALSE;
X    cv = zrel(z1, z2);    /* can't divide smaller number */
X    if (cv <= 0)
X        return (cv == 0);
X    /*
X     * Now do the real work.  Divisor divides dividend if the gcd of the
X     * two numbers equals the divisor.
X     */
X    zgcd(z1, z2, &temp);
X    cv = zcmp(z2, temp);
X    zfree(temp);
X    return (cv == 0);
X}
X
X
X/*
X * Compute the logical OR of two numbers
X */
Xvoid
Xzor(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    register HALF *sp, *dp;
X    long len;
X    ZVALUE bz, lz, dest;
X
X    if (z1.len >= z2.len) {
X        bz = z1;
X        lz = z2;
X    } else {
X        bz = z2;
X        lz = z1;
X    }
X    dest.len = bz.len;
X    dest.v = alloc(dest.len);
X    dest.sign = 0;
X    zcopyval(bz, dest);
X    len = lz.len;
X    sp = lz.v;
X    dp = dest.v;
X    while (len--)
X        *dp++ |= *sp++;
X    *res = dest;
X}
X
X
X/*
X * Compute the logical AND of two numbers.
X */
Xvoid
Xzand(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    register HALF *h1, *h2, *hd;
X    register long len;
X    ZVALUE dest;
X
X    len = ((z1.len <= z2.len) ? z1.len : z2.len);
X    h1 = &z1.v[len-1];
X    h2 = &z2.v[len-1];
X    while ((len > 1) && ((*h1 & *h2) == 0)) {
X        h1--;
X        h2--;
X        len--;
X    }
X    dest.len = len;
X    dest.v = alloc(len);
X    dest.sign = 0;
X    h1 = z1.v;
X    h2 = z2.v;
X    hd = dest.v;
X    while (len--)
X        *hd++ = (*h1++ & *h2++);
X    *res = dest;
X}
X
X
X/*
X * Compute the logical XOR of two numbers.
X */
Xvoid
Xzxor(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    register HALF *sp, *dp;
X    long len;
X    ZVALUE bz, lz, dest;
X
X    if (z1.len == z2.len) {
X        for (len = z1.len; ((len > 1) && (z1.v[len-1] == z2.v[len-1])); len--) ;
X        z1.len = len;
X        z2.len = len;
X    }
X    if (z1.len >= z2.len) {
X        bz = z1;
X        lz = z2;
X    } else {
X        bz = z2;
X        lz = z1;
X    }
X    dest.len = bz.len;
X    dest.v = alloc(dest.len);
X    dest.sign = 0;
X    zcopyval(bz, dest);
X    len = lz.len;
X    sp = lz.v;
X    dp = dest.v;
X    while (len--)
X        *dp++ ^= *sp++;
X    *res = dest;
X}
X
X
X/*
X * Shift a number left (or right) by the specified number of bits.
X * Positive shift means to the left.  When shifting right, rightmost
X * bits are lost.  The sign of the number is preserved.
X */
Xvoid
Xzshift(z, n, res)
X    ZVALUE z, *res;
X    long n;
X{
X    ZVALUE ans;
X    long hc;        /* number of halfwords shift is by */
X
X    if (ziszero(z)) {
X        *res = _zero_;
X        return;
X    }
X    if (n == 0) {
X        zcopy(z, res);
X        return;
X    }
X    /*
X     * If shift value is negative, then shift right.
X     * Check for large shifts, and handle word-sized shifts quickly.
X     */
X    if (n < 0) {
X        n = -n;
X        if ((n < 0) || (n >= (z.len * BASEB))) {
X            *res = _zero_;
X            return;
X        }
X        hc = n / BASEB;
X        n %= BASEB;
X        z.v += hc;
X        z.len -= hc;
X        ans.len = z.len;
X        ans.v = alloc(ans.len);
X        ans.sign = z.sign;
X        zcopyval(z, ans);
X        if (n > 0) {
X            zshiftr(ans, n);
X            ztrim(&ans);
X        }
X        if (ziszero(ans)) {
X            zfree(ans);
X            ans = _zero_;
X        }
X        *res = ans;
X        return;
X    }
X    /*
X     * Shift value is positive, so shift leftwards.
X     * Check specially for a shift of the value 1, since this is common.
X     * Also handle word-sized shifts quickly.
X     */
X    if (zisunit(z)) {
X        zbitvalue(n, res);
X        res->sign = z.sign;
X        return;
X    }
X    hc = n / BASEB;
X    n %= BASEB;
X    ans.len = z.len + hc + 1;
X    ans.v = alloc(ans.len);
X    ans.sign = z.sign;
X    if (hc > 0)
X        memset((char *) ans.v, 0, hc * sizeof(HALF));
X    memcpy((char *) (ans.v + hc),
X        (char *) z.v, z.len * sizeof(HALF));
X    ans.v[ans.len - 1] = 0;
X    if (n > 0) {
X        ans.v += hc;
X        ans.len -= hc;
X        zshiftl(ans, n);
X        ans.v -= hc;
X        ans.len += hc;
X    }
X    ztrim(&ans);
X    *res = ans;
X}
X
X
X/*
X * Return the position of the lowest bit which is set in the binary
X * representation of a number (counting from zero).  This is the highest
X * power of two which evenly divides the number.
X */
Xlong
Xzlowbit(z)
X    ZVALUE z;
X{
X    register HALF *zp;
X    long n;
X    HALF dataval;
X    HALF *bitval;
X
X    n = 0;
X    for (zp = z.v; *zp == 0; zp++)
X        if (++n >= z.len)
X            return 0;
X    dataval = *zp;
X    bitval = bitmask;
X    while ((*(bitval++) & dataval) == 0) {
X    }
X    return (n*BASEB)+(bitval-bitmask-1);
X}
X
X
X/*
X * Return the position of the highest bit which is set in the binary
X * representation of a number (counting from zero).  This is the highest power
X * of two which is less than or equal to the number (which is assumed nonzero).
X */
Xlong
Xzhighbit(z)
X    ZVALUE z;
X{
X    HALF dataval;
X    HALF *bitval;
X
X    dataval = z.v[z.len-1];
X    if (dataval) {
X        bitval = bitmask+BASEB;
X        while ((*(--bitval) & dataval) == 0) {
X        }
X    }
X    return (z.len*BASEB)+(bitval-bitmask-BASEB);
X}
X
X
X#if 0
X/*
X * Reverse the bits of a particular range of bits of a number.
X *
X * This function returns an integer with bits a thru b swapped.
X * That is, bit a is swapped with bit b, bit a+1 is swapped with b-1,
X * and so on.
X *
X * As a special case, if the ending bit position is < 0, is it taken to
X * mean the highest bit set.  Thus zbitrev(0, -1, z, &res) will
X * perform a complete bit reverse of the number 'z'.
X *
X * As a special case, if the starting bit position is < 0, is it taken to
X * mean the lowest bit set.  Thus zbitrev(-1, -1, z, &res) is the
X * same as zbitrev(lowbit(z), highbit(z), z, &res).
X *
X * Note that the low order bit number is taken to be 0.  Also, bitrev
X * ignores the sign of the number.
X *
X * Bits beyond the highest bit are taken to be zero.  Thus the calling
X * bitrev(0, 100, _one_, &res) will result in a value of 2^100.
X */
Xvoid
Xzbitrev(low, high, z, res)
X    long low;    /* lowest bit to reverse, <0 => lowbit(z) */
X    long high;    /* highest bit to reverse, <0 => highbit(z) */
X    ZVALUE z;    /* value to bit reverse */
X    ZVALUE *res;    /* resulting bit reverse number */
X{
X}
X#endif
X
X
X/*
X * Return whether or not the specifed bit number is set in a number.
X * Rightmost bit of a number is bit 0.
X */
XBOOL
Xzisset(z, n)
X    ZVALUE z;
X    long n;
X{
X    if ((n < 0) || ((n / BASEB) >= z.len))
X        return FALSE;
X    return ((z.v[n / BASEB] & (((HALF) 1) << (n % BASEB))) != 0);
X}
X
X
X/*
X * Check whether or not a number has exactly one bit set, and
X * thus is an exact power of two.  Returns TRUE if so.
X */
XBOOL
Xzisonebit(z)
X    ZVALUE z;
X{
X    register HALF *hp;
X    register LEN len;
X
X    if (ziszero(z) || zisneg(z))
X        return FALSE;
X    hp = z.v;
X    len = z.len;
X    while (len > 4) {
X        len -= 4;
X        if (*hp++ || *hp++ || *hp++ || *hp++)
X            return FALSE;
X    }
X    while (--len > 0) {
X        if (*hp++)
X            return FALSE;
X    }
X    return ((*hp & -*hp) == *hp);        /* NEEDS 2'S COMPLEMENT */
X}
X
X
X/*
X * Check whether or not a number has all of its bits set below some
X * bit position, and thus is one less than an exact power of two.
X * Returns TRUE if so.
X */
XBOOL
Xzisallbits(z)
X    ZVALUE z;
X{
X    register HALF *hp;
X    register LEN len;
X    HALF digit;
X
X    if (ziszero(z) || zisneg(z))
X        return FALSE;
X    hp = z.v;
X    len = z.len;
X    while (len > 4) {
X        len -= 4;
X        if ((*hp++ != BASE1) || (*hp++ != BASE1) ||
X            (*hp++ != BASE1) || (*hp++ != BASE1))
X                return FALSE;
X    }
X    while (--len > 0) {
X        if (*hp++ != BASE1)
X            return FALSE;
X    }
X    digit = *hp + 1;
X    return ((digit & -digit) == digit);    /* NEEDS 2'S COMPLEMENT */
X}
X
X
X/*
X * Return the number whose binary representation contains only one bit which
X * is in the specified position (counting from zero).  This is equivilant
X * to raising two to the given power.
X */
Xvoid
Xzbitvalue(n, res)
X    long n;
X    ZVALUE *res;
X{
X    ZVALUE z;
X
X    if (n < 0) n = 0;
X    z.sign = 0;
X    z.len = (n / BASEB) + 1;
X    z.v = alloc(z.len);
X    zclearval(z);
X    z.v[z.len-1] = (((HALF) 1) << (n % BASEB));
X    *res = z;
X}
X
X
X/*
X * Compare a number against zero.
X * Returns the sgn function of the number (-1, 0, or 1).
X */
XFLAG
Xztest(z)
X    ZVALUE z;
X{
X    register int sign;
X    register HALF *h;
X    register long len;
X
X    sign = 1;
X    if (z.sign)
X        sign = -sign;
X    h = z.v;
X    len = z.len;
X    while (len--)
X        if (*h++)
X            return sign;
X    return 0;
X}
X
X
X/*
X * Compare two numbers to see which is larger.
X * Returns -1 if first number is smaller, 0 if they are equal, and 1 if
X * first number is larger.  This is the same result as ztest(z2-z1).
X */
XFLAG
Xzrel(z1, z2)
X    ZVALUE z1, z2;
X{
X    register HALF *h1, *h2;
X    register long len1, len2;
X    int sign;
X
X    sign = 1;
X    if (z1.sign < z2.sign)
X        return 1;
X    if (z2.sign < z1.sign)
X        return -1;
X    if (z2.sign)
X        sign = -1;
X    len1 = z1.len;
X    len2 = z2.len;
X    h1 = z1.v + z1.len - 1;
X    h2 = z2.v + z2.len - 1;
X    while (len1 > len2) {
X        if (*h1--)
X            return sign;
X        len1--;
X    }
X    while (len2 > len1) {
X        if (*h2--)
X            return -sign;
X        len2--;
X    }
X    while (len1--) {
X        if (*h1-- != *h2--)
X            break;
X    }
X    if ((len1 = *++h1) > (len2 = *++h2))
X        return sign;
X    if (len1 < len2)
X        return -sign;
X    return 0;
X}
X
X
X/*
X * Compare two numbers to see if they are equal or not.
X * Returns TRUE if they differ.
X */
XBOOL
Xzcmp(z1, z2)
X    ZVALUE z1, z2;
X{
X    register HALF *h1, *h2;
X    register long len;
X
X    if ((z1.sign != z2.sign) || (z1.len != z2.len) || (*z1.v != *z2.v))
X        return TRUE;
X    len = z1.len;
X    h1 = z1.v;
X    h2 = z2.v;
X    while (len-- > 0) {
X        if (*h1++ != *h2++)
X            return TRUE;
X    }
X    return FALSE;
X}
X
X
X/*
X * Internal utility subroutines
X */
Xstatic void
Xdadd(z1, z2, y, n)
X    ZVALUE z1, z2;
X    long y, n;
X{
X    HALF *s1, *s2;
X    short carry;
X    long sum;
X
X    s1 = z1.v + y - n;
X    s2 = z2.v;
X    carry = 0;
X    while (n--) {
X        sum = (long)*s1 + (long)*s2 + carry;
X        carry = 0;
X        if (sum >= BASE) {
X            sum -= BASE;
X            carry = 1;
X        }
X        *s1 = (HALF)sum;
X        ++s1;
X        ++s2;
X    }
X    sum = (long)*s1 + carry;
X    *s1 = (HALF)sum;
X}
X
X
X/*
X * Do subtract for divide, returning TRUE if subtraction went negative.
X */
Xstatic BOOL
Xdsub(z1, z2, y, n)
X    ZVALUE z1, z2;
X    long y, n;
X{
X    HALF *s1, *s2, *s3;
X    FULL i1;
X    BOOL neg;
X
X    neg = FALSE;
X    s1 = z1.v + y - n;
X    s2 = z2.v;
X    if (++n > z2.len)
X        n = z2.len;
X    while (n--) {
X        i1 = (FULL) *s1;
X        if (i1 < (FULL) *s2) {
X            s3 = s1 + 1;
X            while (s3 < z1.v + z1.len && !(*s3)) {
X                *s3 = BASE1;
X                ++s3;
X            }
X            if (s3 >= z1.v + z1.len)
X                neg = TRUE;
X            else
X                --(*s3);
X            i1 += BASE;
X        }
X        *s1 = i1 - (FULL) *s2;
X        ++s1;
X        ++s2;
X    }
X    return neg;
X}
X
X
X/*
X * Multiply a number by a single 'digit'.
X * This is meant to be used only by the divide routine, and so the
X * destination area must already be allocated and be large enough.
X */
Xstatic void
Xdmul(z, mul, dest)
X    ZVALUE z;
X    FULL mul;
X    ZVALUE *dest;
X{
X    register HALF *zp, *dp;
X    SIUNION sival;
X    FULL carry;
X    long len;
X
X    dp = dest->v;
X    dest->sign = 0;
X    if (mul == 0) {
X        dest->len = 1;
X        *dp = 0;
X        return;
X    }
X    len = z.len;
X    zp = z.v + len - 1;
X    while ((*zp == 0) && (len > 1)) {
X        len--;
X        zp--;
X    }
X    dest->len = len;
X    zp = z.v;
X    carry = 0;
X    while (len >= 4) {
X        len -= 4;
X        sival.ivalue = (mul * ((FULL) *zp++)) + carry;
X        *dp++ = sival.silow;
X        sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
X        *dp++ = sival.silow;
X        sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
X        *dp++ = sival.silow;
X        sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
X        *dp++ = sival.silow;
X        carry = sival.sihigh;
X    }
X    while (--len >= 0) {
X        sival.ivalue = (mul * ((FULL) *zp++)) + carry;
X        *dp++ = sival.silow;
X        carry = sival.sihigh;
X    }
X    if (carry) {
X        *dp = (HALF)carry;
X        dest->len++;
X    }
X}
X
X
X/*
X * Utility to calculate the gcd of two small integers.
X */
Xlong
Xiigcd(i1, i2)
X    long i1, i2;
X{
X    FULL f1, f2, temp;
X
X    f1 = (FULL) ((i1 >= 0) ? i1 : -i1);
X    f2 = (FULL) ((i2 >= 0) ? i2 : -i2);
X    while (f1) {
X        temp = f2 % f1;
X        f2 = f1;
X        f1 = temp;
X    }
X    return (long) f2;
X}
X
X
Xvoid
Xztrim(z)
X    ZVALUE *z;
X{
X    register HALF *h;
X    register long len;
X
X    h = z->v + z->len - 1;
X    len = z->len;
X    while (*h == 0 && len > 1) {
X        --h;
X        --len;
X    }
X    z->len = len;
X}
X
X
X/*
X * Utility routine to shift right.
X */
Xvoid
Xzshiftr(z, n)
X    ZVALUE z;
X    long n;
X{
X    register HALF *h, *lim;
X    FULL mask, maskt;
X    long len;
X
X    if (n >= BASEB) {
X        len = n / BASEB;
X        h = z.v;
X        lim = z.v + z.len - len;
X        while (h < lim) {
X            h[0] = h[len];
X            ++h;
X        }
X        n -= BASEB * len;
X        lim = z.v + z.len;
X        while (h < lim)
X            *h++ = 0;
X    }
X    if (n) {
X        len = z.len;
X        h = z.v + len - 1;
X        mask = 0;
X        while (len--) {
X            maskt = (((FULL) *h) << (BASEB - n)) & BASE1;
X            *h = (*h >> n) | mask;
X            mask = maskt;
X            --h;
X        }
X    }
X}
X
X
X/*
X * Utility routine to shift left.
X */
Xvoid
Xzshiftl(z, n)
X    ZVALUE z;
X    long n;
X{
X    register HALF *h;
X    FULL mask, i;
X    long len;
X
X    if (n >= BASEB) {
X        len = n / BASEB;
X        h = z.v + z.len - 1;
X        while (!*h)
X            --h;
X        while (h >= z.v) {
X            h[len] = h[0];
X            --h;
X        }
X        n -= BASEB * len;
X        while (len)
X            h[len--] = 0;
X    }
X    if (n > 0) {
X        len = z.len;
X        h = z.v;
X        mask = 0;
X        while (len--) {
X            i = (((FULL) *h) << n) | mask;
X            if (i > BASE1) {
X                mask = i >> BASEB;
X                i &= BASE1;
X            } else
X                mask = 0;
X            *h = (HALF) i;
X            ++h;
X        }
X    }
X}
X
X/*
X * initmasks - init the bitmask rotation arrays
X *
X * bitmask[i]      (1 << (i-1)),  for  -BASEB*4<=i<=BASEB*4
X * rotmask[j][k] (1 << ((j+k-1)%BASEB)),  for  -BASEB*2<=j,k<=BASEB*2
X *
X * The bmask array contains 8 cycles of rotations of a bit mask.
X * We point bitmask and the rotmask pointers into the middle to
X * ensure that we can have a complete two cycle swing forward
X * and backward.
X */
Xvoid
Xinitmasks()
X{
X    int i;
X
X    /*
X     * setup the bmask array
X     */
X    bmask = alloc((long)((8*BASEB)+1));
X    for (i=0; i < (8*BASEB)+1; ++i) {
X        bmask[i] = 1 << (i%BASEB);
X    }
X
X    /*
X     * setup the rmask pointers
X     */
X    rmask = (HALF **)malloc(sizeof(HALF *)*((BASEB*4)+2));
X    for (i = 0; i <= (4*BASEB)+1; ++i) {
X        rmask[i] = &bmask[(2*BASEB)+i];
X    }
X
X#if 0
X    /*
X     * setup the rotmask array, allow for -2*BASEB thru 2*BASEB indexing
X     */
X    rotmask = &rmask[2*BASEB];
X#endif
X
X    /*
X     * setup the bitmask array to allow -4*BASEB thru 4*BASEB indexing
X     */
X    bitmask = &bmask[4*BASEB];
X    return;
X}
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/zmath.c || echo "restore of calc2.9.0/zmath.c fails"
set `wc -c calc2.9.0/zmath.c`;Sum=$1
if test "$Sum" != "32248"
then echo original size 32248, current size $Sum;fi
echo "x - extracting calc2.9.0/zmath.h (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/zmath.h &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Data structure declarations for extended precision integer arithmetic.
X * The assumption made is that a long is 32 bits and shorts are 16 bits,
X * and longs must be addressible on word boundaries.
X */
X
X#ifndef    ZMATH_H
X#define    ZMATH_H
X
X#include <stdio.h>
X#include "alloc.h"
X#include "endian.h"
X
X#include "have_stdlib.h"
X#ifdef HAVE_STDLIB_H
X# include <stdlib.h>
X#endif
X
X
X#ifndef ALLOCTEST
X# if defined(CALC_MALLOC)
X#  define freeh(p) (((void *)p == (void *)_zeroval_) ||            \
X            ((void *)p == (void *)_oneval_) || free((void *)p))
X# else
X#  define freeh(p) { if (((void *)p != (void *)_zeroval_) &&        \
X             ((void *)p != (void *)_oneval_)) free((void *)p); }
X# endif
X#endif
X
X
Xtypedef    short FLAG;            /* small value (e.g. comparison) */
Xtypedef int BOOL;            /* TRUE or FALSE value */
Xtypedef unsigned long HASH;        /* hash value */
X
X
X#if !defined(TRUE)
X#define    TRUE    ((BOOL) 1)            /* booleans */
X#endif
X#if !defined(FALSE)
X#define    FALSE    ((BOOL) 0)
X#endif
X
X
X/*
X * NOTE: FULL must be twice the storage size of a HALF
X *     LEN storage size must be <= FULL storage size
X */
Xtypedef unsigned short HALF;        /* unit of number storage */
Xtypedef unsigned long FULL;        /* double unit of number storage */
Xtypedef long LEN;            /* unit of length storage */
X
X#define BASE    ((FULL) 65536)        /* base for calculations (2^16) */
X#define BASE1    ((FULL) (BASE - 1))    /* one less than base */
X#define BASEB    16            /* number of bits in base */
X#define    BASEDIG    5            /* number of digits in base */
X#define    MAXHALF    ((FULL) 0x7fff)        /* largest positive half value */
X#define    MAXFULL    ((FULL) 0x7fffffff)    /* largest positive full value */
X#define    TOPHALF    ((FULL) 0x8000)        /* highest bit in half value */
X#define    TOPFULL    ((FULL) 0x80000000)    /* highest bit in full value */
X#define MAXLEN    ((LEN)    0x7fffffff)    /* longest value allowed */
X#define    MAXREDC    5            /* number of entries in REDC cache */
X#define    SQ_ALG2    20            /* size for alternative squaring */
X#define    MUL_ALG2 20            /* size for alternative multiply */
X#define    POW_ALG2 40            /* size for using REDC for powers */
X#define    REDC_ALG2 50            /* size for using alternative REDC */
X
X
X
Xtypedef union {
X    FULL    ivalue;
X    struct {
X        HALF Svalue1;
X        HALF Svalue2;
X    } sis;
X} SIUNION;
X
X
X#if !defined(BYTE_ORDER)
X#include <machine/endian.h>
X#endif
X
X#if !defined(LITTLE_ENDIAN)
X#define LITTLE_ENDIAN    1234    /* Least Significant Byte first */
X#endif
X#if !defined(BIG_ENDIAN)
X#define BIG_ENDIAN    4321    /* Most Significant Byte first */
X#endif
X/* PDP_ENDIAN - LSB in word, MSW in long is not supported */
X
X#if BYTE_ORDER == LITTLE_ENDIAN
X# define silow    sis.Svalue1    /* low order half of full value */
X# define sihigh    sis.Svalue2    /* high order half of full value */
X#else
X# if BYTE_ORDER == BIG_ENDIAN
X#  define silow    sis.Svalue2    /* low order half of full value */
X#  define sihigh sis.Svalue1    /* high order half of full value */
X# else
X   :@</*/>@:    BYTE_ORDER must be BIG_ENDIAN or LITTLE_ENDIAN    :@</*/>@:
X# endif
X#endif
X
X
Xtypedef struct {
X    HALF    *v;        /* pointer to array of values */
X    LEN    len;        /* number of values in array */
X    BOOL    sign;        /* sign, nonzero is negative */
X} ZVALUE;
X
X
X
X/*
X * Function prototypes for integer math routines.
X */
X#if defined(__STDC__)
X#define MATH_PROTO(a) a
X#else
X#define MATH_PROTO(a) ()
X#endif
X
Xextern HALF * alloc MATH_PROTO((LEN len));
X#ifdef    ALLOCTEST
Xextern void freeh MATH_PROTO((HALF *));
X#endif
X
X
X/*
X * Input, output, and conversion routines.
X */
Xextern void zcopy MATH_PROTO((ZVALUE z, ZVALUE *res));
Xextern void itoz MATH_PROTO((long i, ZVALUE *res));
Xextern void atoz MATH_PROTO((char *s, ZVALUE *res));
Xextern long ztoi MATH_PROTO((ZVALUE z));
Xextern void zprintval MATH_PROTO((ZVALUE z, long decimals, long width));
Xextern void zprintx MATH_PROTO((ZVALUE z, long width));
Xextern void zprintb MATH_PROTO((ZVALUE z, long width));
Xextern void zprinto MATH_PROTO((ZVALUE z, long width));
X
X
X/*
X * Basic numeric routines.
X */
Xextern void zmuli MATH_PROTO((ZVALUE z, long n, ZVALUE *res));
Xextern long zdivi MATH_PROTO((ZVALUE z, long n, ZVALUE *res));
Xextern long zmodi MATH_PROTO((ZVALUE z, long n));
Xextern void zadd MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zsub MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zmul MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zdiv MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res, ZVALUE *rem));
Xextern void zquo MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *rem));
Xextern BOOL zdivides MATH_PROTO((ZVALUE z1, ZVALUE z2));
Xextern void zor MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zand MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zxor MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zshift MATH_PROTO((ZVALUE z, long n, ZVALUE *res));
Xextern void zsquare MATH_PROTO((ZVALUE z, ZVALUE *res));
Xextern long zlowbit MATH_PROTO((ZVALUE z));
Xextern long zhighbit MATH_PROTO((ZVALUE z));
Xextern void zbitvalue MATH_PROTO((long n, ZVALUE *res));
Xextern BOOL zisset MATH_PROTO((ZVALUE z, long n));
Xextern BOOL zisonebit MATH_PROTO((ZVALUE z));
Xextern BOOL zisallbits MATH_PROTO((ZVALUE z));
Xextern FLAG ztest MATH_PROTO((ZVALUE z));
Xextern FLAG zrel MATH_PROTO((ZVALUE z1, ZVALUE z2));
Xextern BOOL zcmp MATH_PROTO((ZVALUE z1, ZVALUE z2));
X
X
X/*
X * More complicated numeric functions.
X */
Xextern long iigcd MATH_PROTO((long i1, long i2));
Xextern void zgcd MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zlcm MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zreduce MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *z1res, ZVALUE *z2res));
Xextern void zfact MATH_PROTO((ZVALUE z, ZVALUE *dest));
Xextern void zpfact MATH_PROTO((ZVALUE z, ZVALUE *dest));
Xextern void zlcmfact MATH_PROTO((ZVALUE z, ZVALUE *dest));
Xextern void zperm MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zcomb MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern BOOL zprimetest MATH_PROTO((ZVALUE z, long count));
Xextern FLAG zjacobi MATH_PROTO((ZVALUE z1, ZVALUE z2));
Xextern void zfib MATH_PROTO((ZVALUE z, ZVALUE *res));
Xextern void zpowi MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void ztenpow MATH_PROTO((long power, ZVALUE *res));
Xextern void zpowermod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
Xextern BOOL zmodinv MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern BOOL zrelprime MATH_PROTO((ZVALUE z1, ZVALUE z2));
Xextern long zlog MATH_PROTO((ZVALUE z1, ZVALUE z2));
Xextern long zlog10 MATH_PROTO((ZVALUE z));
Xextern long zdivcount MATH_PROTO((ZVALUE z1, ZVALUE z2));
Xextern long zfacrem MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *rem));
Xextern void zgcdrem MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern long zlowfactor MATH_PROTO((ZVALUE z, long count));
Xextern long zdigits MATH_PROTO((ZVALUE z1));
Xextern FLAG zdigit MATH_PROTO((ZVALUE z1, long n));
Xextern BOOL zsqrt MATH_PROTO((ZVALUE z1, ZVALUE *dest));
Xextern void zroot MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *dest));
Xextern BOOL zissquare MATH_PROTO((ZVALUE z));
Xextern HASH zhash MATH_PROTO((ZVALUE z));
X
X#if 0
Xextern void zapprox MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res1, ZVALUE *res2));
X#endif
X
X
X#if 0
Xextern void zmulmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
Xextern void zsquaremod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));
Xextern void zsubmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
SHAR_EOF
echo "End of part 12"
echo "File calc2.9.0/zmath.h is continued in part 13"
echo "13" > s2_seq_.tmp
exit 0