Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/unix/volume26/calc/part10

Text File | 1992-05-09 | 53.6 KB | 2,519 lines

Newsgroups: comp.sources.unix
From: [e-mail] (David I. Bell)
Subject: v26i036: CALC - An arbitrary precision C-like calculator, Part10/21
Sender: [e-mail]
Approved: [e-mail]

Submitted-By: [e-mail] (David I. Bell)
Posting-Number: Volume 26, Issue 36
Archive-Name: calc/part10

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 10 (of 21)."
# Contents:  io.c qfunc.c
# Wrapped by dbell@elm on Tue Feb 25 15:21:08 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'io.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'io.c'\"
else
echo shar: Extracting \"'io.c'\" \(26180 characters\)
sed "s/^X//" >'io.c' <<'END_OF_FILE'
X/*
X * Copyright (c) 1992 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 extended precision numbers
X */
X
X#include <stdio.h>
X#include "stdarg.h"
X#include "math.h"
X
X#define    OUTBUFSIZE    200        /* realloc size for output buffers */
X
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 long    scalefactor;
Xstatic long    outbufsize;
Xstatic long    outbufused;
Xstatic BOOL    outputisstring;
Xstatic ZVALUE    scalenumber = { 0, 0, 0 };
X
X#if 0
Xstatic long    etoalen;
Xstatic char    *etoabuf = NULL;
X#endif
X
Xstatic void zprintx proto((ZVALUE z, long width));
Xstatic void zprintb proto((ZVALUE z, long width));
Xstatic void zprinto proto((ZVALUE z, long width));
Xstatic void zprintval proto((ZVALUE z, long decimals, long width));
X
Xstatic void atoz proto((char * s, ZVALUE * res));
X
Xstatic void qprintff();
Xstatic void qprintfd();
Xstatic void qprintfe();
Xstatic void qprintfr();
Xstatic void qprintfo();
Xstatic void qprintfb();
X#ifdef CODE
Xextern void qprintfx();
X#else
Xstatic void qprintfx();
X#endif
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            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            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 * 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
Xdivertio()
X{
X    register IOSTATE *sp;
X
X    sp = (IOSTATE *) malloc(sizeof(IOSTATE));
X    if (sp == NULL)
X        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        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 *
Xgetdivertedio()
X{
X    register IOSTATE *sp;
X    char *cp;
X
X    sp = oldiostates;
X    if (sp == NULL)
X        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
Xcleardiversions()
X{
X    while (oldiostates)
X        free(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
Xsetfp(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 */
Xvoid
Xsetmode(newmode)
X{
X    if ((newmode <= MODE_DEFAULT) || (newmode > MODE_MAX))
X        error("Setting illegal output mode");
X    _outmode_ = newmode;
X}
X
X
X/*
X * Set the number of digits for float or exponential output.
X */
Xvoid
Xsetdigits(newdigits)
X    long newdigits;
X{
X    if (newdigits < 0)
X        error("Setting illegal number of digits");
X    _outdigits_ = newdigits;
X}
X
X
X/*
X * Print a formatted string containing arbitrary numbers, similar to printf.
X * ALL numeric arguments to this routine are rational NUMBERs.
X * Various forms of printing such numbers are supplied, in addition
X * to strings and characters.  Output can actually be to any FILE
X * stream or a string.
X */
X#ifdef VARARGS
X# define VA_ALIST1 fmt, va_alist
X# define VA_DCL1 char *fmt; va_dcl
X#else
X# ifdef __STDC__
X#  define VA_ALIST1 char *fmt, ...
X#  define VA_DCL1
X# else
X#  define VA_ALIST1 fmt
X#  define VA_DCL1 char *fmt;
X# endif
X#endif
X/*VARARGS*/
Xvoid
Xqprintf(VA_ALIST1)
X    VA_DCL1
X{
X    va_list ap;
X    NUMBER *q;
X    int ch, sign;
X    long width, precision;
X
X#ifdef VARARGS
X    va_start(ap);
X#else
X    va_start(ap, fmt);
X#endif
X    while ((ch = *fmt++) != '\0') {
X        if (ch == '\\') {
X            ch = *fmt++;
X            switch (ch) {
X                case 'n': ch = '\n'; break;
X                case 'r': ch = '\r'; break;
X                case 't': ch = '\t'; break;
X                case 'f': ch = '\f'; break;
X                case 'v': ch = '\v'; break;
X                case 'b': ch = '\b'; break;
X                case 0:
X                    va_end(ap);
X                    return;
X            }
X            PUTCHAR(ch);
X            continue;
X        }
X        if (ch != '%') {
X            PUTCHAR(ch);
X            continue;
X        }
X        ch = *fmt++;
X        width = 0; precision = 8; sign = 1;
Xpercent:    ;
X        switch (ch) {
X            case 'd':
X                q = va_arg(ap, NUMBER *);
X                qprintfd(q, width);
X                break;
X            case 'f':
X                q = va_arg(ap, NUMBER *);
X                qprintff(q, width, precision);
X                break;
X            case 'e':
X                q = va_arg(ap, NUMBER *);
X                qprintfe(q, width, precision);
X                break;
X            case 'r':
X            case 'R':
X                q = va_arg(ap, NUMBER *);
X                qprintfr(q, width, (BOOL) (ch == 'R'));
X                break;
X            case 'N':
X                q = va_arg(ap, NUMBER *);
X                zprintval(q->num, 0L, width);
X                break;
X            case 'D':
X                q = va_arg(ap, NUMBER *);
X                zprintval(q->den, 0L, width);
X                break;
X            case 'o':
X                q = va_arg(ap, NUMBER *);
X                qprintfo(q, width);
X                break;
X            case 'x':
X                q = va_arg(ap, NUMBER *);
X                qprintfx(q, width);
X                break;
X            case 'b':
X                q = va_arg(ap, NUMBER *);
X                qprintfb(q, width);
X                break;
X            case 's':
X                PUTSTR(va_arg(ap, char *));
X                break;
X            case 'c':
X                PUTCHAR(va_arg(ap, int));
X                break;
X            case 0:
X                va_end(ap);
X                return;
X            case '-':
X                sign = -1;
X                ch = *fmt++;
X            default:
X        if (('0' <= ch && ch <= '9') || ch == '.' || ch == '*') {
X            if (ch == '*') {
X                q = va_arg(ap, NUMBER *);
X                width = sign * qtoi(q);
X                ch = *fmt++;
X            } else if (ch != '.') {
X                width = ch - '0';
X                while ('0' <= (ch = *fmt++) && ch <= '9')
X                    width = width * 10 + ch - '0';
X                width *= sign;
X            }
X            if (ch == '.') {
X                if ((ch = *fmt++) == '*') {
X                    q = va_arg(ap, NUMBER *);
X                    precision = qtoi(q);
X                    ch = *fmt++;
X                } else {
X                    precision = 0;
X                    while ('0' <= (ch = *fmt++) && ch <= '9')
X                        precision = precision * 10 + ch - '0';
X                }
X            }
X            goto percent;
X        }
X        }
X    }
X    va_end(ap);
X}
X
X
X#if 0
X/*
X * Read a number from the specified FILE stream (NULL means stdin).
X * The number can be an integer, a fraction, a real number, an
X * exponential number, or a hex, octal or binary number.  Leading blanks
X * are skipped.  Illegal numbers return NULL.  Unrecognized characters
X * remain to be read on the line.
X *    q = qreadval(fp);
X */
XNUMBER *
Xqreadval(fp)
X    FILE *fp;        /* file stream to read from (or NULL) */
X{
X    NUMBER *r;        /* returned number */
X    char *cp;         /* current buffer location */
X    long savecc;        /* characters saved in buffer */
X    long scancc;         /* characters parsed correctly */
X    int ch;            /* current character */
X
X    if (fp == NULL)
X        fp = stdin;
X    if (etoabuf == NULL) {
X        etoabuf = (char *)malloc(OUTBUFSIZE + 2);
X        if (etoabuf == NULL)
X            return NULL;
X        etoalen = OUTBUFSIZE;
X    }
X    cp = etoabuf;
X    ch = fgetc(fp);
X    while ((ch == ' ') || (ch == '\t'))
X        ch = fgetc(fp);
X    savecc = 0;
X    for (;;) {
X        if (ch == EOF)
X            return NULL;
X        if (savecc >= etoalen)
X        {
X            cp = (char *)realloc(etoabuf, etoalen + OUTBUFSIZE + 2);
X            if (cp == NULL)
X                return NULL;
X            etoabuf = cp;
X            etoalen += OUTBUFSIZE;
X            cp += savecc;
X        }
X        *cp++ = (char)ch;
X        *cp = '\0';
X        scancc = qparse(etoabuf, QPF_SLASH);
X        if (scancc != ++savecc)
X            break;
X        ch = fgetc(fp);
X    }
X    ungetc(ch, fp);
X    if (scancc < 0)
X        return NULL;
X    r = atoq(etoabuf);
X    if (iszero(r->den)) {
X        qfree(r);
X        r = NULL;
X    }
X    return r;
X}
X#endif
X
X
X/*
X * Print a complex number in rational representation.
X * Example:  2/3-4i/5
X */
Xvoid
Xcprintfr(c)
X    COMPLEX *c;
X{
X    NUMBER *r;
X    NUMBER *i;
X
X    r = c->real;
X    i = c->imag;
X    if (!qiszero(r) || qiszero(i))
X        qprintfr(r, 0L, FALSE);
X    if (qiszero(i))
X        return;
X    if (!qiszero(r) && !qisneg(i))
X        PUTCHAR('+');
X    zprintval(i->num, 0L, 0L);
X    PUTCHAR('i');
X    if (qisfrac(i)) {
X        PUTCHAR('/');
X        zprintval(i->den, 0L, 0L);
X    }
X}
X
X
X/*
X * Print a number in the specified output mode.
X * If MODE_DEFAULT is given, then the default output mode is used.
X * Any approximate output is flagged with a leading tilde.
X * Integers are always printed as themselves.
X */
Xvoid
Xqprintnum(q, outmode)
X    NUMBER *q;
X{
X    NUMBER tmpval;
X    long prec, exp;
X
X    if (outmode == MODE_DEFAULT)
X        outmode = _outmode_;
X    if ((outmode == MODE_FRAC) || ((outmode == MODE_REAL) && qisint(q))) {
X        qprintfr(q, 0L, FALSE);
X        return;
X    }
X    switch (outmode) {
X        case MODE_INT:
X            if (qisfrac(q))
X                PUTCHAR('~');
X            qprintfd(q, 0L);
X            break;
X
X        case MODE_REAL:
X            prec = qplaces(q);
X            if ((prec < 0) || (prec > _outdigits_)) {
X                prec = _outdigits_;
X                PUTCHAR('~');
X            }
X            qprintff(q, 0L, prec);
X            break;
X
X        case MODE_EXP:
X            if (qiszero(q)) {
X                PUTCHAR('0');
X                return;
X            }
X            tmpval = *q;
X            tmpval.num.sign = 0;
X            exp = qilog10(&tmpval);
X            if (exp == 0) {        /* in range to output as real */
X                qprintnum(q, MODE_REAL);
X                return;
X            }
X            tmpval.num = _one_;
X            tmpval.den = _one_;
X            if (exp > 0)
X                ztenpow(exp, &tmpval.den);
X            else
X                ztenpow(-exp, &tmpval.num);
X            q = qmul(q, &tmpval);
X            freeh(tmpval.num.v);
X            freeh(tmpval.den.v);
X            qprintnum(q, MODE_REAL);
X            qfree(q);
X            PRINTF1("e%ld", exp);
X            break;
X
X        case MODE_HEX:
X            qprintfx(q, 0L);
X            break;
X
X        case MODE_OCTAL:
X            qprintfo(q, 0L);
X            break;
X
X        case MODE_BINARY:
X            qprintfb(q, 0L);
X            break;
X
X        default:
X            error("Bad mode for print");
X    }
X}
X
X
X/*
X * Print a number in floating point representation.
X * Example:  193.784
X */
Xstatic void
Xqprintff(q, width, precision)
X    NUMBER *q;
X    long width;
X    long precision;
X{
X    ZVALUE z, z1;
X
X    if (precision != scalefactor) {
X        if (scalenumber.v)
X            freeh(scalenumber.v);
X        ztenpow(precision, &scalenumber);
X        scalefactor = precision;
X    }
X    if (scalenumber.v)
X        zmul(q->num, scalenumber, &z);
X    else
X        z = q->num;
X    if (qisfrac(q)) {
X        zquo(z, q->den, &z1);
X        if (z.v != q->num.v)
X            freeh(z.v);
X        z = z1;
X    }
X    if (qisneg(q) && iszero(z))
X        PUTCHAR('-');
X    zprintval(z, precision, width);
X    if (z.v != q->num.v)
X        freeh(z.v);
X}
X
X
X/*
X * Print a number in exponential notation.
X * Example: 4.1856e34
X */
X/*ARGSUSED*/
Xstatic void
Xqprintfe(q, width, precision)
X    register NUMBER *q;
X    long width;
X    long precision;
X{
X    long exponent;
X    NUMBER q2;
X    ZVALUE num, den, tenpow, tmp;
X
X    if (qiszero(q)) {
X        PUTSTR("0.0");
X        return;
X    }
X    num = q->num;
X    den = q->den;
X    num.sign = 0;
X    exponent = zdigits(num) - zdigits(den);
X    if (exponent > 0) {
X        ztenpow(exponent, &tenpow);
X        zmul(den, tenpow, &tmp);
X        freeh(tenpow.v);
X        den = tmp;
X    }
X    if (exponent < 0) {
X        ztenpow(-exponent, &tenpow);
X        zmul(num, tenpow, &tmp);
X        freeh(tenpow.v);
X        num = tmp;
X    }
X    if (zrel(num, den) < 0) {
X        zmuli(num, 10L, &tmp);
X        if (num.v != q->num.v)
X            freeh(num.v);
X        num = tmp;
X        exponent--;
X    }
X    q2.num = num;
X    q2.den = den;
X    q2.num.sign = q->num.sign;
X    qprintff(&q2, 0L, precision);
X    if (exponent)
X        PRINTF1("e%ld", exponent);
X    if (num.v != q->num.v)
X        freeh(num.v);
X    if (den.v != q->den.v)
X        freeh(den.v);
X}
X
X
X/*
X * Print a number in rational representation.
X * Example: 397/37
X */
Xstatic void
Xqprintfr(q, width, force)
X    NUMBER *q;
X    long width;
X    BOOL force;
X{
X    zprintval(q->num, 0L, width);
X    if (force || qisfrac(q)) {
X        PUTCHAR('/');
X        zprintval(q->den, 0L, width);
X    }
X}
X
X
X/*
X * Print a number as an integer (truncating fractional part).
X * Example: 958421
X */
Xstatic void
Xqprintfd(q, width)
X    NUMBER *q;
X    long width;
X{
X    ZVALUE z;
X
X    if (qisfrac(q)) {
X        zquo(q->num, q->den, &z);
X        zprintval(z, 0L, width);
X        freeh(z.v);
X    } else
X        zprintval(q->num, 0L, width);
X}
X
X
X/*
X * Print a number in hex.
X * This prints the numerator and denominator in hex.
X */
X#ifndef CODE
Xstatic
X#endif
Xvoid
Xqprintfx(q, width)
X    NUMBER *q;
X    long width;
X{
X    zprintx(q->num, width);
X    if (qisfrac(q)) {
X        PUTCHAR('/');
X        zprintx(q->den, 0L);
X    }
X}
X
X
X/*
X * Print a number in binary.
X * This prints the numerator and denominator in binary.
X */
Xstatic void
Xqprintfb(q, width)
X    NUMBER *q;
X    long width;
X{
X    zprintb(q->num, width);
X    if (qisfrac(q)) {
X        PUTCHAR('/');
X        zprintb(q->den, 0L);
X    }
X}
X
X
X/*
X * Print a number in octal.
X * This prints the numerator and denominator in octal.
X */
Xstatic void
Xqprintfo(q, width)
X    NUMBER *q;
X    long width;
X{
X    zprinto(q->num, width);
X    if (qisfrac(q)) {
X        PUTCHAR('/');
X        zprinto(q->den, 0L);
X    }
X}
X
X
X/*
X * Print an integer value as a hex number.
X * The special characters 0x appear to indicate the number is hex.
X */
X/*ARGSUSED*/
Xstatic void
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
X    len = z.len - 1;
X    if (isneg(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*/
Xstatic void
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
X    len = z.len - 1;
X    if (isneg(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*/
Xstatic void
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
X    if (isneg(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 */
Xstatic void
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 (!iszero(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                freeh(leftnums[n].v);
X            n--;
X        }
X        freeh(leftnums[n].v);
X        leftnums[n] = rightnums[n];
X        rightnums[n].len = 0;
X    }
X}
X
X
X/*
X * Convert a string to a number in rational, floating point,
X * exponential notation, hex, or octal.
X *    q = atoq(string);
X */
XNUMBER *
Xatoq(s)
X    register char *s;
X{
X    register NUMBER *q;
X    register char *t;
X    ZVALUE div, newnum, newden, tmp;
X    long decimals, exp;
X    BOOL hex, negexp;
X
X    q = qalloc();
X    decimals = 0;
X    exp = 0;
X    negexp = FALSE;
X    hex = FALSE;
X    t = s;
X    if ((*t == '+') || (*t == '-'))
X        t++;
X    if ((*t == '0') && ((t[1] == 'x') || (t[1] == 'X'))) {
X        hex = TRUE;
X        t += 2;
X    }
X    while (((*t >= '0') && (*t <= '9')) || (hex &&
X        (((*t >= 'a') && (*t <= 'f')) || ((*t >= 'A') && (*t <= 'F')))))
X            t++;
X    if (*t == '/') {
X        t++;
X        atoz(t, &q->den);
X    } else if ((*t == '.') || (*t == 'e') || (*t == 'E')) {
X        if (*t == '.') {
X            t++;
X            while ((*t >= '0') && (*t <= '9')) {
X                t++;
X                decimals++;
X            }
X        }
X        /*
X         * Parse exponent if any
X         */
X        if ((*t == 'e') || (*t == 'E')) {
X            t++;
X            if (*t == '+')
X                t++;
X            else if (*t == '-') {
X                negexp = TRUE;
X                t++;
X            }
X            while ((*t >= '0') && (*t <= '9')) {
X                exp = (exp * 10) + *t++ - '0';
X                if (exp > 1000000)
X                    error("Exponent too large");
X            }
X        }
X        ztenpow(decimals, &q->den);
X    }
X    atoz(s, &q->num);
X    if (qiszero(q)) {
X        qfree(q);
X        return qlink(&_qzero_);
X    }
X    /*
X     * Apply the exponential if any
X     */
X    if (exp) {
X        ztenpow(exp, &tmp);
X        if (negexp) {
X            zmul(q->den, tmp, &newden);
X            freeh(q->den.v);
X            q->den = newden;
X        } else {
X            zmul(q->num, tmp, &newnum);
X            freeh(q->num.v);
X            q->num = newnum;
X        }
X        freeh(tmp.v);
X    }
X    /*
X     * Reduce the fraction to lowest terms
X     */
X    if (isunit(q->num) || isunit(q->den))
X        return q;
X    zgcd(q->num, q->den, &div);
X    if (isunit(div))
X        return q;
X    zquo(q->num, div, &newnum);
X    freeh(q->num.v);
X    zquo(q->den, div, &newden);
X    freeh(q->den.v);
X    q->num = newnum;
X    q->den = newden;
X    return q;
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 */
Xstatic void
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        freeh(z.v);
X        zadd(ztmp, digit, &z);
X        freeh(ztmp.v);
X    }
X    trim(&z);
X    if (minus && !iszero(z))
X        z.sign = 1;
X    *res = z;
X}
X
X
X/*
X * Parse a number in any of the various legal forms, and return the count
X * of characters that are part of a legal number.  Numbers can be either a
X * decimal integer, possibly two decimal integers separated with a slash, a
X * floating point or exponential number, a hex number beginning with "0x",
X * a binary number beginning with "0b", or an octal number beginning with "0".
X * The flags argument modifies the end of number testing for ease in handling
X * fractions or complex numbers.  Minus one is returned if the number format
X * is definitely illegal.
X */
Xlong
Xqparse(cp, flags)
X    register char *cp;
X{
X    char *oldcp;
X
X    oldcp = cp;
X    if ((*cp == '+') || (*cp == '-'))
X        cp++;
X    if ((*cp == '+') || (*cp == '-'))
X        return -1;
X    if ((*cp == '0') && ((cp[1] == 'x') || (cp[1] == 'X'))) {    /* hex */
X        cp += 2;
X        while (((*cp >= '0') && (*cp <= '9')) ||
X            ((*cp >= 'a') && (*cp <= 'f')) ||
X            ((*cp >= 'A') && (*cp <= 'F')))
X                cp++;
X        goto done;
X    }
X    if ((*cp == '0') && ((cp[1] == 'b') || (cp[1] == 'B'))) {    /* binary */
X        cp += 2;
X        while ((*cp == '0') || (*cp == '1'))
X            cp++;
X        goto done;
X    }
X    if ((*cp == '0') && (cp[1] >= '0') && (cp[1] <= '9')) { /* octal */
X        while ((*cp >= '0') && (*cp <= '7'))
X            cp++;
X        goto done;
X    }
X    /*
X     * Number is decimal, but can still be a fraction or real or exponential.
X     */
X    while ((*cp >= '0') && (*cp <= '9'))
X        cp++;
X    if (*cp == '/' && flags & QPF_SLASH) {    /* fraction */
X        cp++;
X        while ((*cp >= '0') && (*cp <= '9'))
X            cp++;
X        goto done;
X    }
X    if (*cp == '.') {    /* floating point */
X        cp++;
X        while ((*cp >= '0') && (*cp <= '9'))
X            cp++;
X    }
X    if ((*cp == 'e') || (*cp == 'E')) {    /* exponential */
X        cp++;
X        if ((*cp == '+') || (*cp == '-'))
X            cp++;
X        if ((*cp == '+') || (*cp == '-'))
X            return -1;
X        while ((*cp >= '0') && (*cp <= '9'))
X            cp++;
X    }
X
Xdone:
X    if (((*cp == 'i') || (*cp == 'I')) && (flags & QPF_IMAG))
X        cp++;
X    if ((*cp == '.') || ((*cp == '/') && (flags & QPF_SLASH)) ||
X        ((*cp >= '0') && (*cp <= '9')) ||
X        ((*cp >= 'a') && (*cp <= 'z')) ||
X        ((*cp >= 'A') && (*cp <= 'Z')))
X            return -1;
X    return (cp - oldcp);
X}
X
X/* END CODE */
END_OF_FILE
if test 26180 -ne `wc -c <'io.c'`; then
    echo shar: \"'io.c'\" unpacked with wrong size!
fi
# end of 'io.c'
fi
if test -f 'qfunc.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'qfunc.c'\"
else
echo shar: Extracting \"'qfunc.c'\" \(24256 characters\)
sed "s/^X//" >'qfunc.c' <<'END_OF_FILE'
X/*
X * Copyright (c) 1992 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 rational arithmetic non-primitive functions
X */
X
X#include "math.h"
X
X
XNUMBER *_epsilon_;    /* default precision for real functions */
Xlong _epsilonprec_;    /* bits of precision for epsilon */
X
X#if 0
Xstatic char *abortmsg = "Calculation aborted";
Xstatic char *memmsg = "Not enough memory";
X#endif
X
X
X/*
X * Set the default precision for real calculations.
X * The precision must be between zero and one.
X */
Xvoid
Xsetepsilon(q)
X    NUMBER *q;        /* number to be set as the new epsilon */
X{
X    NUMBER *old;
X
X    if (qisneg(q) || qiszero(q) || (qreli(q, 1L) >= 0))
X        error("Epsilon value must be between zero and one");
X    old = _epsilon_;
X    _epsilonprec_ = qprecision(q);
X    _epsilon_ = qlink(q);
X    if (old)
X        qfree(old);
X}
X
X
X/*
X * Return the inverse of one number modulo another.
X * That is, find x such that:
X *    Ax = 1 (mod B)
X * Returns zero if the numbers are not relatively prime (temporary hack).
X */
XNUMBER *
Xqminv(q1, q2)
X    NUMBER *q1, *q2;
X{
X    NUMBER *r;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for minv");
X    r = qalloc();
X    if (zmodinv(q1->num, q2->num, &r->num)) {
X        qfree(r);
X        return qlink(&_qzero_);
X    }
X    return r;
X}
X
X
X/*
X * Return the modulo of one number raised to another.
X * Here q1 is the number to be raised, q2 is the power to raise it to,
X * and q3 is the number to take the modulo with the result.
X * The second and third numbers are assumed nonnegative.
X * Returned answer is non-negative.
X *    q4 = qpowermod(q1, q2, q3);
X */
XNUMBER *
Xqpowermod(q1, q2, q3)
X    NUMBER *q1, *q2, *q3;
X{
X    NUMBER *r;
X
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        error("Non-integers for powermod");
X    r = qalloc();
X    zpowermod(q1->num, q2->num, q3->num, &r->num);
X    return r;
X}
X
X
X/*
X * Return the power function of one number with another.
X * The power must be integral.
X *    q3 = qpowi(q1, q2);
X */
XNUMBER *
Xqpowi(q1, q2)
X    NUMBER *q1, *q2;
X{
X    register NUMBER *r;
X    BOOL invert, sign;
X    ZVALUE num, den, z2;
X
X    if (qisfrac(q2))
X        error("Raising number to fractional power");
X    num = q1->num;
X    den = q1->den;
X    z2 = q2->num;
X    sign = (num.sign && isodd(z2));
X    invert = z2.sign;
X    num.sign = 0;
X    z2.sign = 0;
X    /*
X    * Check for trivial cases first.
X    */
X    if (iszero(num)) {    /* zero raised to a power */
X        if (invert || iszero(z2))
X            error("Zero raised to non-positive power");
X        return qlink(&_qzero_);
X    }
X    if (isunit(num) && isunit(den)) {    /* 1 or -1 raised to a power */
X        r = (sign ? q1 : &_qone_);
X        r->links++;
X        return r;
X    }
X    if (iszero(z2))    /* raising to zeroth power */
X        return qlink(&_qone_);
X    if (isunit(z2)) {    /* raising to power 1 or -1 */
X        if (invert)
X            return qinv(q1);
X        return qlink(q1);
X    }
X    /*
X     * Not a trivial case.  Do the real work.
X     */
X    r = qalloc();
X    if (!isunit(num))
X        zpowi(num, z2, &r->num);
X    if (!isunit(den))
X        zpowi(den, z2, &r->den);
X    if (invert) {
X        z2 = r->num;
X        r->num = r->den;
X        r->den = z2;
X    }
X    r->num.sign = sign;
X    return r;
X}
X
X
X/*
X * Given the legs of a right triangle, compute its hypothenuse within
X * the specified error.  This is sqrt(a^2 + b^2).
X */
XNUMBER *
Xqhypot(q1, q2, epsilon)
X    NUMBER *q1, *q2, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *tmp3;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Bad epsilon value for hypot");
X    if (qiszero(q1))
X        return qabs(q2);
X    if (qiszero(q2))
X        return qabs(q1);
X    tmp1 = qsquare(q1);
X    tmp2 = qsquare(q2);
X    tmp3 = qadd(tmp1, tmp2);
X    qfree(tmp1);
X    qfree(tmp2);
X    tmp1 = qsqrt(tmp3, epsilon);
X    qfree(tmp3);
X    return tmp1;
X}
X
X
X/*
X * Given one leg of a right triangle with unit hypothenuse, calculate
X * the other leg within the specified error.  This is sqrt(1 - a^2).
X * If the wantneg flag is nonzero, then negative square root is returned.
X */
XNUMBER *
Xqlegtoleg(q, epsilon, wantneg)
X    NUMBER *q, *epsilon;
X    BOOL wantneg;
X{
X    NUMBER *qt, *res, qtmp;
X    ZVALUE num, ztmp1, ztmp2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Bad epsilon value for legtoleg");
X    if (qisunit(q))
X        return qlink(&_qzero_);
X    if (qiszero(q)) {
X        if (wantneg)
X            return qlink(&_qnegone_);
X        return qlink(&_qone_);
X    }
X    num = q->num;
X    num.sign = 0;
X    if (zrel(num, q->den) >= 0)
X        error("Leg too large in legtoleg");
X    zsquare(q->den, &ztmp1);
X    zsquare(num, &ztmp2);
X    zsub(ztmp1, ztmp2, &qtmp.num);
X    freeh(ztmp1.v);
X    freeh(ztmp2.v);
X    qtmp.den = _one_;
X    qt = qsqrt(&qtmp, epsilon);
X    freeh(qtmp.num.v);
X    qtmp.num = q->den;
X    res = qdiv(qt, &qtmp);
X    qfree(qt);
X    qt = qbappr(res, epsilon);
X    qfree(res);
X    if (wantneg) {
X        res = qneg(qt);
X        qfree(qt);
X        qt = res;
X    }
X    return qt;
X}
X
X
X/*
X * Compute the square root of a number to within the specified error.
X * If the number is an exact square, the exact result is returned.
X *    q3 = qsqrt(q1, q2);
X */
XNUMBER *
Xqsqrt(q1, epsilon)
X    register NUMBER *q1, *epsilon;
X{
X    long bits, bits2;    /* number of bits of precision */
X    int exact;
X    NUMBER *r;
X    ZVALUE t1, t2;
X
X    if (qisneg(q1))
X        error("Square root of negative number");
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Bad epsilon value for sqrt");
X    if (qiszero(q1))
X        return qlink(&_qzero_);
X    if (qisunit(q1))
X        return qlink(&_qone_);
X    if (qiszero(epsilon) && qisint(q1) && istiny(q1->num) && (*q1->num.v <= 3))
X        return qlink(&_qone_);
X    bits = zhighbit(epsilon->den) - zhighbit(epsilon->num) + 1;
X    if (bits < 0)
X        bits = 0;
X    bits2 = zhighbit(q1->num) - zhighbit(q1->den) + 1;
X    if (bits2 > 0)
X        bits += bits2;
X    r = qalloc();
X    zshift(q1->num, bits * 2, &t2);
X    zmul(q1->den, t2, &t1);
X    freeh(t2.v);
X    exact = zsqrt(t1, &t2);
X    freeh(t1.v);
X    if (exact) {
X        zshift(q1->den, bits, &t1);
X        zreduce(t2, t1, &r->num, &r->den);
X    } else {
X        zquo(t2, q1->den, &t1);
X        freeh(t2.v);
X        zbitvalue(bits, &t2);
X        zreduce(t1, t2, &r->num, &r->den);
X    }
X    freeh(t1.v);
X    freeh(t2.v);
X    if (qiszero(r)) {
X        qfree(r);
X        r = qlink(&_qzero_);
X    }
X    return r;
X}
X
X
X/*
X * Calculate the integral part of the square root of a number.
X * Example:  qisqrt(13) = 3.
X */
XNUMBER *
Xqisqrt(q)
X    NUMBER *q;
X{
X    NUMBER *r;
X    ZVALUE tmp;
X
X    if (qisneg(q))
X        error("Square root of negative number");
X    if (qiszero(q))
X        return qlink(&_qzero_);
X    if (qisint(q) && istiny(q->num) && (z1tol(q->num) < 4))
X        return qlink(&_qone_);
X    r = qalloc();
X    if (qisint(q)) {
X        (void) zsqrt(q->num, &r->num);
X        return r;
X    }
X    zquo(q->num, q->den, &tmp);
X    (void) zsqrt(tmp, &r->num);
X    freeh(tmp.v);
X    return r;
X}
X
X
X/*
X * Return whether or not a number is an exact square.
X */
XBOOL
Xqissquare(q)
X    NUMBER *q;
X{
X    BOOL flag;
X
X    flag = zissquare(q->num);
X    if (qisint(q) || !flag)
X        return flag;
X    return zissquare(q->den);
X}
X
X
X/*
X * Compute the greatest integer of the Kth root of a number.
X * Example:  qiroot(85, 3) = 4.
X */
XNUMBER *
Xqiroot(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    NUMBER *r;
X    ZVALUE tmp;
X
X    if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
X        error("Taking number to bad root value");
X    if (qiszero(q1))
X        return qlink(&_qzero_);
X    if (qisone(q1) || qisone(q2))
X        return qlink(q1);
X    if (qistwo(q2))
X        return qisqrt(q1);
X    r = qalloc();
X    if (qisint(q1)) {
X        zroot(q1->num, q2->num, &r->num);
X        return r;
X    }
X    zquo(q1->num, q1->den, &tmp);
X    zroot(tmp, q2->num, &r->num);
X    freeh(tmp.v);
X    return r;
X}
X
X
X/*
X * Return the greatest integer of the base 2 log of a number.
X * This is the number such that  1 <= x / log2(x) < 2.
X * Examples:  qilog2(8) = 3, qilog2(1.3) = 1, qilog2(1/7) = -3.
X */
Xlong
Xqilog2(q)
X    NUMBER *q;        /* number to take log of */
X{
X    long n;            /* power of two */
X    int c;            /* result of comparison */
X    ZVALUE tmp;        /* temporary value */
X
X    if (qisneg(q) || qiszero(q))
X        error("Non-positive number for log2");
X    if (qisint(q))
X        return zhighbit(q->num);
X    n = zhighbit(q->num) - zhighbit(q->den);
X    if (n == 0)
X        c = zrel(q->num, q->den);
X    else if (n > 0) {
X        zshift(q->den, n, &tmp);
X        c = zrel(q->num, tmp);
X    } else {
X        zshift(q->num, n, &tmp);
X        c = zrel(tmp, q->den);
X    }
X    if (n)
X        freeh(tmp.v);
X    if (c < 0)
X        n--;
X    return n;
X}
X
X
X/*
X * Return the greatest integer of the base 10 log of a number.
X * This is the number such that  1 <= x / log10(x) < 10.
X * Examples:  qilog10(100) = 2, qilog10(12.3) = 1, qilog10(.023) = -2.
X */
Xlong
Xqilog10(q)
X    NUMBER *q;        /* number to take log of */
X{
X    long n;            /* log value */
X    ZVALUE temp;        /* temporary value */
X
X    if (qisneg(q) || qiszero(q))
X        error("Non-positive number for log10");
X    if (qisint(q))
X        return zlog10(q->num);
X    /*
X     * The number is not an integer.
X     * Compute the result if the number is greater than one.
X     */
X    if ((q->num.len > q->den.len) ||
X        ((q->num.len == q->den.len) && (zrel(q->num, q->den) > 0))) {
X            zquo(q->num, q->den, &temp);
X            n = zlog10(temp);
X            freeh(temp.v);
X            return n;
X    }
X    /*
X     * Here if the number is less than one.
X     * If the number is the inverse of a power of ten, then the obvious answer
X     * will be off by one.  Subtracting one if the number is the inverse of an
X     * integer will fix it.
X     */
X    if (isunit(q->num))
X        zsub(q->den, _one_, &temp);
X    else
X        zquo(q->den, q->num, &temp);
X    n = -zlog10(temp) - 1;
X    freeh(temp.v);
X    return n;
X}
X
X
X/*
X * Return the number of digits in a number, ignoring the sign.
X * For fractions, this is the number of digits of its greatest integer.
X * Examples: qdigits(3456) = 4, qdigits(-23.45) = 2, qdigits(.0120) = 1.
X */
Xlong
Xqdigits(q)
X    NUMBER *q;        /* number to count digits of */
X{
X    long n;            /* number of digits */
X    ZVALUE temp;        /* temporary value */
X
X    if (qisint(q))
X        return zdigits(q->num);
X    zquo(q->num, q->den, &temp);
X    n = zdigits(temp);
X    freeh(temp.v);
X    return n;
X}
X
X
X/*
X * Return the digit at the specified decimal place of a number represented
X * in floating point.  The lowest digit of the integral part of a number
X * is the zeroth decimal place.  Negative decimal places indicate digits
X * to the right of the decimal point.  Examples: qdigit(1234.5678, 1) = 3,
X * qdigit(1234.5678, -3) = 7.
X */
XFLAG
Xqdigit(q, n)
X    NUMBER *q;
X    long n;
X{
X    ZVALUE tenpow, tmp1, tmp2;
X    FLAG res;
X
X    /*
X     * Zero number or negative decimal place of integer is trivial.
X     */
X    if (qiszero(q) || (qisint(q) && (n < 0)))
X        return 0;
X    /*
X     * For non-negative decimal places, answer is easy.
X     */
X    if (n >= 0) {
X        if (qisint(q))
X            return zdigit(q->num, n);
X        zquo(q->num, q->den, &tmp1);
X        res = zdigit(tmp1, n);
X        freeh(tmp1.v);
X        return res;
X    }
X    /*
X     * Fractional value and want negative digit, must work harder.
X     */
X    ztenpow(-n, &tenpow);
X    zmul(q->num, tenpow, &tmp1);
X    freeh(tenpow.v);
X    zquo(tmp1, q->den, &tmp2);
X    res = zmodi(tmp2, 10L);
X    freeh(tmp1.v);
X    freeh(tmp2.v);
X    return res;
X}
X
X
X/*
X * Return whether or not a bit is set at the specified bit position in a
X * number.  The lowest bit of the integral part of a number is the zeroth
X * bit position.  Negative bit positions indicate bits to the right of the
X * binary decimal point.  Examples: qdigit(17.1, 0) = 1, qdigit(17.1, -1) = 0.
X */
XBOOL
Xqisset(q, n)
X    NUMBER *q;
X    long n;
X{
X    NUMBER *qtmp1, *qtmp2;
X    ZVALUE ztmp;
X    BOOL res;
X
X    /*
X     * Zero number or negative bit position place of integer is trivial.
X     */
X    if (qiszero(q) || (qisint(q) && (n < 0)))
X        return FALSE;
X    /*
X     * For non-negative bit positions, answer is easy.
X     */
X    if (n >= 0) {
X        if (qisint(q))
X            return zisset(q->num, n);
X        zquo(q->num, q->den, &ztmp);
X        res = zisset(ztmp, n);
X        freeh(ztmp.v);
X        return res;
X    }
X    /*
X     * Fractional value and want negative bit position, must work harder.
X     */
X    qtmp1 = qscale(q, -n);
X    qtmp2 = qint(qtmp1);
X    qfree(qtmp1);
X    res = ((qtmp2->num.v[0] & 0x01) != 0);
X    qfree(qtmp2);
X    return res;
X}
X
X
X/*
X * Compute the factorial of an integer.
X *    q2 = qfact(q1);
X */
XNUMBER *
Xqfact(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (qisfrac(q))
X        error("Non-integral factorial");
X    if (qiszero(q) || isone(q->num))
X        return qlink(&_qone_);
X    r = qalloc();
X    zfact(q->num, &r->num);
X    return r;
X}
X
X
X/*
X * Compute the product of the primes less than or equal to a number.
X *    q2 = qpfact(q1);
X */
XNUMBER *
Xqpfact(q)
X    register NUMBER *q;
X{
X    NUMBER *r;
X
X    if (qisfrac(q))
X        error("Non-integral factorial");
X    r = qalloc();
X    zpfact(q->num, &r->num);
X    return r;
X}
X
X
X/*
X * Compute the lcm of all the numbers less than or equal to a number.
X *    q2 = qlcmfact(q1);
X */
XNUMBER *
Xqlcmfact(q)
X    register NUMBER *q;
X{
X    NUMBER *r;
X
X    if (qisfrac(q))
X        error("Non-integral lcmfact");
X    r = qalloc();
X    zlcmfact(q->num, &r->num);
X    return r;
X}
X
X
X/*
X * Compute the permutation function  M! / (M - N)!.
X */
XNUMBER *
Xqperm(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    register NUMBER *r;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integral arguments for permutation");
X    r = qalloc();
X    zperm(q1->num, q2->num, &r->num);
X    return r;
X}
X
X
X/*
X * Compute the combinatorial function  M! / (N! * (M - N)!).
X */
XNUMBER *
Xqcomb(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    register NUMBER *r;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integral arguments for combinatorial");
X    r = qalloc();
X    zcomb(q1->num, q2->num, &r->num);
X    return r;
X}
X
X
X/*
X * Compute the Jacobi function (a / b).
X * -1 => a is not quadratic residue mod b
X * 1 => b is composite, or a is quad residue of b
X * 0 => b is even or b < 0
X */
XNUMBER *
Xqjacobi(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integral arguments for jacobi");
X    return itoq((long) zjacobi(q1->num, q2->num));
X}
X
X
X/*
X * Compute the Fibonacci number F(n).
X */
XNUMBER *
Xqfib(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (qisfrac(q))
X        error("Non-integral Fibonacci number");
X    r = qalloc();
X    zfib(q->num, &r->num);
X    return r;
X}
X
X
X/*
X * Truncate a number to the specified number of decimal places.
X * Specifying zero places makes the result identical to qint.
X * Example: qtrunc(2/3, 3) = .666
X */
XNUMBER *
Xqtrunc(q1, q2)
X    NUMBER *q1, *q2;
X{
X    long places;
X    NUMBER *r;
X    ZVALUE tenpow, tmp1, tmp2;
X
X    if (qisfrac(q2) || qisneg(q2) || !istiny(q2->num))
X        error("Bad number of places for qtrunc");
X    if (qisint(q1))
X        return qlink(q1);
X    r = qalloc();
X    places = z1tol(q2->num);
X    /*
X     * Ok, produce the result.
X     * First see if we want no places, in which case just take integer part.
X     */
X    if (places == 0) {
X        zquo(q1->num, q1->den, &r->num);
X        return r;
X    }
X    ztenpow(places, &tenpow);
X    zmul(q1->num, tenpow, &tmp1);
X    zquo(tmp1, q1->den, &tmp2);
X    freeh(tmp1.v);
X    if (iszero(tmp2)) {
X        freeh(tmp2.v);
X        return qlink(&_qzero_);
X    }
X    /*
X     * Now reduce the result to the lowest common denominator.
X     */
X    zgcd(tmp2, tenpow, &tmp1);
X    if (isunit(tmp1)) {
X        freeh(tmp1.v);
X        r->num = tmp2;
X        r->den = tenpow;
X        return r;
X    }
X    zquo(tmp2, tmp1, &r->num);
X    zquo(tenpow, tmp1, &r->den);
X    freeh(tmp1.v);
X    freeh(tmp2.v);
X    freeh(tenpow.v);
X    return r;
X}
X
X
X/*
X * Round a number to the specified number of decimal places.
X * Zero decimal places means round to the nearest integer.
X * Example: qround(2/3, 3) = .667
X */
XNUMBER *
Xqround(q, places)
X    NUMBER *q;        /* number to be rounded */
X    long places;        /* number of decimal places to round to */
X{
X    NUMBER *r;
X    ZVALUE tenpow, roundval, tmp1, tmp2;
X
X    if (places < 0)
X        error("Negative places for qround");
X    if (qisint(q))
X        return qlink(q);
X    /*
X     * Calculate one half of the denominator, ignoring fractional results.
X     * This is the value we will add in order to cause rounding.
X     */
X    zshift(q->den, -1L, &roundval);
X    roundval.sign = q->num.sign;
X    /*
X     * Ok, now do the actual work to produce the result.
X     */
X    r = qalloc();
X    ztenpow(places, &tenpow);
X    zmul(q->num, tenpow, &tmp2);
X    zadd(tmp2, roundval, &tmp1);
X    freeh(tmp2.v);
X    freeh(roundval.v);
X    zquo(tmp1, q->den, &tmp2);
X    freeh(tmp1.v);
X    if (iszero(tmp2)) {
X        freeh(tmp2.v);
X        return qlink(&_qzero_);
X    }
X    /*
X     * Now reduce the result to the lowest common denominator.
X     */
X    zgcd(tmp2, tenpow, &tmp1);
X    if (isunit(tmp1)) {
X        freeh(tmp1.v);
X        r->num = tmp2;
X        r->den = tenpow;
X        return r;
X    }
X    zquo(tmp2, tmp1, &r->num);
X    zquo(tenpow, tmp1, &r->den);
X    freeh(tmp1.v);
X    freeh(tmp2.v);
X    freeh(tenpow.v);
X    return r;
X}
X
X
X/*
X * Truncate a number to the specified number of binary places.
X * Specifying zero places makes the result identical to qint.
X */
XNUMBER *
Xqbtrunc(q1, q2)
X    NUMBER *q1, *q2;
X{
X    long places, twopow;
X    NUMBER *r;
X    ZVALUE tmp1, tmp2;
X
X    if (qisfrac(q2) || qisneg(q2) || !istiny(q2->num))
X        error("Bad number of places for qtrunc");
X    if (qisint(q1))
X        return qlink(q1);
X    r = qalloc();
X    places = z1tol(q2->num);
X    /*
X     * Ok, produce the result.
X     * First see if we want no places, in which case just take integer part.
X     */
X    if (places == 0) {
X        zquo(q1->num, q1->den, &r->num);
X        return r;
X    }
X    zshift(q1->num, places, &tmp1);
X    zquo(tmp1, q1->den, &tmp2);
X    freeh(tmp1.v);
X    if (iszero(tmp2)) {
X        freeh(tmp2.v);
X        return qlink(&_qzero_);
X    }
X    /*
X     * Now reduce the result to the lowest common denominator.
X     */
X    if (isodd(tmp2)) {
X        r->num = tmp2;
X        zbitvalue(places, &r->den);
X        return r;
X    }
X    twopow = zlowbit(tmp2);
X    if (twopow > places)
X        twopow = places;
X    places -= twopow;
X    zshift(tmp2, -twopow, &r->num);
X    freeh(tmp2.v);
X    zbitvalue(places, &r->den);
X    return r;
X}
X
X
X/*
X * Round a number to the specified number of binary places.
X * Zero binary places means round to the nearest integer.
X */
XNUMBER *
Xqbround(q, places)
X    NUMBER *q;        /* number to be rounded */
X    long places;        /* number of binary places to round to */
X{
X    long twopow;
X    NUMBER *r;
X    ZVALUE roundval, tmp1, tmp2;
X
X    if (places < 0)
X        error("Negative places for qbround");
X    if (qisint(q))
X        return qlink(q);
X    r = qalloc();
X    /*
X     * Calculate one half of the denominator, ignoring fractional results.
X     * This is the value we will add in order to cause rounding.
X     */
X    zshift(q->den, -1L, &roundval);
X    roundval.sign = q->num.sign;
X    /*
X     * Ok, produce the result.
X     */
X    zshift(q->num, places, &tmp1);
X    zadd(tmp1, roundval, &tmp2);
X    freeh(roundval.v);
X    freeh(tmp1.v);
X    zquo(tmp2, q->den, &tmp1);
X    freeh(tmp2.v);
X    if (iszero(tmp1)) {
X        freeh(tmp1.v);
X        return qlink(&_qzero_);
X    }
X    /*
X     * Now reduce the result to the lowest common denominator.
X     */
X    if (isodd(tmp1)) {
X        r->num = tmp1;
X        zbitvalue(places, &r->den);
X        return r;
X    }
X    twopow = zlowbit(tmp1);
X    if (twopow > places)
X        twopow = places;
X    places -= twopow;
X    zshift(tmp1, -twopow, &r->num);
X    freeh(tmp1.v);
X    zbitvalue(places, &r->den);
X    return r;
X}
X
X
X/*
X * Approximate a number by using binary rounding with the minimum number
X * of binary places so that the resulting number is within the specified
X * epsilon of the original number.
X */
XNUMBER *
Xqbappr(q, e)
X    NUMBER *q, *e;
X{
X    long bits;
X
X    if (qisneg(e) || qiszero(e))
X        error("Bad epsilon value for qbappr");
X    if (e == _epsilon_)
X        bits = _epsilonprec_ + 1;
X    else
X        bits = qprecision(e) + 1;
X    return qbround(q, bits);
X}
X
X
X/*
X * Approximate a number using continued fractions until the approximation
X * error is less than the specified value.  If a NULL pointer is given
X * for the error value, then the closest simpler fraction is returned.
X */
XNUMBER *
Xqcfappr(q, e)
X    NUMBER *q, *e;
X{
X    NUMBER qtest, *qtmp;
X    ZVALUE u1, u2, u3, v1, v2, v3, t1, t2, t3, qq, tt;
X    int i;
X    BOOL haveeps;
X
X    haveeps = TRUE;
X    if (e == NULL) {
X        haveeps = FALSE;
X        e = &_qzero_;
X    }
X    if (qisneg(e))
X        error("Negative epsilon for cfappr");
X    if (qisint(q) || isunit(q->num) || (haveeps && qiszero(e)))
X        return qlink(q);
X    u1 = _one_;
X    u2 = _zero_;
X    u3 = q->num;
X    u3.sign = 0;
X    v1 = _zero_;
X    v2 = _one_;
X    v3 = q->den;
X    while (!iszero(v3)) {
X        if (!iszero(u2) && !iszero(u1)) {
X            qtest.num = u2;
X            qtest.den = u1;
X            qtest.den.sign = 0;
X            qtest.num.sign = q->num.sign;
X            qtmp = qsub(q, &qtest);
X            qtest = *qtmp;
X            qtest.num.sign = 0;
X            i = qrel(&qtest, e);
X            qfree(qtmp);
X            if (i <= 0)
X                break;
X        }
X        zquo(u3, v3, &qq);
X        zmul(qq, v1, &tt); zsub(u1, tt, &t1); freeh(tt.v);
X        zmul(qq, v2, &tt); zsub(u2, tt, &t2); freeh(tt.v);
X        zmul(qq, v3, &tt); zsub(u3, tt, &t3); freeh(tt.v);
X        freeh(qq.v); freeh(u1.v); freeh(u2.v);
X        if ((u3.v != q->num.v) && (u3.v != q->den.v))
X            freeh(u3.v);
X        u1 = v1; u2 = v2; u3 = v3;
X        v1 = t1; v2 = t2; v3 = t3;
X    }
X    if (u3.v != q->den.v)
X        freeh(u3.v);
X    freeh(v1.v);
X    freeh(v2.v);
X    i = iszero(v3);
X    freeh(v3.v);
X    if (i && haveeps) {
X        freeh(u1.v);
X        freeh(u2.v);
X        return qlink(q);
X    }
X    qtest.num = u2;
X    qtest.den = u1;
X    qtest.den.sign = 0;
X    qtest.num.sign = q->num.sign;
X    qtmp = qcopy(&qtest);
X    freeh(u1.v);
X    freeh(u2.v);
X    return qtmp;
X}
X
X
X/*
X * Return an indication on whether or not two fractions are approximately
X * equal within the specified epsilon. Returns negative if the absolute value
X * of the difference between the two numbers is less than epsilon, zero if
X * the difference is equal to epsilon, and positive if the difference is
X * greater than epsilon.
X */
XFLAG
Xqnear(q1, q2, epsilon)
X    NUMBER *q1, *q2, *epsilon;
X{
X    int res;
X    NUMBER qtemp, *qq;
X
X    if (qisneg(epsilon))
X        error("Negative epsilon for near");
X    if (q1 == q2) {
X        if (qiszero(epsilon))
X            return 0;
X        return -1;
X    }
X    if (qiszero(epsilon))
X        return qcmp(q1, q2);
X    if (qiszero(q2)) {
X        qtemp = *q1;
X        qtemp.num.sign = 0;
X        return qrel(&qtemp, epsilon);
X    }
X    if (qiszero(q1)) {
X        qtemp = *q2;
X        qtemp.num.sign = 0;
X        return qrel(&qtemp, epsilon);
X    }
X    qq = qsub(q1, q2);
X    qtemp = *qq;
X    qtemp.num.sign = 0;
X    res = qrel(&qtemp, epsilon);
X    qfree(qq);
X    return res;
X}
X
X
X/*
X * Compute the gcd (greatest common divisor) of two numbers.
X *    q3 = qgcd(q1, q2);
X */
XNUMBER *
Xqgcd(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    ZVALUE z;
X    NUMBER *q;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for gcd");
X    zgcd(q1->num, q2->num, &z);
X    if (isunit(z)) {
X        freeh(z.v);
X        return qlink(&_qone_);
X    }
X    q = qalloc();
X    q->num = z;
X    return q;
X}
X
X
X/*
X * Compute the lcm (least common denominator) of two numbers.
X *    q3 = qlcm(q1, q2);
X */
XNUMBER *
Xqlcm(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    NUMBER *q;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for lcm");
X    if (qisunit(q1))
X        return qlink(q2);
X    if (qisunit(q2))
X        return qlink(q1);
X    q = qalloc();
X    zlcm(q1->num, q2->num, &q->num);
X    return q;
X}
X
X
X/*
X * Remove all occurances of the specified factor from a number.
X * Returned number is always positive.
X */
XNUMBER *
Xqfacrem(q1, q2)
X    NUMBER *q1, *q2;
X{
X    long count;
X    ZVALUE tmp;
X    NUMBER *r;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for factor removal");
X    count = zfacrem(q1->num, q2->num, &tmp);
X    if (isunit(tmp)) {
X        freeh(tmp.v);
X        return qlink(&_qone_);
X    }
X    if (count == 0) {
X        freeh(tmp.v);
X        return qlink(q1);
X    }
X    r = qalloc();
X    r->num = tmp;
X    return r;
X}
X
X
X/*
X * Divide one number by the gcd of it with another number repeatedly until
X * the number is relatively prime.
X */
XNUMBER *
Xqgcdrem(q1, q2)
X    NUMBER *q1, *q2;
X{
X    ZVALUE tmp;
X    NUMBER *r;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for gcdrem");
X    zgcdrem(q1->num, q2->num, &tmp);
X    if (isunit(tmp)) {
X        freeh(tmp.v);
X        return qlink(&_qone_);
X    }
X    if (zcmp(q1->num, tmp) == 0) {
X        freeh(tmp.v);
X        return qlink(q1);
X    }
X    r = qalloc();
X    r->num = tmp;
X    return r;
X}
X
X
X/*
X * Return the lowest prime factor of a number.
X * Search is conducted for the specified number of primes.
X * Returns one if no factor was found.
X */
XNUMBER *
Xqlowfactor(q1, q2)
X    NUMBER *q1, *q2;
X{
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for lowfactor");
X    return itoq(zlowfactor(q1->num, ztoi(q2->num)));
X}
X
X
X/*
X * Return the number of places after the decimal point needed to exactly
X * represent the specified number as a real number.  Integers return zero,
X * and non-terminating decimals return minus one.  Examples:
X *    qplaces(1/7)=-1, qplaces(3/10)= 1, qplaces(1/8)=3, qplaces(4)=0.
X */
Xlong
Xqplaces(q)
X    NUMBER *q;
X{
X    long twopow, fivepow;
X    HALF fiveval[2];
X    ZVALUE five;
X    ZVALUE tmp;
X
X    if (qisint(q))    /* no decimal places if number is integer */
X        return 0;
X    /*
X     * The number of decimal places of a fraction in lowest terms is finite
X     * if an only if the denominator is of the form 2^A * 5^B, and then the
X     * number of decimal places is equal to MAX(A, B).
X     */
X    five.sign = 0;
X    five.len = 1;
X    five.v = fiveval;
X    fiveval[0] = 5;
X    fivepow = zfacrem(q->den, five, &tmp);
X    if (!zisonebit(tmp)) {
X        freeh(tmp.v);
X        return -1;
X    }
X    twopow = zlowbit(tmp);
X    freeh(tmp.v);
X    if (twopow < fivepow)
X        twopow = fivepow;
X    return twopow;
X}
X
X
X/*
X * Perform a probabilistic primality test (algorithm P in Knuth).
X * Returns FALSE if definitely not prime, or TRUE if probably prime.
X * Second arg determines how many times to check for primality.
X */
XBOOL
Xqprimetest(q1, q2)
X    NUMBER *q1, *q2;
X{
X    if (qisfrac(q1) || qisfrac(q2) || qisneg(q2))
X        error("Bad arguments for qprimetest");
X    return zprimetest(q1->num, qtoi(q2));
X}
X
X/* END CODE */
END_OF_FILE
if test 24256 -ne `wc -c <'qfunc.c'`; then
    echo shar: \"'qfunc.c'\" unpacked with wrong size!
fi
# end of 'qfunc.c'
fi
echo shar: End of archive 10 \(of 21\).
cp /dev/null ark10isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 21 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0