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

Text File | 1992-05-09 | 43.8 KB | 2,181 lines

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

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

#! /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 9 (of 21)."
# Contents:  qmath.c qtrans.c
# Wrapped by dbell@elm on Tue Feb 25 15:21:04 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'qmath.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'qmath.c'\"
else
echo shar: Extracting \"'qmath.c'\" \(20144 characters\)
sed "s/^X//" >'qmath.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 primitive routines
X */
X
X#include <stdio.h>
X#include "math.h"
X
X
XNUMBER _qzero_ =    { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
XNUMBER _qone_ =        { { _oneval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
Xstatic NUMBER _qtwo_ =    { { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
Xstatic NUMBER _qten_ =    { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
XNUMBER _qnegone_ =    { { _oneval_, 1, 1 }, { _oneval_, 1, 0 }, 1 };
XNUMBER _qonehalf_ =    { { _oneval_, 1, 0 }, { _twoval_, 1, 0 }, 1 };
X
X
X#if 0
Xstatic char *abortmsg = "Calculation aborted";
X#endif
Xstatic char *memmsg = "Not enough memory";
X
X
X/*
X * Create another copy of a number.
X *    q2 = qcopy(q1);
X */
XNUMBER *
Xqcopy(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    r = qalloc();
X    r->num.sign = q->num.sign;
X    if (!isunit(q->num)) {
X        r->num.len = q->num.len;
X        r->num.v = alloc(r->num.len);
X        copyval(q->num, r->num);
X    }
X    if (!isunit(q->den)) {
X        r->den.len = q->den.len;
X        r->den.v = alloc(r->den.len);
X        copyval(q->den, r->den);
X    }
X    return r;
X}
X
X
X/*
X * Convert a number to a normal integer.
X *    i = qtoi(q);
X */
Xlong
Xqtoi(q)
X    register NUMBER *q;
X{
X    long i;
X    ZVALUE res;
X
X    if (qisint(q))
X        return ztoi(q->num);
X    zquo(q->num, q->den, &res);
X    i = ztoi(res);
X    freeh(res.v);
X    return i;
X}
X
X
X/*
X * Convert a normal integer into a number.
X *    q = itoq(i);
X */
XNUMBER *
Xitoq(i)
X    long i;
X{
X    register NUMBER *q;
X
X    if ((i >= -1) && (i <= 10)) {
X        switch ((int) i) {
X            case 0: q = &_qzero_; break;
X            case 1: q = &_qone_; break;
X            case 2: q = &_qtwo_; break;
X            case 10: q = &_qten_; break;
X            case -1: q = &_qnegone_; break;
X            default: q = NULL;
X        }
X        if (q)
X            return qlink(q);
X    }
X    q = qalloc();
X    itoz(i, &q->num);
X    return q;
X}
X
X
X/*
X * Create a number from the given integral numerator and denominator.
X *    q = iitoq(inum, iden);
X */
XNUMBER *
Xiitoq(inum, iden)
X    long inum, iden;
X{
X    register NUMBER *q;
X    long d;
X    BOOL sign;
X
X    if (iden == 0)
X        error("Division by zero");
X    if (inum == 0)
X        return qlink(&_qzero_);
X    sign = 0;
X    if (inum < 0) {
X        sign = 1;
X        inum = -inum;
X    }
X    if (iden < 0) {
X        sign = 1 - sign;
X        iden = -iden;
X    }
X    d = iigcd(inum, iden);
X    inum /= d;
X    iden /= d;
X    if (iden == 1)
X        return itoq(sign ? -inum : inum);
X    q = qalloc();
X    if (inum != 1)
X        itoz(inum, &q->num);
X    itoz(iden, &q->den);
X    q->num.sign = sign;
X    return q;
X}
X
X
X/*
X * Add two numbers to each other.
X *    q3 = qadd(q1, q2);
X */
XNUMBER *
Xqadd(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    NUMBER *r;
X    ZVALUE t1, t2, temp, d1, d2, vpd1, upd1;
X
X    r = qalloc();
X    /*
X     * If either number is an integer, then the result is easy.
X     */
X    if (qisint(q1) && qisint(q2)) {
X        zadd(q1->num, q2->num, &r->num);
X        return r;
X    }
X    if (qisint(q2)) {
X        zmul(q1->den, q2->num, &temp);
X        zadd(q1->num, temp, &r->num);
X        freeh(temp.v);
X        zcopy(q1->den, &r->den);
X        return r;
X    }
X    if (qisint(q1)) {
X        zmul(q2->den, q1->num, &temp);
X        zadd(q2->num, temp, &r->num);
X        freeh(temp.v);
X        zcopy(q2->den, &r->den);
X        return r;
X    }
X    /*
X     * Both arguments are true fractions, so we need more work.
X     * If the denominators are relatively prime, then the answer is the
X     * straightforward cross product result with no need for reduction.
X     */
X    zgcd(q1->den, q2->den, &d1);
X    if (isunit(d1)) {
X        freeh(d1.v);
X        zmul(q1->num, q2->den, &t1);
X        zmul(q1->den, q2->num, &t2);
X        zadd(t1, t2, &r->num);
X        freeh(t1.v);
X        freeh(t2.v);
X        zmul(q1->den, q2->den, &r->den);
X        return r;
X    }
X    /*
X     * The calculation is now more complicated.
X     * See Knuth Vol 2 for details.
X     */
X    zquo(q2->den, d1, &vpd1);
X    zquo(q1->den, d1, &upd1);
X    zmul(q1->num, vpd1, &t1);
X    zmul(q2->num, upd1, &t2);
X    zadd(t1, t2, &temp);
X    freeh(t1.v);
X    freeh(t2.v);
X    freeh(vpd1.v);
X    zgcd(temp, d1, &d2);
X    freeh(d1.v);
X    if (isunit(d2)) {
X        freeh(d2.v);
X        r->num = temp;
X        zmul(upd1, q2->den, &r->den);
X        freeh(upd1.v);
X        return r;
X    }
X    zquo(temp, d2, &r->num);
X    freeh(temp.v);
X    zquo(q2->den, d2, &temp);
X    freeh(d2.v);
X    zmul(temp, upd1, &r->den);
X    freeh(temp.v);
X    freeh(upd1.v);
X    return r;
X}
X
X
X/*
X * Subtract one number from another.
X *    q3 = qsub(q1, q2);
X */
XNUMBER *
Xqsub(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    NUMBER t, *r;
X
X    if (q1 == q2)
X        return qlink(&_qzero_);
X    if (qisint(q1) && qisint(q2)) {
X        r = qalloc();
X        zsub(q1->num, q2->num, &r->num);
X        return r;
X    }
X    t = *q2;
X    t.num.sign = !t.num.sign;
X    return qadd(q1, &t);
X}
X
X
X/*
X * Increment a number by one.
X */
XNUMBER *
Xqinc(q)
X    NUMBER *q;
X{
X    NUMBER *r;
X
X    r = qalloc();
X    if (qisint(q)) {
X        zadd(q->num, _one_, &r->num);
X        return r;
X    }
X    zadd(q->num, q->den, &r->num);
X    zcopy(q->den, &r->den);
X    return r;
X}
X
X
X/*
X * Decrement a number by one.
X */
XNUMBER *
Xqdec(q)
X    NUMBER *q;
X{
X    NUMBER *r;
X
X    r = qalloc();
X    if (qisint(q)) {
X        zsub(q->num, _one_, &r->num);
X        return r;
X    }
X    zsub(q->num, q->den, &r->num);
X    zcopy(q->den, &r->den);
X    return r;
X}
X
X
X#ifdef CODE
X/*
X * Add a normal small integer value to an arbitrary number.
X */
XNUMBER *
Xqaddi(q1, n)
X    NUMBER *q1;
X    long n;
X{
X    NUMBER addnum;        /* temporary number */
X    HALF addval[2];        /* value of small number */
X    BOOL neg;        /* TRUE if number is neg */
X
X    if (n == 0)
X        return qlink(q1);
X    if (n == 1)
X        return qinc(q1);
X    if (n == -1)
X        return qdec(q1);
X    if (qiszero(q1))
X        return itoq(n);
X    addnum.num.sign = 0;
X    addnum.num.len = 1;
X    addnum.num.v = addval;
X    addnum.den = _one_;
X    neg = (n < 0);
X    if (neg)
X        n = -n;
X    addval[0] = (HALF) n;
X    n = (((FULL) n) >> BASEB);
X    if (n) {
X        addval[1] = (HALF) n;
X        addnum.num.len = 2;
X    }
X    if (neg)
X        return qsub(q1, &addnum);
X    else
X        return qadd(q1, &addnum);
X}
X#endif
X
X
X/*
X * Multiply two numbers.
X *    q3 = qmul(q1, q2);
X */
XNUMBER *
Xqmul(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    NUMBER *r;            /* returned value */
X    ZVALUE n1, n2, d1, d2;        /* numerators and denominators */
X    ZVALUE tmp;
X
X    if (qiszero(q1) || qiszero(q2))
X        return qlink(&_qzero_);
X    if (qisone(q1))
X        return qlink(q2);
X    if (qisone(q2))
X        return qlink(q1);
X    if (qisint(q1) && qisint(q2)) {    /* easy results if integers */
X        r = qalloc();
X        zmul(q1->num, q2->num, &r->num);
X        return r;
X    }
X    n1 = q1->num;
X    n2 = q2->num;
X    d1 = q1->den;
X    d2 = q2->den;
X    if (iszero(d1) || iszero(d2))
X        error("Division by zero");
X    if (iszero(n1) || iszero(n2))
X        return qlink(&_qzero_);
X    if (!isunit(n1) && !isunit(d2)) {    /* possibly reduce */
X        zgcd(n1, d2, &tmp);
X        if (!isunit(tmp)) {
X            zquo(q1->num, tmp, &n1);
X            zquo(q2->den, tmp, &d2);
X        }
X        freeh(tmp.v);
X    }
X    if (!isunit(n2) && !isunit(d1)) {    /* again possibly reduce */
X        zgcd(n2, d1, &tmp);
X        if (!isunit(tmp)) {
X            zquo(q2->num, tmp, &n2);
X            zquo(q1->den, tmp, &d1);
X        }
X        freeh(tmp.v);
X    }
X    r = qalloc();
X    zmul(n1, n2, &r->num);
X    zmul(d1, d2, &r->den);
X    if (q1->num.v != n1.v)
X        freeh(n1.v);
X    if (q1->den.v != d1.v)
X        freeh(d1.v);
X    if (q2->num.v != n2.v)
X        freeh(n2.v);
X    if (q2->den.v != d2.v)
X        freeh(d2.v);
X    return r;
X}
X
X
X#if 0
X/*
X * Multiply a number by a small integer.
X *    q2 = qmuli(q1, n);
X */
XNUMBER *
Xqmuli(q, n)
X    NUMBER *q;
X    long n;
X{
X    NUMBER *r;
X    long d;            /* gcd of multiplier and denominator */
X    int sign;
X
X    if ((n == 0) || qiszero(q))
X        return qlink(&_qzero_);
X    if (n == 1)
X        return qlink(q);
X    r = qalloc();
X    if (qisint(q)) {
X        zmuli(q->num, n, &r->num);
X        return r;
X    }
X    sign = 1;
X    if (n < 0) {
X        n = -n;
X        sign = -1;
X    }
X    d = zmodi(q->den, n);
X    d = iigcd(d, n);
X    zmuli(q->num, (n * sign) / d, &r->num);
X    (void) zdivi(q->den, d, &r->den);
X    return r;
X}
X#endif
X
X
X/*
X * Divide two numbers (as fractions).
X *    q3 = qdiv(q1, q2);
X */
XNUMBER *
Xqdiv(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    NUMBER temp;
X
X    if (q1 == q2) {
X        if (qiszero(q2))
X            error("Division by zero");
X        return qlink(&_qone_);
X    }
X    if (qisone(q1))
X        return qinv(q2);
X    temp.num = q2->den;
X    temp.den = q2->num;
X    temp.num.sign = temp.den.sign;
X    temp.den.sign = 0;
X    temp.links = 1;
X    return qmul(q1, &temp);
X}
X
X
X/*
X * Divide a number by a small integer.
X *    q2 = qdivi(q1, n);
X */
XNUMBER *
Xqdivi(q, n)
X    NUMBER *q;
X    long n;
X{
X    NUMBER *r;
X    long d;            /* gcd of divisor and numerator */
X    int sign;
X
X    if (n == 0)
X        error("Division by zero");
X    if ((n == 1) || qiszero(q))
X        return qlink(q);
X    sign = 1;
X    if (n < 0) {
X        n = -n;
X        sign = -1;
X    }
X    r = qalloc();
X    d = zmodi(q->num, n);
X    d = iigcd(d, n);
X    (void) zdivi(q->num, d * sign, &r->num);
X    zmuli(q->den, n / d, &r->den);
X    return r;
X}
X
X
X/*
X * Return the quotient when one number is divided by another.
X * This works for fractional values also, and in all cases:
X *    qquo(q1, q2) = int(q1 / q2).
X */
XNUMBER *
Xqquo(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    ZVALUE num, den, res;
X    NUMBER *q;
X
X    if (isunit(q1->num))
X        num = q2->den;
X    else if (isunit(q2->den))
X        num = q1->num;
X    else
X        zmul(q1->num, q2->den, &num);
X    if (isunit(q1->den))
X        den = q2->num;
X    else if (isunit(q2->num))
X        den = q1->den;
X    else
X        zmul(q1->den, q2->num, &den);
X    zquo(num, den, &res);
X    if ((num.v != q2->den.v) && (num.v != q1->num.v))
X        freeh(num.v);
X    if ((den.v != q2->num.v) && (den.v != q1->den.v))
X        freeh(den.v);
X    if (iszero(res)) {
X        freeh(res.v);
X        return qlink(&_qzero_);
X    }
X    if (isunit(res)) {
X        q = (res.sign ? &_qnegone_ : &_qone_);
X        freeh(res.v);
X        return qlink(q);
X    }
X    q = qalloc();
X    q->num = res;
X    return q;
X}
X
X
X/*
X * Return the absolute value of a number.
X *    q2 = qabs(q1);
X */
XNUMBER *
Xqabs(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (q->num.sign == 0)
X        return qlink(q);
X    r = qalloc();
X    if (!isunit(q->num))
X        zcopy(q->num, &r->num);
X    if (!isunit(q->den))
X        zcopy(q->den, &r->den);
X    r->num.sign = 0;
X    return r;
X}
X
X
X/*
X * Negate a number.
X *    q2 = qneg(q1);
X */
XNUMBER *
Xqneg(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (qiszero(q))
X        return qlink(&_qzero_);
X    r = qalloc();
X    if (!isunit(q->num))
X        zcopy(q->num, &r->num);
X    if (!isunit(q->den))
X        zcopy(q->den, &r->den);
X    r->num.sign = !q->num.sign;
X    return r;
X}
X
X
X/*
X * Return the sign of a number (-1, 0 or 1)
X */
XNUMBER *
Xqsign(q)
X    NUMBER *q;
X{
X    if (qiszero(q))
X        return qlink(&_qzero_);
X    if (qisneg(q))
X        return qlink(&_qnegone_);
X    return qlink(&_qone_);
X}
X
X
X/*
X * Invert a number.
X *    q2 = qinv(q1);
X */
XNUMBER *
Xqinv(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (qisunit(q)) {
X        r = (qisneg(q) ? &_qnegone_ : &_qone_);
X        return qlink(r);
X    }
X    if (qiszero(q))
X        error("Division by zero");
X    r = qalloc();
X    if (!isunit(q->num))
X        zcopy(q->num, &r->den);
X    if (!isunit(q->den))
X        zcopy(q->den, &r->num);
X    r->num.sign = q->num.sign;
X    r->den.sign = 0;
X    return r;
X}
X
X
X/*
X * Return just the numerator of a number.
X *    q2 = qnum(q1);
X */
XNUMBER *
Xqnum(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (qisint(q))
X        return qlink(q);
X    if (isunit(q->num)) {
X        r = (qisneg(q) ? &_qnegone_ : &_qone_);
X        return qlink(r);
X    }
X    r = qalloc();
X    zcopy(q->num, &r->num);
X    return r;
X}
X
X
X/*
X * Return just the denominator of a number.
X *    q2 = qden(q1);
X */
XNUMBER *
Xqden(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (qisint(q))
X        return qlink(&_qone_);
X    r = qalloc();
X    zcopy(q->den, &r->num);
X    return r;
X}
X
X
X/*
X * Return the fractional part of a number.
X *    q2 = qfrac(q1);
X */
XNUMBER *
Xqfrac(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (qisint(q))
X        return qlink(&_qzero_);
X    if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
X        (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
X            return qlink(q);
X    r = qalloc();
X    zmod(q->num, q->den, &r->num);
X    zcopy(q->den, &r->den);
X    r->num.sign = q->num.sign;
X    return r;
X}
X
X
X/*
X * Return the integral part of a number.
X *    q2 = qint(q1);
X */
XNUMBER *
Xqint(q)
X    register NUMBER *q;
X{
X    register NUMBER *r;
X
X    if (qisint(q))
X        return qlink(q);
X    if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
X        (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
X            return qlink(&_qzero_);
X    r = qalloc();
X    zquo(q->num, q->den, &r->num);
X    return r;
X}
X
X
X/*
X * Compute the square of a number.
X */
XNUMBER *
Xqsquare(q)
X    register NUMBER *q;
X{
X    ZVALUE num, den;
X
X    if (qiszero(q))
X        return qlink(&_qzero_);
X    if (qisunit(q))
X        return qlink(&_qone_);
X    num = q->num;
X    den = q->den;
X    q = qalloc();
X    if (!isunit(num))
X        zsquare(num, &q->num);
X    if (!isunit(den))
X        zsquare(den, &q->den);
X    return q;
X}
X
X
X/*
X * Shift an integer by a given number of bits. This multiplies the number
X * by the appropriate power of two.  Positive numbers shift left, negative
X * ones shift right.  Low bits are truncated when shifting right.
X */
XNUMBER *
Xqshift(q, n)
X    NUMBER *q;
X    long n;
X{
X    register NUMBER *r;
X
X    if (qisfrac(q))
X        error("Shift of non-integer");
X    if (qiszero(q) || (n == 0))
X        return qlink(q);
X    if (n <= -(q->num.len * BASEB))
X        return qlink(&_qzero_);
X    r = qalloc();
X    zshift(q->num, n, &r->num);
X    return r;
X}
X
X
X/*
X * Scale a number by a power of two, as in:
X *    ans = q * 2^n.
X * This is similar to shifting, except that fractions work.
X */
XNUMBER *
Xqscale(q, pow)
X    NUMBER *q;
X    long pow;
X{
X    long numshift, denshift, tmp;
X    NUMBER *r;
X
X    if (qiszero(q) || (pow == 0))
X        return qlink(q);
X    if ((pow > 1000000L) || (pow < -1000000L))
X        error("Very large scale value");
X    numshift = isodd(q->num) ? 0 : zlowbit(q->num);
X    denshift = isodd(q->den) ? 0 : zlowbit(q->den);
X    if (pow > 0) {
X        tmp = pow;
X        if (tmp > denshift)
X        tmp = denshift;
X        denshift = -tmp;
X        numshift = (pow - tmp);
X    } else {
X        pow = -pow;
X        tmp = pow;
X        if (tmp > numshift)
X            tmp = numshift;
X        numshift = -tmp;
X        denshift = (pow - tmp);
X    }
X    r = qalloc();
X    if (numshift)
X        zshift(q->num, numshift, &r->num);
X    else
X        zcopy(q->num, &r->num);
X    if (denshift)
X        zshift(q->den, denshift, &r->den);
X    else
X        zcopy(q->den, &r->den);
X    return r;
X}
X
X
X/*
X * Return the minimum of two numbers.
X */
XNUMBER *
Xqmin(q1, q2)
X    NUMBER *q1, *q2;
X{
X    if (qrel(q1, q2) > 0)
X        q1 = q2;
X    return qlink(q1);
X}
X
X
X/*
X * Return the maximum of two numbers.
X */
XNUMBER *
Xqmax(q1, q2)
X    NUMBER *q1, *q2;
X{
X    if (qrel(q1, q2) < 0)
X        q1 = q2;
X    return qlink(q1);
X}
X
X
X/*
X * Perform the logical OR of two integers.
X */
XNUMBER *
Xqor(q1, q2)
X    NUMBER *q1, *q2;
X{
X    register NUMBER *r;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for logical or");
X    if ((q1 == q2) || qiszero(q2))
X        return qlink(q1);
X    if (qiszero(q1))
X        return qlink(q2);
X    r = qalloc();
X    zor(q1->num, q2->num, &r->num);
X    return r;
X}
X
X
X/*
X * Perform the logical AND of two integers.
X */
XNUMBER *
Xqand(q1, q2)
X    NUMBER *q1, *q2;
X{
X    register NUMBER *r;
X    ZVALUE res;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for logical and");
X    if (q1 == q2)
X        return qlink(q1);
X    if (qiszero(q1) || qiszero(q2))
X        return qlink(&_qzero_);
X    zand(q1->num, q2->num, &res);
X    if (iszero(res)) {
X        freeh(res.v);
X        return qlink(&_qzero_);
X    }
X    r = qalloc();
X    r->num = res;
X    return r;
X}
X
X
X/*
X * Perform the logical XOR of two integers.
X */
XNUMBER *
Xqxor(q1, q2)
X    NUMBER *q1, *q2;
X{
X    register NUMBER *r;
X    ZVALUE res;
X
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for logical xor");
X    if (q1 == q2)
X        return qlink(&_qzero_);
X    if (qiszero(q1))
X        return qlink(q2);
X    if (qiszero(q2))
X        return qlink(q1);
X    zxor(q1->num, q2->num, &res);
X    if (iszero(res)) {
X        freeh(res.v);
X        return qlink(&_qzero_);
X    }
X    r = qalloc();
X    r->num = res;
X    return r;
X}
X
X
X#if 0
X/*
X * Return the number whose binary representation only has the specified
X * bit set (counting from zero).  This thus produces a given power of two.
X */
XNUMBER *
Xqbitvalue(n)
X    long n;
X{
X    register NUMBER *r;
X
X    if (n <= 0)
X        return qlink(&_qone_);
X    r = qalloc();
X    zbitvalue(n, &r->num);
X    return r;
X}
X
X
X/*
X * Test to see if the specified bit of a number is on (counted from zero).
X * Returns TRUE if the bit is set, or FALSE if it is not.
X *    i = qbittest(q, n);
X */
XBOOL
Xqbittest(q, n)
X    register NUMBER *q;
X    long n;
X{
X    int x, y;
X
X    if ((n < 0) || (n >= (q->num.len * BASEB)))
X        return FALSE;
X    x = q->num.v[n / BASEB];
X    y = (1 << (n % BASEB));
X    return ((x & y) != 0);
X}
X#endif
X
X
X/*
X * Return the precision of a number (usually for examining an epsilon value).
X * This is the largest power of two whose reciprocal is not smaller in absolute
X * value than the specified number.  For example, qbitprec(1/100) = 6.
X * Numbers larger than one have a precision of zero.
X */
Xlong
Xqprecision(q)
X    NUMBER *q;
X{
X    long r;
X
X    if (qisint(q))
X        return 0;
X    if (isunit(q->num))
X        return zhighbit(q->den);
X    r = zhighbit(q->den) - zhighbit(q->num) - 1;
X    if (r < 0)
X        r = 0;
X    return r;
X}
X
X
X#if 0
X/*
X * Return an integer indicating the sign of a number (-1, 0, or 1).
X *    i = qtst(q);
X */
XFLAG
Xqtest(q)
X    register NUMBER *q;
X{
X    if (!ztest(q->num))
X        return 0;
X    if (q->num.sign)
X        return -1;
X    return 1;
X}
X#endif
X
X
X/*
X * Compare two numbers and return an integer indicating their relative size.
X *    i = qrel(q1, q2);
X */
XFLAG
Xqrel(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    ZVALUE z1, z2;
X    long wc1, wc2;
X    int sign;
X    int z1f = 0, z2f = 0;
X
X    if (q1 == q2)
X        return 0;
X    sign = q2->num.sign - q1->num.sign;
X    if (sign)
X        return sign;
X    if (qiszero(q2))
X        return !qiszero(q1);
X    if (qiszero(q1))
X        return -1;
X    /*
X     * Make a quick comparison by calculating the number of words resulting as
X     * if we multiplied through by the denominators, and then comparing the
X     * word counts.
X     */
X    sign = 1;
X    if (qisneg(q1))
X        sign = -1;
X    wc1 = q1->num.len + q2->den.len;
X    wc2 = q2->num.len + q1->den.len;
X    if (wc1 < wc2 - 1)
X        return -sign;
X    if (wc2 < wc1 - 1)
X        return sign;
X    /*
X     * Quick check failed, must actually do the full comparison.
X     */
X    if (isunit(q2->den))
X        z1 = q1->num;
X    else if (isone(q1->num))
X        z1 = q2->den;
X    else {
X        z1f = 1;
X        zmul(q1->num, q2->den, &z1);
X    }
X    if (isunit(q1->den))
X        z2 = q2->num;
X    else if (isone(q2->num))
X        z2 = q1->den;
X    else {
X        z2f = 1;
X        zmul(q2->num, q1->den, &z2);
X    }
X    sign = zrel(z1, z2);
X    if (z1f)
X        freeh(z1.v);
X    if (z2f)
X        freeh(z2.v);
X    return sign;
X}
X
X
X/*
X * Compare two numbers to see if they are equal.
X * This differs from qrel in that the numbers are not ordered.
X * Returns TRUE if they differ.
X */
XBOOL
Xqcmp(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    if (q1 == q2)
X        return FALSE;
X    if ((q1->num.sign != q2->num.sign) || (q1->num.len != q2->num.len) ||
X        (q2->den.len != q2->den.len) || (*q1->num.v != *q2->num.v) ||
X        (*q1->den.v != *q2->den.v))
X            return TRUE;
X    if (zcmp(q1->num, q2->num))
X        return TRUE;
X    if (qisint(q1))
X        return FALSE;
X    return zcmp(q1->den, q2->den);
X}
X
X
X/*
X * Compare a number against a normal small integer.
X * Returns 1, 0, or -1, according to whether the first number is greater,
X * equal, or less than the second number.
X *    n = qreli(q, n);
X */
XFLAG
Xqreli(q, n)
X    NUMBER *q;
X    long n;
X{
X    int sign;
X    ZVALUE num;
X    HALF h2[2];
X    NUMBER q2;
X
X    sign = ztest(q->num);        /* do trivial sign checks */
X    if (sign == 0) {
X        if (n > 0)
X            return -1;
X        return (n < 0);
X    }
X    if ((sign < 0) && (n >= 0))
X        return -1;
X    if ((sign > 0) && (n <= 0))
X        return 1;
X    n *= sign;
X    if (n == 1) {            /* quick check against 1 or -1 */
X        num = q->num;
X        num.sign = 0;
X        return (sign * zrel(num, q->den));
X    }
X    num.sign = (sign < 0);
X    num.len = 1 + (n >= BASE);
X    num.v = h2;
X    h2[0] = (n & BASE1);
X    h2[1] = (n >> BASEB);
X    if (isunit(q->den))    /* integer compare if no denominator */
X        return zrel(q->num, num);
X    q2.num = num;
X    q2.den = _one_;
X    q2.links = 1;
X    return qrel(q, &q2);    /* full fractional compare */
X}
X
X
X/*
X * Compare a number against a small integer to see if they are equal.
X * Returns TRUE if they differ.
X */
XBOOL
Xqcmpi(q, n)
X    NUMBER *q;
X    long n;
X{
X    long len;
X
X    len = q->num.len;
X    if ((len > 2) || qisfrac(q) || (q->num.sign != (n < 0)))
X        return TRUE;
X    if (n < 0)
X        n = -n;
X    if (((HALF)(n)) != q->num.v[0])
X        return TRUE;
X    n = ((FULL) n) >> BASEB;
X    return (((n != 0) != (len == 2)) || (n != q->num.v[1]));
X}
X
X
X/*
X * Number node allocation routines
X */
X
X#define    NNALLOC    1000
X
Xunion allocNode {
X    NUMBER    num;
X    union allocNode    *link;
X};
X
Xstatic union allocNode    *freeNum;
X
X
XNUMBER *
Xqalloc()
X{
X    register union allocNode *temp;
X
X    if (!freeNum) {
X        freeNum = (union allocNode *)
X        malloc(sizeof (NUMBER) * NNALLOC);
X        if (freeNum == NULL)
X            error(memmsg);
X        temp = freeNum;
X        while (temp != freeNum + NNALLOC - 2) {
X            temp->link = temp+1;
X            ++temp;
X        }
X    }
X    temp = freeNum;
X    freeNum = temp->link;
X    temp->num.links = 1;
X    temp->num.num = _one_;
X    temp->num.den = _one_;
X    return &temp->num;
X}
X
X
Xvoid
Xqfreenum(q)
X    register NUMBER *q;
X{
X    union allocNode *a;
X
X    if (q == NULL)
X        return;
X    freeh(q->num.v);
X    freeh(q->den.v);
X    a = (union allocNode *) q;
X    a->link = freeNum;
X    freeNum = a;
X}
X
X/* END CODE */
END_OF_FILE
if test 20144 -ne `wc -c <'qmath.c'`; then
    echo shar: \"'qmath.c'\" unpacked with wrong size!
fi
# end of 'qmath.c'
fi
if test -f 'qtrans.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'qtrans.c'\"
else
echo shar: Extracting \"'qtrans.c'\" \(20565 characters\)
sed "s/^X//" >'qtrans.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 * Transcendental functions for real numbers.
X * These are sin, cos, exp, ln, power, cosh, sinh.
X */
X
X#include "math.h"
X
XBOOL _sinisneg_;    /* whether sin(x) < 0 (set by cos(x)) */
X
X
X/*
X * Calculate the cosine of a number with an accuracy within epsilon.
X * This also saves the sign of the corresponding sin function.
X */
XNUMBER *
Xqcos(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *term, *sum, *qsq, *epsilon2, *tmp;
X    FULL n, i;
X    long scale, bits, bits2;
X
X    _sinisneg_ = qisneg(q);
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for cosine");
X    if (qiszero(q))
X        return qlink(&_qone_);
X    bits = qprecision(epsilon) + 1;
X    epsilon = qscale(epsilon, -4L);
X    /*
X     * If the argument is larger than one, then divide it by a power of two
X     * so that it is one or less.  This will make the series converge quickly.
X     * We will extrapolate the result for the original argument afterwards.
X     */
X    scale = zhighbit(q->num) - zhighbit(q->den) + 1;
X    if (scale < 0)
X        scale = 0;
X    if (scale > 0) {
X        q = qscale(q, -scale);
X        tmp = qscale(epsilon, -scale);
X        qfree(epsilon);
X        epsilon = tmp;
X    }
X    epsilon2 = qscale(epsilon, -4L);
X    qfree(epsilon);
X    bits2 = qprecision(epsilon2) + 10;
X    /*
X     * Now use the Taylor series expansion to calculate the cosine.
X     * Keep using approximations so that the fractions don't get too large.
X     */
X    qsq = qsquare(q);
X    if (scale > 0)
X        qfree(q);
X    term = qlink(&_qone_);
X    sum = qlink(&_qone_);
X    n = 0;
X    while (qrel(term, epsilon2) > 0) {
X        i = ++n;
X        i *= ++n;
X        tmp = qmul(term, qsq);
X        qfree(term);
X        term = qdivi(tmp, (long) i);
X        qfree(tmp);
X        tmp = qbround(term, bits2);
X        qfree(term);
X        term = tmp;
X        if (n & 2)
X            tmp = qsub(sum, term);
X        else
X            tmp = qadd(sum, term);
X        qfree(sum);
X        sum = qbround(tmp, bits2);
X        qfree(tmp);
X    }
X    qfree(term);
X    qfree(qsq);
X    qfree(epsilon2);
X    /*
X     * Now scale back up to the original value of x by using the formula:
X     *    cos(2 * x) = 2 * (cos(x) ^ 2) - 1.
X     */
X    while (--scale >= 0) {
X        if (qisneg(sum))
X            _sinisneg_ = !_sinisneg_;
X        tmp = qsquare(sum);
X        qfree(sum);
X        sum = qscale(tmp, 1L);
X        qfree(tmp);
X        tmp = qdec(sum);
X        qfree(sum);
X        sum = qbround(tmp, bits2);
X        qfree(tmp);
X    }
X    tmp = qbround(sum, bits);
X    qfree(sum);
X    return tmp;
X}
X
X
X/*
X * Calculate the sine of a number with an accuracy within epsilon.
X * This is calculated using the formula:
X *    sin(x)^2 + cos(x)^2 = 1.
X * The only tricky bit is resolving the sign of the result.
X * Future: Use sin(3*x) = 3*sin(x) - 4*sin(x)^3.
X */
XNUMBER *
Xqsin(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *epsilon2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for sine");
X    if (qiszero(q))
X        return qlink(q);
X    epsilon2 = qscale(epsilon, -4L);
X    tmp1 = qcos(q, epsilon2);
X    qfree(epsilon2);
X    tmp2 = qlegtoleg(tmp1, epsilon, _sinisneg_);
X    qfree(tmp1);
X    return tmp2;
X}
X
X
X/*
X * Calculate the tangent function.
X */
XNUMBER *
Xqtan(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *cosval, *sinval, *epsilon2, *tmp, *res;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for tangent");
X    if (qiszero(q))
X        return qlink(q);
X    epsilon2 = qscale(epsilon, -4L);
X    cosval = qcos(q, epsilon2);
X    sinval = qlegtoleg(cosval, epsilon2, _sinisneg_);
X    qfree(epsilon2);
X    tmp = qdiv(sinval, cosval);
X    qfree(cosval);
X    qfree(sinval);
X    res = qbround(tmp, qprecision(epsilon) + 1);
X    qfree(tmp);
X    return res;
X}
X
X
X/*
X * Calculate the arcsine function.
X * The result is in the range -pi/2 to pi/2.
X */
XNUMBER *
Xqasin(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *sum, *term, *epsilon2, *qsq, *tmp;
X    FULL n, i;
X    long bits, bits2;
X    int neg;
X    NUMBER mulnum;
X    HALF numval[2];
X    HALF denval[2];
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for arcsine");
X    if (qiszero(q))
X        return qlink(&_qzero_);
X    if ((qrel(q, &_qone_) > 0) || (qrel(q, &_qnegone_) < 0))
X        error("Argument too large for asin");
X    neg = qisneg(q);
X    q = qabs(q);
X    epsilon = qscale(epsilon, -4L);
X    epsilon2 = qscale(epsilon, -4L);
X    mulnum.num.sign = 0;
X    mulnum.num.len = 1;
X    mulnum.num.v = numval;
X    mulnum.den.sign = 0;
X    mulnum.den.len = 1;
X    mulnum.den.v = denval;
X    /*
X     * If the argument is too near one (we use .5) then reduce the
X     * argument to a more accurate range using the formula:
X     *    asin(x) = 2 * asin(sqrt((1 - sqrt(1 - x^2)) / 2)).
X     */
X    if (qrel(q, &_qonehalf_) > 0) {
X        sum = qlegtoleg(q, epsilon2, FALSE);
X        qfree(q);
X        tmp = qsub(&_qone_, sum);
X        qfree(sum);
X        sum = qscale(tmp, -1L);
X        qfree(tmp);
X        tmp = qsqrt(sum, epsilon2);
X        qfree(sum);
X        qfree(epsilon2);
X        sum = qasin(tmp, epsilon);
X        qfree(tmp);
X        qfree(epsilon);
X        tmp = qscale(sum, 1L);
X        qfree(sum);
X        if (neg) {
X            sum = qneg(tmp);
X            qfree(tmp);
X            tmp = sum;
X        }
X        return tmp;
X    }
X    /*
X     * Argument is between zero and .5, so use the series.
X     */
X    epsilon = qscale(epsilon, -4L);
X    epsilon2 = qscale(epsilon, -4L);
X    bits = qprecision(epsilon) + 1;
X    bits2 = bits + 10;
X    sum = qlink(q);
X    term = qlink(q);
X    qsq = qsquare(q);
X    qfree(q);
X    n = 1;
X    while (qrel(term, epsilon2) > 0) {
X        i = n * n;
X        numval[0] = i & BASE1;
X        if (i >= BASE) {
X            numval[1] = i / BASE;
X            mulnum.den.len = 2;
X        }
X        i = (n + 1) * (n + 2);
X        denval[0] = i & BASE1;
X        if (i >= BASE) {
X            denval[1] = i / BASE;
X            mulnum.den.len = 2;
X        }
X        tmp = qmul(term, qsq);
X        qfree(term);
X        term = qmul(tmp, &mulnum);
X        qfree(tmp);
X        tmp = qbround(term, bits2);
X        qfree(term);
X        term = tmp;
X        tmp = qadd(sum, term);
X        qfree(sum);
X        sum = qbround(tmp, bits2);
X        qfree(tmp);
X        n += 2;
X    }
X    qfree(epsilon);
X    qfree(epsilon2);
X    qfree(term);
X    qfree(qsq);
X    tmp = qbround(sum, bits);
X    qfree(sum);
X    if (neg) {
X        term = qneg(tmp);
X        qfree(tmp);
X        tmp = term;
X    }
X    return tmp;
X}
X
X
X/*
X * Calculate the acos function.
X * The result is in the range 0 to pi.
X */
XNUMBER *
Xqacos(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *epsilon2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for arccosine");
X    if (qisone(q))
X        return qlink(&_qzero_);
X    if ((qrel(q, &_qone_) > 0) || (qrel(q, &_qnegone_) < 0))
X        error("Argument too large for acos");
X    /*
X     * Calculate the result using the formula:
X     *    acos(x) = asin(sqrt(1 - x^2)).
X     */
X    epsilon2 = qscale(epsilon, -8L);
X    tmp1 = qlegtoleg(q, epsilon2, FALSE);
X    qfree(epsilon2);
X    tmp2 = qasin(tmp1, epsilon);
X    qfree(tmp1);
X    return tmp2;
X}
X
X
X/*
X * Calculate the arctangent function with a accuracy less than epsilon.
X * This uses the formula:
X *    atan(x) = asin(sqrt(x^2 / (x^2 + 1))).
X */
XNUMBER *
Xqatan(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for arctangent");
X    if (qiszero(q))
X        return qlink(&_qzero_);
X    tmp1 = qsquare(q);
X    tmp2 = qinc(tmp1);
X    tmp3 = qdiv(tmp1, tmp2);
X    qfree(tmp1);
X    qfree(tmp2);
X    epsilon2 = qscale(epsilon, -8L);
X    tmp1 = qsqrt(tmp3, epsilon2);
X    qfree(epsilon2);
X    qfree(tmp3);
X    tmp2 = qasin(tmp1, epsilon);
X    qfree(tmp1);
X    if (qisneg(q)) {
X        tmp1 = qneg(tmp2);
X        qfree(tmp2);
X        tmp2 = tmp1;
X    }
X    return tmp2;
X}
X
X
X/*
X * Calculate the angle which is determined by the point (x,y).
X * This is the same as arctan for non-negative x, but gives the correct
X * value for negative x.  By convention, y is the first argument.
X * For example, qatan2(1, -1) = 3/4 * pi.
X */
XNUMBER *
Xqatan2(qy, qx, epsilon)
X    NUMBER *qy, *qx, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *epsilon2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for atan2");
X    if (qiszero(qy) && qiszero(qx))
X        error("Zero coordinates for atan2");
X    /*
X     * If the point is on the negative real axis, then the answer is pi.
X     */
X    if (qiszero(qy) && qisneg(qx))
X        return qpi(epsilon);
X    /*
X     * If the point is in the right half plane, then use the normal atan.
X     */
X    if (!qisneg(qx) && !qiszero(qx)) {
X        if (qiszero(qy))
X            return qlink(&_qzero_);
X        tmp1 = qdiv(qy, qx);
X        tmp2 = qatan(tmp1, epsilon);
X        qfree(tmp1);
X        return tmp2;
X    }
X    /*
X     * The point is in the left half plane.  Calculate the angle by finding
X     * the atan of half the angle using the formula:
X     *    atan2(y,x) = 2 * atan((sqrt(x^2 + y^2) - x) / y).
X     */
X    epsilon2 = qscale(epsilon, -4L);
X    tmp1 = qhypot(qx, qy, epsilon2);
X    tmp2 = qsub(tmp1, qx);
X    qfree(tmp1);
X    tmp1 = qdiv(tmp2, qy);
X    qfree(tmp2);
X    tmp2 = qatan(tmp1, epsilon2);
X    qfree(tmp1);
X    qfree(epsilon2);
X    tmp1 = qscale(tmp2, 1L);
X    qfree(tmp2);
X    return tmp1;
X}
X
X
X/*
X * Calculate the value of pi to within the required epsilon.
X * This uses the following formula which only needs integer calculations
X * except for the final operation:
X *    pi = 1 / SUMOF(comb(2 * N, N) ^ 3 * (42 * N + 5) / 2 ^ (12 * N + 4)),
X * where the summation runs from N=0.  This formula gives about 6 bits of
X * accuracy per term.  Since the denominator for each term is a power of two,
X * we can simply use shifts to sum the terms.  The combinatorial numbers
X * in the formula are calculated recursively using the formula:
X *    comb(2*(N+1), N+1) = 2 * comb(2 * N, N) * (2 * N + 1) / N.
X */
XNUMBER *
Xqpi(epsilon)
X    NUMBER *epsilon;
X{
X    ZVALUE comb;            /* current combinatorial value */
X    ZVALUE sum;            /* current sum */
X    ZVALUE tmp1, tmp2;
X    NUMBER *r, *t1, qtmp;
X    long shift;            /* current shift of result */
X    long N;                /* current term number */
X    long bits;            /* needed number of bits of precision */
X    long t;
X
X    if (qiszero(epsilon) || qisneg(epsilon))
X        error("Bad epsilon value for pi");
X    bits = qprecision(epsilon) + 4;
X    comb = _one_;
X    itoz(5L, &sum);
X    N = 0;
X    shift = 4;
X    do {
X        t = 1 + (++N & 0x1);
X        (void) zdivi(comb, N / (3 - t), &tmp1);
X        freeh(comb.v);
X        zmuli(tmp1, t * (2 * N - 1), &comb);
X        freeh(tmp1.v);
X        zsquare(comb, &tmp1);
X        zmul(comb, tmp1, &tmp2);
X        freeh(tmp1.v);
X        zmuli(tmp2, 42 * N + 5, &tmp1);
X        freeh(tmp2.v);
X        zshift(sum, 12L, &tmp2);
X        freeh(sum.v);
X        zadd(tmp1, tmp2, &sum);
X        t = zhighbit(tmp1);
X        freeh(tmp1.v);
X        freeh(tmp2.v);
X        shift += 12;
X    } while ((shift - t) < bits);
X    qtmp.num = _one_;
X    qtmp.den = sum;
X    t1 = qscale(&qtmp, shift);
X    freeh(sum.v);
X    r = qbround(t1, bits);
X    qfree(t1);
X    return r;
X}
X
X
X/*
X * Calculate the exponential function with a relative accuracy less than
X * epsilon.
X */
XNUMBER *
Xqexp(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    long scale;
X    FULL n;
X    long bits, bits2;
X    NUMBER *sum, *term, *qs, *epsilon2, *tmp;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for exp");
X    if (qiszero(q))
X        return qlink(&_qone_);
X    epsilon = qscale(epsilon, -4L);
X    /*
X     * If the argument is larger than one, then divide it by a power of two
X     * so that it is one or less.  This will make the series converge quickly.
X     * We will extrapolate the result for the original argument afterwards.
X     * Also make the argument non-negative.
X     */
X    qs = qabs(q);
X    scale = zhighbit(q->num) - zhighbit(q->den) + 1;
X    if (scale < 0)
X        scale = 0;
X    if (scale > 0) {
X        if (scale > 100000)
X            error("Very large argument for exp");
X        tmp = qscale(qs, -scale);
X        qfree(qs);
X        qs = tmp;
X        tmp = qscale(epsilon, -scale);
X        qfree(epsilon);
X        epsilon = tmp;
X    }
X    epsilon2 = qscale(epsilon, -4L);
X    bits = qprecision(epsilon) + 1;
X    bits2 = bits + 10;
X    qfree(epsilon);
X    /*
X     * Now use the Taylor series expansion to calculate the exponential.
X     * Keep using approximations so that the fractions don't get too large.
X     */
X    sum = qlink(&_qone_);
X    term = qlink(&_qone_);
X    n = 0;
X    while (qrel(term, epsilon2) > 0) {
X        n++;
X        tmp = qmul(term, qs);
X        qfree(term);
X        term = qdivi(tmp, (long) n);
X        qfree(tmp);
X        tmp = qbround(term, bits2);
X        qfree(term);
X        term = tmp;
X        tmp = qadd(sum, term);
X        qfree(sum);
X        sum = qbround(tmp, bits2);
X        qfree(tmp);
X    }
X    qfree(term);
X    qfree(qs);
X    qfree(epsilon2);
X    /*
X     * Now repeatedly square the answer to get the final result.
X     * Then invert it if the original argument was negative.
X     */
X    while (--scale >= 0) {
X        tmp = qsquare(sum);
X        qfree(sum);
X        sum = qbround(tmp, bits2);
X        qfree(tmp);
X    }
X    tmp = qbround(sum, bits);
X    qfree(sum);
X    if (qisneg(q)) {
X        sum = qinv(tmp);
X        qfree(tmp);
X        tmp = sum;
X    }
X    return tmp;
X}
X
X
X/*
X * Calculate the natural logarithm of a number accurate to the specified
X * epsilon.
X */
XNUMBER *
Xqln(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *term, *term2, *sum, *epsilon2, *tmp1, *tmp2;
X    NUMBER *minr, *maxr;
X    long shift, bits, bits2;
X    int neg;
X    FULL n;
X
X    if (qisneg(q) || qiszero(q))
X        error("log of non-positive number");
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon for ln");
X    epsilon = qscale(epsilon, -4L);
X    epsilon2 = qscale(epsilon, -4L);
X    bits = qprecision(epsilon) + 1;
X    bits2 = bits + 10;
X    qfree(epsilon);
X    /*
X     * Scale down the logarithm to the convergent range by taking square
X     * roots.  The result will be modified to get the final value later.
X     */
X    q = qlink(q);
X    minr = iitoq(7L, 10L);
X    maxr = iitoq(7L, 5L);
X    shift = 1;
X    while ((qrel(q, minr) < 0) || (qrel(q, maxr) > 0)) {
X        tmp1 = qsqrt(q, epsilon2);
X        qfree(q);
X        q = tmp1;
X        shift++;
X    }
X    qfree(minr);
X    qfree(maxr);
X    /*
X     * Calculate a value which will always converge using the formula:
X     *    ln((1+x)/(1-x)) = ln(1+x) - ln(1-x).
X     */
X    tmp1 = qdec(q);
X    tmp2 = qinc(q);
X    term = qdiv(tmp1, tmp2);
X    qfree(tmp1);
X    qfree(tmp2);
X    qfree(q);
X    neg = 0;
X    if (qisneg(term)) {
X        neg = 1;
X        tmp1 = qabs(term);
X        qfree(term);
X        term = tmp1;
X    }
X    /*
X     * Now use the Taylor series expansion to calculate the result.
X     */
X    n = 1;
X    term2 = qsquare(term);
X    sum = qlink(term);
X    while (qrel(term, epsilon2) > 0) {
X        n += 2;
X        tmp1 = qmul(term, term2);
X        qfree(term);
X        term = qbround(tmp1, bits2);
X        qfree(tmp1);
X        tmp1 = qdivi(term, (long) n);
X        tmp2 = qadd(sum, tmp1);
X        qfree(tmp1);
X        qfree(sum);
X        sum = qbround(tmp2, bits2);
X    }
X    qfree(epsilon2);
X    qfree(term);
X    qfree(term2);
X    if (neg) {
X        tmp1 = qneg(sum);
X        qfree(sum);
X        sum = tmp1;
X    }
X    /*
X     * Calculate the final result by multiplying by the proper power of two
X     * to undo the square roots done at the top.
X     */
X    tmp1 = qscale(sum, shift);
X    qfree(sum);
X    sum = qbround(tmp1, bits);
X    qfree(tmp1);
X    return sum;
X}
X
X
X/*
X * Calculate the result of raising one number to the power of another.
X * The result is calculated to within the specified relative error.
X */
XNUMBER *
Xqpower(q1, q2, epsilon)
X    NUMBER *q1, *q2, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *epsilon2;
X
X    if (qisint(q2))
X        return qpowi(q1, q2);
X    epsilon2 = qscale(epsilon, -4L);
X    tmp1 = qln(q1, epsilon2);
X    tmp2 = qmul(tmp1, q2);
X    qfree(tmp1);
X    tmp1 = qexp(tmp2, epsilon);
X    qfree(tmp2);
X    qfree(epsilon2);
X    return tmp1;
X}
X
X
X/*
X * Calculate the Kth root of a number to within the specified accuracy.
X */
XNUMBER *
Xqroot(q1, q2, epsilon)
X    NUMBER *q1, *q2, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *epsilon2;
X    int neg;
X
X    if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
X        error("Taking bad root of number");
X    if (qiszero(q1) || qisone(q1) || qisone(q2))
X        return qlink(q1);
X    if (qistwo(q2))
X        return qsqrt(q1, epsilon);
X    neg = qisneg(q1);
X    if (neg) {
X        if (iseven(q2->num))
X            error("Taking even root of negative number");
X        q1 = qabs(q1);
X    }
X    epsilon2 = qscale(epsilon, -4L);
X    tmp1 = qln(q1, epsilon2);
X    tmp2 = qdiv(tmp1, q2);
X    qfree(tmp1);
X    tmp1 = qexp(tmp2, epsilon);
X    qfree(tmp2);
X    qfree(epsilon2);
X    if (neg) {
X        tmp2 = qneg(tmp1);
X        qfree(tmp1);
X        tmp1 = tmp2;
X    }
X    return tmp1;
X}
X
X
X/*
X * Calculate the hyperbolic cosine function with a relative accuracy less
X * than epsilon.  This is defined by:
X *    cosh(x) = (exp(x) + exp(-x)) / 2.
X */
XNUMBER *
Xqcosh(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    long scale;
X    FULL n;
X    FULL m;
X    long bits, bits2;
X    NUMBER *sum, *term, *qs, *epsilon2, *tmp;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for exp");
X    if (qiszero(q))
X        return qlink(&_qone_);
X    epsilon = qscale(epsilon, -4L);
X    /*
X     * If the argument is larger than one, then divide it by a power of two
X     * so that it is one or less.  This will make the series converge quickly.
X     * We will extrapolate the result for the original argument afterwards.
X     */
X    qs = qabs(q);
X    scale = zhighbit(q->num) - zhighbit(q->den) + 1;
X    if (scale < 0)
X        scale = 0;
X    if (scale > 0) {
X        if (scale > 100000)
X            error("Very large argument for exp");
X        tmp = qscale(qs, -scale);
X        qfree(qs);
X        qs = tmp;
X        tmp = qscale(epsilon, -scale);
X        qfree(epsilon);
X        epsilon = tmp;
X    }
X    epsilon2 = qscale(epsilon, -4L);
X    bits = qprecision(epsilon) + 1;
X    bits2 = bits + 10;
X    qfree(epsilon);
X    tmp = qsquare(qs);
X    qfree(qs);
X    qs = tmp;
X    /*
X     * Now use the Taylor series expansion to calculate the exponential.
X     * Keep using approximations so that the fractions don't get too large.
X     */
X    sum = qlink(&_qone_);
X    term = qlink(&_qone_);
X    n = 0;
X    while (qrel(term, epsilon2) > 0) {
X        m = ++n;
X        m *= ++n;
X        tmp = qmul(term, qs);
X        qfree(term);
X        term = qdivi(tmp, (long) m);
X        qfree(tmp);
X        tmp = qbround(term, bits2);
X        qfree(term);
X        term = tmp;
X        tmp = qadd(sum, term);
X        qfree(sum);
X        sum = qbround(tmp, bits2);
X        qfree(tmp);
X    }
X    qfree(term);
X    qfree(qs);
X    qfree(epsilon2);
X    /*
X     * Now bring the number back up into range to get the final result.
X     * This uses the formula:
X     *    cosh(2 * x) = 2 * cosh(x)^2 - 1.
X     */
X    while (--scale >= 0) {
X        tmp = qsquare(sum);
X        qfree(sum);
X        sum = qscale(tmp, 1L);
X        qfree(tmp);
X        tmp = qdec(sum);
X        qfree(sum);
X        sum = qbround(tmp, bits2);
X        qfree(tmp);
X    }
X    tmp = qbround(sum, bits);
X    qfree(sum);
X    return tmp;
X}
X
X
X/*
X * Calculate the hyperbolic sine with an accurary less than epsilon.
X * This is calculated using the formula:
X *    cosh(x)^2 - sinh(x)^2 = 1.
X */
XNUMBER *
Xqsinh(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *tmp1, *tmp2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for sinh");
X    if (qiszero(q))
X        return qlink(q);
X    epsilon = qscale(epsilon, -4L);
X    tmp1 = qcosh(q, epsilon);
X    tmp2 = qsquare(tmp1);
X    qfree(tmp1);
X    tmp1 = qdec(tmp2);
X    qfree(tmp2);
X    tmp2 = qsqrt(tmp1, epsilon);
X    qfree(tmp1);
X    if (qisneg(q)) {
X        tmp1 = qneg(tmp2);
X        qfree(tmp2);
X        tmp2 = tmp1;
X    }
X    qfree(epsilon);
X    return tmp2;
X}
X
X
X/*
X * Calculate the hyperbolic tangent with an accurary less than epsilon.
X * This is calculated using the formula:
X *    tanh(x) = sinh(x) / cosh(x).
X */
XNUMBER *
Xqtanh(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *coshval;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for tanh");
X    if (qiszero(q))
X        return qlink(q);
X    epsilon = qscale(epsilon, -4L);
X    coshval = qcosh(q, epsilon);
X    tmp2 = qsquare(coshval);
X    tmp1 = qdec(tmp2);
X    qfree(tmp2);
X    tmp2 = qsqrt(tmp1, epsilon);
X    qfree(tmp1);
X    if (qisneg(q)) {
X        tmp1 = qneg(tmp2);
X        qfree(tmp2);
X        tmp2 = tmp1;
X    }
X    qfree(epsilon);
X    tmp1 = qdiv(tmp2, coshval);
X    qfree(tmp2);
X    qfree(coshval);
X    return tmp1;
X}
X
X
X/*
X * Compute the hyperbolic arccosine within the specified accuracy.
X * This is calculated using the formula:
X *    acosh(x) = ln(x + sqrt(x^2 - 1)).
X */
XNUMBER *
Xqacosh(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *epsilon2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for acosh");
X    if (qisone(q))
X        return qlink(&_qzero_);
X    if (qreli(q, 1L) < 0)
X        error("Argument less than one for acosh");
X    epsilon2 = qscale(epsilon, -8L);
X    tmp1 = qsquare(q);
X    tmp2 = qdec(tmp1);
X    qfree(tmp1);
X    tmp1 = qsqrt(tmp2, epsilon2);
X    qfree(tmp2);
X    tmp2 = qadd(tmp1, q);
X    qfree(tmp1);
X    tmp1 = qln(tmp2, epsilon);
X    qfree(tmp2);
X    qfree(epsilon2);
X    return tmp1;
X}
X
X
X/*
X * Compute the hyperbolic arcsine within the specified accuracy.
X * This is calculated using the formula:
X *    asinh(x) = ln(x + sqrt(x^2 + 1)).
X */
XNUMBER *
Xqasinh(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *epsilon2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for asinh");
X    if (qiszero(q))
X        return qlink(&_qzero_);
X    epsilon2 = qscale(epsilon, -8L);
X    tmp1 = qsquare(q);
X    tmp2 = qinc(tmp1);
X    qfree(tmp1);
X    tmp1 = qsqrt(tmp2, epsilon2);
X    qfree(tmp2);
X    tmp2 = qadd(tmp1, q);
X    qfree(tmp1);
X    tmp1 = qln(tmp2, epsilon);
X    qfree(tmp2);
X    qfree(epsilon2);
X    return tmp1;
X}
X
X
X/*
X * Compute the hyperbolic arctangent within the specified accuracy.
X * This is calculated using the formula:
X *    atanh(x) = ln((1 + u) / (1 - u)) / 2.
X */
XNUMBER *
Xqatanh(q, epsilon)
X    NUMBER *q, *epsilon;
X{
X    NUMBER *tmp1, *tmp2, *tmp3;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        error("Illegal epsilon value for atanh");
X    if (qiszero(q))
X        return qlink(&_qzero_);
X    if ((qreli(q, 1L) > 0) || (qreli(q, -1L) < 0))
X        error("Argument not between -1 and 1 for atanh");
X    tmp1 = qinc(q);
X    tmp2 = qsub(&_qone_, q);
X    tmp3 = qdiv(tmp1, tmp2);
X    qfree(tmp1);
X    qfree(tmp2);
X    tmp1 = qln(tmp3, epsilon);
X    qfree(tmp3);
X    tmp2 = qscale(tmp1, -1L);
X    qfree(tmp1);
X    return tmp2;
X}
X
X/* END CODE */
END_OF_FILE
if test 20565 -ne `wc -c <'qtrans.c'`; then
    echo shar: \"'qtrans.c'\" unpacked with wrong size!
fi
# end of 'qtrans.c'
fi
echo shar: End of archive 9 \(of 21\).
cp /dev/null ark9isdone
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