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

Text File | 1992-05-09 | 56.8 KB | 2,761 lines

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

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

#! /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 5 (of 21)."
# Contents:  comfunc.c commath.c file.c listfunc.c qmod.c
# Wrapped by dbell@elm on Tue Feb 25 15:20:59 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'comfunc.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'comfunc.c'\"
else
echo shar: Extracting \"'comfunc.c'\" \(10332 characters\)
sed "s/^X//" >'comfunc.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 complex arithmetic non-primitive routines
X */
X
X#include "calc.h"
X
X
X/*
X * Round a complex number to the specified number of decimal places.
X * This simply means to round each of the components of the number.
X * Zero decimal places means round to the nearest complex integer.
X */
XCOMPLEX *
Xcround(c, places)
X    COMPLEX *c;
X    long places;
X{
X    COMPLEX *res;        /* result */
X
X    res = comalloc();
X    res->real = qround(c->real, places);
X    res->imag = qround(c->imag, places);
X    return res;
X}
X
X
X/*
X * Round a complex number to the specified number of binary decimal places.
X * This simply means to round each of the components of the number.
X * Zero binary places means round to the nearest complex integer.
X */
XCOMPLEX *
Xcbround(c, places)
X    COMPLEX *c;
X    long places;
X{
X    COMPLEX *res;        /* result */
X
X    res = comalloc();
X    res->real = qbround(c->real, places);
X    res->imag = qbround(c->imag, places);
X    return res;
X}
X
X
X/*
X * Compute the result of raising a complex number to an integer power.
X */
XCOMPLEX *
Xcpowi(c, q)
X    COMPLEX *c;        /* complex number to be raised */
X    NUMBER *q;        /* power to raise it to */
X{
X    COMPLEX *tmp, *res;    /* temporary values */
X    long power;        /* power to raise to */
X    unsigned long bit;    /* current bit value */
X    int sign;
X
X    if (qisfrac(q))
X        error("Raising number to non-integral power");
X    if (isbig(q->num))
X        error("Raising number to very large power");
X    power = (istiny(q->num) ? z1tol(q->num) : z2tol(q->num));
X    if (ciszero(c) && (power == 0))
X        error("Raising zero to zeroth power");
X    sign = 1;
X    if (qisneg(q))
X        sign = -1;
X    /*
X     * Handle some low powers specially
X     */
X    if (power <= 4) {
X        switch ((int) (power * sign)) {
X            case 0:
X                return clink(&_cone_);
X            case 1:
X                return clink(c);
X            case -1:
X                return cinv(c);
X            case 2:
X                return csquare(c);
X            case -2:
X                tmp = csquare(c);
X                res = cinv(tmp);
X                comfree(tmp);
X                return res;
X            case 3:
X                tmp = csquare(c);
X                res = cmul(c, tmp);
X                comfree(tmp);
X                return res;
X            case 4:
X                tmp = csquare(c);
X                res = csquare(tmp);
X                comfree(tmp);
X                return res;
X        }
X    }
X    /*
X     * Compute the power by squaring and multiplying.
X     * This uses the left to right method of power raising.
X     */
X    bit = TOPFULL;
X    while ((bit & power) == 0)
X        bit >>= 1L;
X    bit >>= 1L;
X    res = csquare(c);
X    if (bit & power) {
X        tmp = cmul(res, c);
X        comfree(res);
X        res = tmp;
X    }
X    bit >>= 1L;
X    while (bit) {
X        tmp = csquare(res);
X        comfree(res);
X        res = tmp;
X        if (bit & power) {
X            tmp = cmul(res, c);
X            comfree(res);
X            res = tmp;
X        }
X        bit >>= 1L;
X    }
X    if (sign < 0) {
X        tmp = cinv(res);
X        comfree(res);
X        res = tmp;
X    }
X    return res;
X}
X
X
X/*
X * Calculate the square root of a complex number, with each component
X * within the specified error. This uses the formula:
X *    sqrt(a+bi) = sqrt((a+sqrt(a^2+b^2))/2) + sqrt((-a+sqrt(a^2+b^2))/2)i.
X */
XCOMPLEX *
Xcsqrt(c, epsilon)
X    COMPLEX *c;
X    NUMBER *epsilon;
X{
X    COMPLEX *r;
X    NUMBER *hypot, *tmp1, *tmp2;
X
X    if (ciszero(c) || cisone(c))
X        return clink(c);
X    r = comalloc();
X    if (cisreal(c)) {
X        if (!qisneg(c->real)) {
X            r->real = qsqrt(c->real, epsilon);
X            return r;
X        }
X        tmp1 = qneg(c->real);
X        r->imag = qsqrt(tmp1, epsilon);
X        qfree(tmp1);
X        return r;
X    }
X    hypot = qhypot(c->real, c->imag, epsilon);
X    tmp1 = qadd(hypot, c->real);
X    tmp2 = qscale(tmp1, -1L);
X    qfree(tmp1);
X    r->real = qsqrt(tmp2, epsilon);
X    qfree(tmp2);
X    tmp1 = qsub(hypot, c->real);
X    qfree(hypot);
X    tmp2 = qscale(tmp1, -1L);
X    qfree(tmp1);
X    tmp1 = qsqrt(tmp2, epsilon);
X    qfree(tmp2);
X    if (qisneg(c->imag)) {
X        tmp2 = qneg(tmp1);
X        qfree(tmp1);
X        tmp1 = tmp2;
X    }
X    r->imag = tmp1;
X    return r;
X}
X
X
X/*
X * Take the Nth root of a complex number, where N is a positive integer.
X * Each component of the result is within the specified error.
X */
XCOMPLEX *
Xcroot(c, q, epsilon)
X    COMPLEX *c;
X    NUMBER *q, *epsilon;
X{
X    COMPLEX *r;
X    NUMBER *a2pb2, *root, *tmp1, *tmp2, *epsilon2;
X
X    if (qisneg(q) || qiszero(q) || qisfrac(q))
X        error("Taking bad root of complex number");
X    if (cisone(c) || qisone(q))
X        return clink(c);
X    if (qistwo(q))
X        return csqrt(c, epsilon);
X    r = comalloc();
X    if (cisreal(c) && !qisneg(c->real)) {
X        r->real = qroot(c->real, q, epsilon);
X        return r;
X    }
X    /*
X     * Calculate the root using the formula:
X     *    croot(a + bi, n) =
X     *        cpolar(qroot(a^2 + b^2, 2 * n), qatan2(b, a) / n).
X     */
X    epsilon2 = qscale(epsilon, -8L);
X    tmp1 = qsquare(c->real);
X    tmp2 = qsquare(c->imag);
X    a2pb2 = qadd(tmp1, tmp2);
X    qfree(tmp1);
X    qfree(tmp2);
X    tmp1 = qscale(q, 1L);
X    root = qroot(a2pb2, tmp1, epsilon2);
X    qfree(a2pb2);
X    qfree(tmp1);
X    tmp1 = qatan2(c->imag, c->real, epsilon2);
X    qfree(epsilon2);
X    tmp2 = qdiv(tmp1, q);
X    qfree(tmp1);
X    r = cpolar(root, tmp2, epsilon);
X    qfree(root);
X    qfree(tmp2);
X    return r;
X}
X
X
X/*
X * Calculate the complex exponential function to the desired accuracy.
X * We use the formula:
X *    exp(a + bi) = exp(a) * (cos(b) + i * sin(b)).
X */
XCOMPLEX *
Xcexp(c, epsilon)
X    COMPLEX *c;
X    NUMBER *epsilon;
X{
X    COMPLEX *r;
X    NUMBER *tmp1, *tmp2, *epsilon2;
X
X    if (ciszero(c))
X        return clink(&_cone_);
X    r = comalloc();
X    if (cisreal(c)) {
X        r->real = qexp(c->real, epsilon);
X        return r;
X    }
X    epsilon2 = qscale(epsilon, -2L);
X    r->real = qcos(c->imag, epsilon2);
X    r->imag = qlegtoleg(r->real, epsilon2, _sinisneg_);
X    if (qiszero(c->real)) {
X        qfree(epsilon2);
X        return r;
X    }
X    tmp1 = qexp(c->real, epsilon2);
X    qfree(epsilon2);
X    tmp2 = qmul(r->real, tmp1);
X    qfree(r->real);
X    r->real = tmp2;
X    tmp2 = qmul(r->imag, tmp1);
X    qfree(r->imag);
X    qfree(tmp1);
X    r->imag = tmp2;
X    return r;
X}
X
X
X/*
X * Calculate the natural logarithm of a complex number within the specified
X * error.  We use the formula:
X *    ln(a + bi) = ln(a^2 + b^2) / 2 + i * atan2(b, a).
X */
XCOMPLEX *
Xcln(c, epsilon)
X    COMPLEX *c;
X    NUMBER *epsilon;
X{
X    COMPLEX *r;
X    NUMBER *a2b2, *tmp1, *tmp2;
X
X    if (ciszero(c))
X        error("Logarithm of zero");
X    if (cisone(c))
X        return clink(&_czero_);
X    r = comalloc();
X    if (cisreal(c) && !qisneg(c->real)) {
X        r->real = qln(c->real, epsilon);
X        return r;
X    }
X    tmp1 = qsquare(c->real);
X    tmp2 = qsquare(c->imag);
X    a2b2 = qadd(tmp1, tmp2);
X    qfree(tmp1);
X    qfree(tmp2);
X    tmp1 = qln(a2b2, epsilon);
X    qfree(a2b2);
X    r->real = qscale(tmp1, -1L);
X    qfree(tmp1);
X    r->imag = qatan2(c->imag, c->real, epsilon);
X    return r;
X}
X
X
X/*
X * Calculate the complex cosine within the specified accuracy.
X * This uses the formula:
X *    cos(a + bi) = cos(a) * cosh(b) - sin(a) * sinh(b) * i.
X */
XCOMPLEX *
Xccos(c, epsilon)
X    COMPLEX *c;
X    NUMBER *epsilon;
X{
X    COMPLEX *r;
X    NUMBER *cosval, *coshval, *tmp1, *tmp2, *tmp3, *epsilon2;
X    int negimag;
X
X    if (ciszero(c))
X        return clink(&_cone_);
X    r = comalloc();
X    if (cisreal(c)) {
X        r->real = qcos(c->real, epsilon);
X        return r;
X    }
X    if (qiszero(c->real)) {
X        r->real = qcosh(c->imag, epsilon);
X        return r;
X    }
X    epsilon2 = qscale(epsilon, -2L);
X    coshval = qcosh(c->imag, epsilon2);
X    cosval = qcos(c->real, epsilon2);
X    negimag = !_sinisneg_;
X    if (qisneg(c->imag))
X        negimag = !negimag;
X    r->real = qmul(cosval, coshval);
X    /*
X     * Calculate the imaginary part using the formula:
X     *    sin(a) * sinh(b) = sqrt((1 - a^2) * (b^2 - 1)).
X     */
X    tmp1 = qsquare(cosval);
X    qfree(cosval);
X    tmp2 = qdec(tmp1);
X    qfree(tmp1);
X    tmp1 = qneg(tmp2);
X    qfree(tmp2);
X    tmp2 = qsquare(coshval);
X    qfree(coshval);
X    tmp3 = qdec(tmp2);
X    qfree(tmp2);
X    tmp2 = qmul(tmp1, tmp3);
X    qfree(tmp1);
X    qfree(tmp3);
X    r->imag = qsqrt(tmp2, epsilon2);
X    qfree(tmp2);
X    qfree(epsilon2);
X    if (negimag) {
X        tmp1 = qneg(r->imag);
X        qfree(r->imag);
X        r->imag = tmp1;
X    }
X    return r;
X}
X
X
X/*
X * Calculate the complex sine within the specified accuracy.
X * This uses the formula:
X *    sin(a + bi) = sin(a) * cosh(b) + cos(a) * sinh(b) * i.
X */
XCOMPLEX *
Xcsin(c, epsilon)
X    COMPLEX *c;
X    NUMBER *epsilon;
X{
X    COMPLEX *r;
X
X    NUMBER *cosval, *coshval, *tmp1, *tmp2, *epsilon2;
X
X    if (ciszero(c))
X        return clink(&_czero_);
X    r = comalloc();
X    if (cisreal(c)) {
X        r->real = qsin(c->real, epsilon);
X        return r;
X    }
X    if (qiszero(c->real)) {
X        r->imag = qsinh(c->imag, epsilon);
X        return r;
X    }
X    epsilon2 = qscale(epsilon, -2L);
X    coshval = qcosh(c->imag, epsilon2);
X    cosval = qcos(c->real, epsilon2);
X    tmp1 = qlegtoleg(cosval, epsilon2, _sinisneg_);
X    r->real = qmul(tmp1, coshval);
X    qfree(tmp1);
X    tmp1 = qsquare(coshval);
X    qfree(coshval);
X    tmp2 = qdec(tmp1);
X    qfree(tmp1);
X    tmp1 = qsqrt(tmp2, epsilon2);
X    qfree(tmp2);
X    r->imag = qmul(tmp1, cosval);
X    qfree(tmp1);
X    qfree(cosval);
X    if (qisneg(c->imag)) {
X        tmp1 = qneg(r->imag);
X        qfree(r->imag);
X        r->imag = tmp1;
X    }
X    return r;
X}
X
X
X/*
X * Convert a number from polar coordinates to normal complex number form
X * within the specified accuracy.  This produces the value:
X *    q1 * cos(q2) + q1 * sin(q2) * i.
X */
XCOMPLEX *
Xcpolar(q1, q2, epsilon)
X    NUMBER *q1, *q2, *epsilon;
X{
X    COMPLEX *r;
X    NUMBER *tmp, *epsilon2;
X    long scale;
X
X    r = comalloc();
X    if (qiszero(q1) || qiszero(q2)) {
X        r->real = qlink(q1);
X        return r;
X    }
X    epsilon2 = epsilon;
X    if (!qisunit(q1)) {
X        scale = zhighbit(q1->num) - zhighbit(q1->den) + 1;
X        if (scale > 0)
X            epsilon2 = qscale(epsilon, -scale);
X    }
X    r->real = qcos(q2, epsilon2);
X    r->imag = qlegtoleg(r->real, epsilon2, _sinisneg_);
X    if (epsilon2 != epsilon)
X        qfree(epsilon2);
X    if (qisone(q1))
X        return r;
X    tmp = qmul(r->real, q1);
X    qfree(r->real);
X    r->real = tmp;
X    tmp = qmul(r->imag, q1);
X    qfree(r->imag);
X    r->imag = tmp;
X    return r;
X}
X
X
X/*
X * Raise one complex number to the power of another one to within the
X * specified error.
X */
XCOMPLEX *
Xcpower(c1, c2, epsilon)
X    COMPLEX *c1, *c2;
X    NUMBER *epsilon;
X{
X    COMPLEX *tmp1, *tmp2;
X    NUMBER *epsilon2;
X
X    if (cisreal(c2) && qisint(c2->real))
X        return cpowi(c1, c2->real);
X    if (cisone(c1) || ciszero(c1))
X        return clink(c1);
X    epsilon2 = qscale(epsilon, -4L);
X    tmp1 = cln(c1, epsilon2);
X    tmp2 = cmul(tmp1, c2);
X    comfree(tmp1);
X    tmp1 = cexp(tmp2, epsilon);
X    comfree(tmp2);
X    qfree(epsilon2);
X    return tmp1;
X}
X
X
X/*
X * Print a complex number.
X */
Xvoid
Xcomprint(c)
X    COMPLEX *c;
X{
X    NUMBER qtmp;
X
X    if (_outmode_ == MODE_FRAC) {
X        cprintfr(c);
X        return;
X    }
X    if (!qiszero(c->real) || qiszero(c->imag))
X        qprintnum(c->real, MODE_DEFAULT);
X    qtmp = c->imag[0];
X    if (qiszero(&qtmp))
X        return;
X    if (!qiszero(c->real) && !qisneg(&qtmp))
X        math_chr('+');
X    if (qisneg(&qtmp)) {
X        math_chr('-');
X        qtmp.num.sign = 0;
X    }
X    qprintnum(&qtmp, MODE_DEFAULT);
X    math_chr('i');
X}
X
X/* END CODE */
END_OF_FILE
if test 10332 -ne `wc -c <'comfunc.c'`; then
    echo shar: \"'comfunc.c'\" unpacked with wrong size!
fi
# end of 'comfunc.c'
fi
if test -f 'commath.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'commath.c'\"
else
echo shar: Extracting \"'commath.c'\" \(9654 characters\)
sed "s/^X//" >'commath.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 complex arithmetic primitive routines
X */
X
X#include <stdio.h>
X#include "math.h"
X
X
XCOMPLEX _czero_ =        { &_qzero_, &_qzero_, 1 };
XCOMPLEX _cone_ =        { &_qone_, &_qzero_, 1 };
Xstatic COMPLEX _cnegone_ =    { &_qnegone_, &_qzero_, 1 };
X
X#if 0
XCOMPLEX _conei_ =    { &_qzero_, &_qone_, 1 };
X#endif
X
X
X/*
X * Free list for complex numbers.
X */
Xstatic FREELIST freelist = {
X    sizeof(COMPLEX),    /* size of an item */
X    100            /* number of free items to keep */
X};
X
X
X/*
X * Add two complex numbers.
X */
XCOMPLEX *
Xcadd(c1, c2)
X    COMPLEX *c1, *c2;
X{
X    COMPLEX *r;
X
X    if (ciszero(c1))
X        return clink(c2);
X    if (ciszero(c2))
X        return clink(c1);
X    r = comalloc();
X    if (!qiszero(c1->real) || !qiszero(c2->real))
X        r->real = qadd(c1->real, c2->real);
X    if (!qiszero(c1->imag) || !qiszero(c2->imag))
X        r->imag = qadd(c1->imag, c2->imag);
X    return r;
X}
X
X
X/*
X * Subtract two complex numbers.
X */
XCOMPLEX *
Xcsub(c1, c2)
X    COMPLEX *c1, *c2;
X{
X    COMPLEX *r;
X
X    if ((c1->real == c2->real) && (c1->imag == c2->imag))
X        return clink(&_czero_);
X    if (ciszero(c2))
X        return clink(c1);
X    r = comalloc();
X    if (!qiszero(c1->real) || !qiszero(c2->real))
X        r->real = qsub(c1->real, c2->real);
X    if (!qiszero(c1->imag) || !qiszero(c2->imag))
X        r->imag = qsub(c1->imag, c2->imag);
X    return r;
X}
X
X
X/*
X * Multiply two complex numbers.
X * This saves one multiplication over the obvious algorithm by
X * trading it for several extra additions, as follows.  Let
X *    q1 = (a + b) * (c + d)
X *    q2 = a * c
X *    q3 = b * d
X * Then (a+bi) * (c+di) = (q2 - q3) + (q1 - q2 - q3)i.
X */
XCOMPLEX *
Xcmul(c1, c2)
X    COMPLEX *c1, *c2;
X{
X    COMPLEX *r;
X    NUMBER *q1, *q2, *q3, *q4;
X
X    if (ciszero(c1) || ciszero(c2))
X        return clink(&_czero_);
X    if (cisone(c1))
X        return clink(c2);
X    if (cisone(c2))
X        return clink(c1);
X    if (cisreal(c2))
X        return cmulq(c1, c2->real);
X    if (cisreal(c1))
X        return cmulq(c2, c1->real);
X    /*
X     * Need to do the full calculation.
X     */
X    r = comalloc();
X    q2 = qadd(c1->real, c1->imag);
X    q3 = qadd(c2->real, c2->imag);
X    q1 = qmul(q2, q3);
X    qfree(q2);
X    qfree(q3);
X    q2 = qmul(c1->real, c2->real);
X    q3 = qmul(c1->imag, c2->imag);
X    q4 = qadd(q2, q3);
X    r->real = qsub(q2, q3);
X    r->imag = qsub(q1, q4);
X    qfree(q1);
X    qfree(q2);
X    qfree(q3);
X    qfree(q4);
X    return r;
X}
X
X
X/*
X * Square a complex number.
X */
XCOMPLEX *
Xcsquare(c)
X    COMPLEX *c;
X{
X    COMPLEX *r;
X    NUMBER *q1, *q2;
X
X    if (ciszero(c))
X        return clink(&_czero_);
X    if (cisrunit(c))
X        return clink(&_cone_);
X    if (cisiunit(c))
X        return clink(&_cnegone_);
X    r = comalloc();
X    if (cisreal(c)) {
X        r->real = qsquare(c->real);
X        return r;
X    }
X    if (cisimag(c)) {
X        q1 = qsquare(c->imag);
X        r->real = qneg(q1);
X        qfree(q1);
X        return r;
X    }
X    q1 = qsquare(c->real);
X    q2 = qsquare(c->imag);
X    r->real = qsub(q1, q2);
X    qfree(q1);
X    qfree(q2);
X    q1 = qmul(c->real, c->imag);
X    r->imag = qscale(q1, 1L);
X    qfree(q1);
X    return r;
X}
X
X
X/*
X * Divide two complex numbers.
X */
XCOMPLEX *
Xcdiv(c1, c2)
X    COMPLEX *c1, *c2;
X{
X    COMPLEX *r;
X    NUMBER *q1, *q2, *q3, *den;
X
X    if (ciszero(c2))
X        error("Division by zero");
X    if ((c1->real == c2->real) && (c1->imag == c2->imag))
X        return clink(&_cone_);
X    r = comalloc();
X    if (cisreal(c1) && cisreal(c2)) {
X        r->real = qdiv(c1->real, c2->real);
X        return r;
X    }
X    if (cisimag(c1) && cisimag(c2)) {
X        r->real = qdiv(c1->imag, c2->imag);
X        return r;
X    }
X    if (cisimag(c1) && cisreal(c2)) {
X        r->imag = qdiv(c1->imag, c2->real);
X        return r;
X    }
X    if (cisreal(c1) && cisimag(c2)) {
X        q1 = qdiv(c1->real, c2->imag);
X        r->imag = qneg(q1);
X        qfree(q1);
X        return r;
X    }
X    if (cisreal(c2)) {
X        r->real = qdiv(c1->real, c2->real);
X        r->imag = qdiv(c1->imag, c2->real);
X        return r;
X    }
X    q1 = qsquare(c2->real);
X    q2 = qsquare(c2->imag);
X    den = qadd(q1, q2);
X    qfree(q1);
X    qfree(q2);
X    q1 = qmul(c1->real, c2->real);
X    q2 = qmul(c1->imag, c2->imag);
X    q3 = qadd(q1, q2);
X    qfree(q1);
X    qfree(q2);
X    r->real = qdiv(q3, den);
X    qfree(q3);
X    q1 = qmul(c1->real, c2->imag);
X    q2 = qmul(c1->imag, c2->real);
X    q3 = qsub(q2, q1);
X    qfree(q1);
X    qfree(q2);
X    r->imag = qdiv(q3, den);
X    qfree(q3);
X    qfree(den);
X    return r;
X}
X
X
X/*
X * Invert a complex number.
X */
XCOMPLEX *
Xcinv(c)
X    COMPLEX *c;
X{
X    COMPLEX *r;
X    NUMBER *q1, *q2, *den;
X
X    if (ciszero(c))
X        error("Inverting zero");
X    r = comalloc();
X    if (cisreal(c)) {
X        r->real = qinv(c->real);
X        return r;
X    }
X    if (cisimag(c)) {
X        q1 = qinv(c->imag);
X        r->imag = qneg(q1);
X        qfree(q1);
X        return r;
X    }
X    q1 = qsquare(c->real);
X    q2 = qsquare(c->imag);
X    den = qadd(q1, q2);
X    qfree(q1);
X    qfree(q2);
X    r->real = qdiv(c->real, den);
X    q1 = qdiv(c->imag, den);
X    r->imag = qneg(q1);
X    qfree(q1);
X    qfree(den);
X    return r;
X}
X
X
X/*
X * Negate a complex number.
X */
XCOMPLEX *
Xcneg(c)
X    COMPLEX *c;
X{
X    COMPLEX *r;
X
X    if (ciszero(c))
X        return clink(&_czero_);
X    r = comalloc();
X    if (!qiszero(c->real))
X        r->real = qneg(c->real);
X    if (!qiszero(c->imag))
X        r->imag = qneg(c->imag);
X    return r;
X}
X
X
X/*
X * Take the integer part of a complex number.
X * This means take the integer part of both components.
X */
XCOMPLEX *
Xcint(c)
X    COMPLEX *c;
X{
X    COMPLEX *r;
X
X    if (cisint(c))
X        return clink(c);
X    r = comalloc();
X    r->real = qint(c->real);
X    r->imag = qint(c->imag);
X    return r;
X}
X
X
X/*
X * Take the fractional part of a complex number.
X * This means take the fractional part of both components.
X */
XCOMPLEX *
Xcfrac(c)
X    COMPLEX *c;
X{
X    COMPLEX *r;
X
X    if (cisint(c))
X        return clink(&_czero_);
X    r = comalloc();
X    r->real = qfrac(c->real);
X    r->imag = qfrac(c->imag);
X    return r;
X}
X
X
X#if 0
X/*
X * Take the conjugate of a complex number.
X * This negates the complex part.
X */
XCOMPLEX *
Xcconj(c)
X    COMPLEX *c;
X{
X    COMPLEX *r;
X
X    if (cisreal(c))
X        return clink(c);
X    r = comalloc();
X    if (!qiszero(c->real))
X        r->real = qlink(c->real);
X    r->imag = qneg(c->imag);
X    return r;
X}
X
X
X/*
X * Return the real part of a complex number.
X */
XCOMPLEX *
Xcreal(c)
X    COMPLEX *c;
X{
X    COMPLEX *r;
X
X    if (cisreal(c))
X        return clink(c);
X    r = comalloc();
X    if (!qiszero(c->real))
X        r->real = qlink(c->real);
X    return r;
X}
X
X
X/*
X * Return the imaginary part of a complex number as a real.
X */
XCOMPLEX *
Xcimag(c)
X    COMPLEX *c;
X{
X    COMPLEX *r;
X
X    if (cisreal(c))
X        return clink(&_czero_);
X    r = comalloc();
X    r->real = qlink(c->imag);
X    return r;
X}
X#endif
X
X
X/*
X * Add a real number to a complex number.
X */
XCOMPLEX *
Xcaddq(c, q)
X    COMPLEX *c;
X    NUMBER *q;
X{
X    COMPLEX *r;
X
X    if (qiszero(q))
X        return clink(c);
X    r = comalloc();
X    r->real = qadd(c->real, q);
X    r->imag = qlink(c->imag);
X    return r;
X}
X
X
X/*
X * Subtract a real number from a complex number.
X */
XCOMPLEX *
Xcsubq(c, q)
X    COMPLEX *c;
X    NUMBER *q;
X{
X    COMPLEX *r;
X
X    if (qiszero(q))
X        return clink(c);
X    r = comalloc();
X    r->real = qsub(c->real, q);
X    r->imag = qlink(c->imag);
X    return r;
X}
X
X
X/*
X * Shift the components of a complex number left by the specified
X * number of bits.  Negative values shift to the right.
X */
XCOMPLEX *
Xcshift(c, n)
X    COMPLEX *c;
X    long n;
X{
X    COMPLEX *r;
X
X    if (ciszero(c) || (n == 0))
X        return clink(c);
X    r = comalloc();
X    r->real = qshift(c->real, n);
X    r->imag = qshift(c->imag, n);
X    return r;
X}
X
X
X/*
X * Scale a complex number by a power of two.
X */
XCOMPLEX *
Xcscale(c, n)
X    COMPLEX *c;
X    long n;
X{
X    COMPLEX *r;
X
X    if (ciszero(c) || (n == 0))
X        return clink(c);
X    r = comalloc();
X    r->real = qscale(c->real, n);
X    r->imag = qscale(c->imag, n);
X    return r;
X}
X
X
X/*
X * Multiply a complex number by a real number.
X */
XCOMPLEX *
Xcmulq(c, q)
X    COMPLEX *c;
X    NUMBER *q;
X{
X    COMPLEX *r;
X
X    if (qiszero(q))
X        return clink(&_czero_);
X    if (qisone(q))
X        return clink(c);
X    if (qisnegone(q))
X        return cneg(c);
X    r = comalloc();
X    r->real = qmul(c->real, q);
X    r->imag = qmul(c->imag, q);
X    return r;
X}
X
X
X/*
X * Divide a complex number by a real number.
X */
XCOMPLEX *
Xcdivq(c, q)
X    COMPLEX *c;
X    NUMBER *q;
X{
X    COMPLEX *r;
X
X    if (qiszero(q))
X        error("Division by zero");
X    if (qisone(q))
X        return clink(c);
X    if (qisnegone(q))
X        return cneg(c);
X    r = comalloc();
X    r->real = qdiv(c->real, q);
X    r->imag = qdiv(c->imag, q);
X    return r;
X}
X
X
X/*
X * Take the integer quotient of a complex number by a real number.
X * This is defined to be the result of doing the quotient for each component.
X */
XCOMPLEX *
Xcquoq(c, q)
X    COMPLEX *c;
X    NUMBER *q;
X{
X    COMPLEX *r;
X
X    if (qiszero(q))
X        error("Division by zero");
X    r = comalloc();
X    r->real = qquo(c->real, q);
X    r->imag = qquo(c->imag, q);
X    return r;
X}
X
X
X/*
X * Take the modulus of a complex number by a real number.
X * This is defined to be the result of doing the modulo for each component.
X */
XCOMPLEX *
Xcmodq(c, q)
X    COMPLEX *c;
X    NUMBER *q;
X{
X    COMPLEX *r;
X
X    if (qiszero(q))
X        error("Division by zero");
X    r = comalloc();
X    r->real = qmod(c->real, q);
X    r->imag = qmod(c->imag, q);
X    return r;
X}
X
X
X#if 0
X/*
X * Construct a complex number given the real and imaginary components.
X */
XCOMPLEX *
Xqqtoc(q1, q2)
X    NUMBER *q1, *q2;
X{
X    COMPLEX *r;
X
X    if (qiszero(q1) && qiszero(q2))
X        return clink(&_czero_);
X    r = comalloc();
X    if (!qiszero(q1))
X        r->real = qlink(q1);
X    if (!qiszero(q2))
X        r->imag = qlink(q2);
X    return r;
X}
X#endif
X
X
X/*
X * Compare two complex numbers for equality, returning FALSE if they are equal,
X * and TRUE if they differ.
X */
XBOOL
Xccmp(c1, c2)
X    COMPLEX *c1, *c2;
X{
X    BOOL i;
X
X    i = qcmp(c1->real, c2->real);
X    if (!i)
X        i = qcmp(c1->imag, c2->imag);
X    return i;
X}
X
X
X/*
X * Allocate a new complex number.
X */
XCOMPLEX *
Xcomalloc()
X{
X    COMPLEX *r;
X
X    r = (COMPLEX *) allocitem(&freelist);
X    if (r == NULL)
X        error("Cannot allocate complex number");
X    r->links = 1;
X    r->real = qlink(&_qzero_);
X    r->imag = qlink(&_qzero_);
X    return r;
X}
X
X
X/*
X * Free a complex number.
X */
Xvoid
Xcomfree(c)
X    COMPLEX *c;
X{
X    if (--(c->links) > 0)
X        return;
X    qfree(c->real);
X    qfree(c->imag);
X    freeitem(&freelist, (FREEITEM *) c);
X}
X
X/* END CODE */
END_OF_FILE
if test 9654 -ne `wc -c <'commath.c'`; then
    echo shar: \"'commath.c'\" unpacked with wrong size!
fi
# end of 'commath.c'
fi
if test -f 'file.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'file.c'\"
else
echo shar: Extracting \"'file.c'\" \(10407 characters\)
sed "s/^X//" >'file.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 * File I/O routines callable by users.
X */
X
X#include "stdarg.h"
X#include "calc.h"
X
X
X#define    READSIZE    1024    /* buffer size for reading */
X
X/*
X * Definition of opened files.
X */
Xtypedef struct {
X    FILEID id;        /* id to identify this file */
X    FILE *fp;        /* real file structure for I/O */
X    char *name;        /* file name */
X    BOOL reading;        /* TRUE if opened for reading */
X    BOOL writing;        /* TRUE if opened for writing */
X    char *mode;        /* open mode */
X} FILEIO;
X
X
X/*
X * Table of opened files.
X * The first three entries always correspond to stdin, stdout, and stderr,
X * and cannot be closed.  Their file ids are always 0, 1, and 2.
X */
Xstatic FILEIO files[MAXFILES] = {
X    FILEID_STDIN,  stdin,  "(stdin)",  TRUE, FALSE, "reading",
X    FILEID_STDOUT, stdout, "(stdout)", FALSE, TRUE, "writing",
X    FILEID_STDERR, stderr, "(stderr)", FALSE, TRUE, "writing"
X};
X
Xstatic FILEID lastid = FILEID_STDERR;        /* last allocated file id */
X
X
X
X/*
X * Open the specified file name for reading or writing as determined by
X * the specified mode ("r", "w", or "a").  Returns a file id which can be
X * used to do I/O to the file, or else FILEID_NONE if the open failed.
X * Aborts with an error if too many files are opened or the mode is illegal.
X */
XFILEID
Xopenid(name, mode)
X    char *name;        /* file name */
X    char *mode;        /* open mode */
X{
X    FILEIO *fiop;        /* file structure */
X    FILEID id;        /* new file id */
X    int count;
X
X    if (((*mode != 'r') && (*mode != 'w') && (*mode != 'a')) || mode[1])
X        error("Illegal mode for fopen");
X
X    count = MAXFILES;
X    do {
X        if (--count < 0)
X            error("Too many open files");
X        id = ++lastid;
X        fiop = &files[id % MAXFILES];
X
X    } while (fiop->reading || fiop->writing);
X
X    fiop->name = (char *)malloc(strlen(name) + 1);
X    if (fiop->name == NULL) {
X        lastid--;
X        error("No memory for filename");
X    }
X    strcpy(fiop->name, name);
X
X    fiop->fp = f_open(name, mode);
X    if (fiop->fp == NULL) {
X        free(fiop->name);
X        fiop->name = NULL;
X        lastid--;
X        return FILEID_NONE;
X    }
X
X    switch (*mode) {
X        case 'r':
X            fiop->mode = "reading";
X            fiop->reading = TRUE;
X            break;
X        case 'w':
X            fiop->mode = "writing";
X            fiop->writing = TRUE;
X            break;
X        case 'a':
X            fiop->mode = "appending";
X            fiop->writing = TRUE;
X            break;
X    }
X
X    fiop->id = id;
X
X    return id;
X}
X
X
X/*
X * Find the file I/O structure for the specified file id, and verify that
X * it is opened in the required manner ('r' for reading or 'w' for writing).
X * If mode is 0, then no open checks are made at all, and NULL is then
X * returned if the id represents a closed file.
X */
Xstatic FILEIO *
Xfindid(id, mode)
X    FILEID id;
X{
X    FILEIO *fiop;        /* file structure */
X    char *msg;
X    BOOL flag;
X
X    if ((id < 0) || (id > lastid))
X        error("Illegal file id");
X
X    fiop = &files[id % MAXFILES];
X
X    switch (mode) {
X        case 'r':
X            msg = "Reading from";
X            flag = fiop->reading;
X            break;
X        case 'w':
X            msg = "Writing to";
X            flag = fiop->writing;
X            break;
X        case 0:
X            msg = NULL;
X            break;
X        default:
X            error("Unknown findid mode");
X    }
X
X    if (fiop->id != id) {
X        if (msg)
X            error("%s closed file", msg);
X        return NULL;
X    }
X
X    if (msg && !flag)
X        error("%s file not opened that way", msg);
X
X    return fiop;
X}
X
X
X/*
X * Return whether or not a file id is valid.  This is used for if tests.
X */
XBOOL
Xvalidid(id)
X    FILEID id;
X{
X    return (findid(id, 0) != NULL);
X}
X
X
X/*
X * Return the file id for the entry in the file table at the specified index.
X * Returns FILEID_NONE if the index is illegal or the file is closed.
X */
XFILEID
Xindexid(index)
X    long index;
X{
X    FILEIO *fiop;        /* file structure */
X
X    if ((index < 0) || (index >= MAXFILES))
X        return FILEID_NONE;
X
X    fiop = &files[index];
X    if (fiop->reading || fiop->writing)
X        return fiop->id;
X
X    return FILEID_NONE;
X}
X
X
X/*
X * Close the specified file id.  Returns TRUE if there was an error.
X * Closing of stdin, stdout, or stderr is illegal, but closing of already
X * closed files is allowed.
X */
XBOOL
Xcloseid(id)
X    FILEID id;
X{
X    FILEIO *fiop;        /* file structure */
X    int err;
X
X    if ((id == FILEID_STDIN) || (id == FILEID_STDOUT) ||
X        (id == FILEID_STDERR))
X            error("Cannot close stdin, stdout, or stderr");
X
X    fiop = findid(id, 0);
X    if (fiop == NULL)
X        return FALSE;
X
X    fiop->id = FILEID_NONE;
X    if (!fiop->reading && !fiop->writing)
X        error("Closing non-opened file");
X    fiop->reading = FALSE;
X    fiop->writing = FALSE;
X
X    if (fiop->name)
X        free(fiop->name);
X    fiop->name = NULL;
X
X    err = ferror(fiop->fp);
X    err |= fclose(fiop->fp);
X    fiop->fp = NULL;
X
X    return (err != 0);
X}
X
X
X/*
X * Return whether or not an error occurred to a file.
X */
XBOOL
Xerrorid(id)
X    FILEID id;
X{
X    FILEIO *fiop;        /* file structure */
X
X    fiop = findid(id, 0);
X    if (fiop == NULL)
X        error("Closed file for ferror");
X    return (ferror(fiop->fp) != 0);
X}
X
X
X/*
X * Return whether or not end of file occurred to a file.
X */
XBOOL
Xeofid(id)
X    FILEID id;
X{
X    FILEIO *fiop;        /* file structure */
X
X    fiop = findid(id, 0);
X    if (fiop == NULL)
X        error("Closed file for feof");
X    return (feof(fiop->fp) != 0);
X}
X
X
X/*
X * Flush output to an opened file.
X */
Xvoid
Xflushid(id)
X    FILEID id;
X{
X    FILEIO *fiop;        /* file structure */
X
X    fiop = findid(id, 'w');
X    fflush(fiop->fp);
X}
X
X
X/*
X * Read the next line from an opened file.
X * Returns a pointer to an allocated string holding the null-terminated
X * line (without any terminating newline), or else a NULL pointer on an
X * end of file or error.
X */
Xvoid
Xreadid(id, retptr)
X    FILEID    id;        /* file to read from */
X    char **retptr;        /* returned pointer to string */
X{
X    FILEIO *fiop;        /* file structure */
X    char *str;        /* current string */
X    int len;        /* current length of string */
X    int totlen;        /* total length of string */
X    char buf[READSIZE];    /* temporary buffer */
X
X    totlen = 0;
X    str = NULL;
X
X    fiop = findid(id, 'r');
X
X    while (fgets(buf, READSIZE, fiop->fp) && buf[0]) {
X        len = strlen(buf);
X        if (totlen)
X            str = (char *)realloc(str, totlen + len + 1);
X        else
X            str = (char *)malloc(len + 1);
X        if (str == NULL)
X            error("No memory in freadline");
X        strcpy(&str[totlen], buf);
X        totlen += len;
X        if (buf[len - 1] == '\n') {
X            str[totlen - 1] = '\0';
X            *retptr = str;
X            return;
X        }
X    }
X    if (totlen && ferror(fiop->fp)) {
X        free(str);
X        str = NULL;
X    }
X    *retptr = str;
X}
X
X
X/*
X * Return the next character from an opened file.
X * Returns EOF if there was an error or end of file.
X */
Xint
Xgetcharid(id)
X    FILEID id;
X{
X    return fgetc(findid(id, 'r')->fp);
X}
X
X
X/*
X * Print out the name of an opened file.
X * If the file has been closed, a null name is printed.
X * If flags contain PRINT_UNAMBIG then extra information is printed
X * identifying the output as a file and some data about it.
X */
Xvoid
Xprintid(id, flags)
X    FILEID id;
X{
X    FILEIO *fiop;        /* file structure */
X    FILE *fp;
X
X    fiop = findid(id, 0);
X    if (fiop == NULL) {
X        math_str((flags & PRINT_UNAMBIG) ? "FILE (closed)" : "\"\"");
X        return;
X    }
X    if ((flags & PRINT_UNAMBIG) == 0) {
X        math_chr('"');
X        math_str(fiop->name);
X        math_chr('"');
X        return;
X    }
X
X    fp = fiop->fp;
X    math_fmt("FILE \"%s\" (%s, pos %ld", fiop->name,  fiop->mode,
X        ftell(fp));
X    if (ferror(fp))
X        math_str(", error");
X    if (feof(fp))
X        math_str(", eof");
X    math_chr(')');
X}
X
X
X/*
X * Print a formatted string similar to printf.  Various formats of output
X * are possible, depending on the format string AND the actual types of the
X * values.  Mismatches do not cause errors, instead something reasonable is
X * printed instead.  The output goes to the file with the specified id.
X */
Xvoid
Xidprintf(id, fmt, count, vals)
X    FILEID id;            /* file id to print to */
X    char *fmt;            /* standard format string */
X    VALUE **vals;            /* table of values to print */
X{
X    FILEIO *fiop;
X    VALUE *vp;
X    char *str;
X    int ch, len;
X    int oldmode, newmode;
X    long olddigits, newdigits;
X    long width, precision;
X    BOOL didneg, didprecision;
X
X    fiop = findid(id, 'w');
X
X    setfp(fiop->fp);
X
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                    setfp(stdout);
X                    return;
X            }
X            math_chr(ch);
X            continue;
X        }
X
X        if (ch != '%') {
X            math_chr(ch);
X            continue;
X        }
X
X        /*
X         * Here to handle formats.
X         */
X        didneg = FALSE;
X        didprecision = FALSE;
X        width = 0;
X        precision = 0;
X
X        ch = *fmt++;
X        if (ch == '-') {
X            didneg = TRUE;
X            ch = *fmt++;
X        }
X        while ((ch >= '0') && (ch <= '9')) {
X            width = width * 10 + (ch - '0');
X            ch = *fmt++;
X        }
X        if (ch == '.') {
X            didprecision = TRUE;
X            ch = *fmt++;
X            while ((ch >= '0') && (ch <= '9')) {
X                precision = precision * 10 + (ch - '0');
X                ch = *fmt++;
X            }
X        }
X        if (ch == 'l')
X            ch = *fmt++;
X
X        oldmode = _outmode_;
X        newmode = oldmode;
X        olddigits = _outdigits_;
X        newdigits = olddigits;
X        if (didprecision)
X            newdigits = precision;
X
X        switch (ch) {
X            case 'd':
X            case 's':
X            case 'c':
X                break;
X            case 'f':
X                newmode = MODE_REAL;
X                break;
X            case 'e':
X                newmode = MODE_EXP;
X                break;
X            case 'r':
X                newmode = MODE_FRAC;
X                break;
X            case 'o':
X                newmode = MODE_OCTAL;
X                break;
X            case 'x':
X                newmode = MODE_HEX;
X                break;
X            case 'b':
X                newmode = MODE_BINARY;
X                break;
X            case 0:
X                setfp(stdout);
X                return;
X            default:
X                math_chr(ch);
X                continue;
X        }
X
X        if (--count < 0)
X            error("Not enough arguments for fprintf");
X        vp = *vals++;
X
X        setdigits(newdigits);
X        setmode(newmode);
X
X        /*
X         * If there is no width specification, or if the type of
X         * value requires multiple lines, then just output the
X         * value directly.
X         */
X        if ((width == 0) ||
X            (vp->v_type == V_MAT) || (vp->v_type == V_LIST))
X        {
X            printvalue(vp, PRINT_NORMAL);
X            setmode(oldmode);
X            setdigits(olddigits);
X            continue;
X        }
X
X        /*
X         * There is a field width.  Collect the output in a string,
X         * print it padded appropriately with spaces, and free it.
X         * However, if the output contains a newline, then ignore
X         * the field width.
X         */
X        divertio();
X        printvalue(vp, PRINT_NORMAL);
X        str = getdivertedio();
X        if (strchr(str, '\n'))
X            width = 0;
X        len = strlen(str);
X        while (!didneg && (width > len)) {
X            width--;
X            math_chr(' ');
X        }
X        math_str(str);
X        free(str);
X        while (didneg && (width > len)) {
X            width--;
X            math_chr(' ');
X        }
X        setmode(oldmode);
X        setdigits(olddigits);
X    }
X    setfp(stdout);
X}
X
X/* END CODE */
END_OF_FILE
if test 10407 -ne `wc -c <'file.c'`; then
    echo shar: \"'file.c'\" unpacked with wrong size!
fi
# end of 'file.c'
fi
if test -f 'listfunc.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'listfunc.c'\"
else
echo shar: Extracting \"'listfunc.c'\" \(11116 characters\)
sed "s/^X//" >'listfunc.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 * List handling routines.
X * Lists can be composed of any types of values, mixed if desired.
X * Lists are doubly linked so that elements can be inserted or
X * deleted efficiently at any point in the list.  A pointer is
X * kept to the most recently indexed element so that sequential
X * accesses are fast.
X */
X
X#include "calc.h"
X
X
Xstatic LISTELEM *elemalloc();
Xstatic LISTELEM *listelement();
Xstatic void elemfree();
Xstatic void removelistelement();
X
X
X/*
X * Free lists for list headers and list elements.
X */
Xstatic FREELIST    headerfreelist = {
X    sizeof(LIST),        /* size of list header */
X    20            /* number of free headers to keep */
X};
X
Xstatic FREELIST elementfreelist = {
X    sizeof(LISTELEM),    /* size of list element */
X    1000            /* number of free list elements to keep */
X};
X
X
X/*
X * Insert an element before the first element of a list.
X */
Xvoid
Xinsertlistfirst(lp, vp)
X    LIST *lp;        /* list to put element onto */
X    VALUE *vp;        /* value to be inserted */
X{
X    LISTELEM *ep;        /* list element */
X
X    ep = elemalloc();
X    copyvalue(vp, &ep->e_value);
X    if (lp->l_count == 0)
X        lp->l_last = ep;
X    else {
X        lp->l_cacheindex++;
X        lp->l_first->e_prev = ep;
X        ep->e_next = lp->l_first;
X    }
X    lp->l_first = ep;
X    lp->l_count++;
X}
X
X
X/*
X * Insert an element after the last element of a list.
X */
Xvoid
Xinsertlistlast(lp, vp)
X    LIST *lp;        /* list to put element onto */
X    VALUE *vp;        /* value to be inserted */
X{
X    LISTELEM *ep;        /* list element */
X
X    ep = elemalloc();
X    copyvalue(vp, &ep->e_value);
X    if (lp->l_count == 0)
X        lp->l_first = ep;
X    else {
X        lp->l_last->e_next = ep;
X        ep->e_prev = lp->l_last;
X    }
X    lp->l_last = ep;
X    lp->l_count++;
X}
X
X
X/*
X * Insert an element into the middle of list at the given index (zero based).
X * The specified index will select the new element, so existing elements
X * at or beyond the index will be shifted down one position.  It is legal
X * to specify an index which is right at the end of the list, in which
X * case the element is appended to the list.
X */
Xvoid
Xinsertlistmiddle(lp, index, vp)
X    LIST *lp;        /* list to put element onto */
X    long index;        /* element number to insert in front of */
X    VALUE *vp;        /* value to be inserted */
X{
X    LISTELEM *ep;        /* list element */
X    LISTELEM *oldep;    /* old list element at desired index */
X
X    if (index == 0) {
X        insertlistfirst(lp, vp);
X        return;
X    }
X    if (index == lp->l_count) {
X        insertlistlast(lp, vp);
X        return;
X    }
X    oldep = NULL;
X    if ((index >= 0) && (index < lp->l_count))
X        oldep = listelement(lp, index);
X    if (oldep == NULL)
X        error("Index out of bounds for list insertion");
X    ep = elemalloc();
X    copyvalue(vp, &ep->e_value);
X    ep->e_next = oldep;
X    ep->e_prev = oldep->e_prev;
X    ep->e_prev->e_next = ep;
X    oldep->e_prev = ep;
X    lp->l_cache = ep;
X    lp->l_cacheindex = index;
X    lp->l_count++;
X}
X
X
X/*
X * Remove the first element from a list, returning its value.
X * Returns the null value if no more elements exist.
X */
Xvoid
Xremovelistfirst(lp, vp)
X    LIST *lp;        /* list to have element removed */
X    VALUE *vp;        /* location of the value */
X{
X    if (lp->l_count == 0) {
X        vp->v_type = V_NULL;
X        return;
X    }
X    *vp = lp->l_first->e_value;
X    lp->l_first->e_value.v_type = V_NULL;
X    removelistelement(lp, lp->l_first);
X}
X
X
X/*
X * Remove the last element from a list, returning its value.
X * Returns the null value if no more elements exist.
X */
Xvoid
Xremovelistlast(lp, vp)
X    LIST *lp;        /* list to have element removed */
X    VALUE *vp;        /* location of the value */
X{
X    if (lp->l_count == 0) {
X        vp->v_type = V_NULL;
X        return;
X    }
X    *vp = lp->l_last->e_value;
X    lp->l_last->e_value.v_type = V_NULL;
X    removelistelement(lp, lp->l_last);
X}
X
X
X/*
X * Remove the element with the given index from a list, returning its value.
X */
Xvoid
Xremovelistmiddle(lp, index, vp)
X    LIST *lp;        /* list to have element removed */
X    long index;        /* list element to be removed */
X    VALUE *vp;        /* location of the value */
X{
X    LISTELEM *ep;        /* element being removed */
X
X    ep = NULL;
X    if ((index >= 0) && (index < lp->l_count))
X        ep = listelement(lp, index);
X    if (ep == NULL)
X        error("Index out of bounds for list deletion");
X    *vp = ep->e_value;
X    ep->e_value.v_type = V_NULL;
X    removelistelement(lp, ep);
X}
X
X
X/*
X * Remove an arbitrary element from a list.
X * The value contained in the element is freed.
X */
Xstatic void
Xremovelistelement(lp, ep)
X    register LIST *lp;        /* list header */
X    register LISTELEM *ep;        /* list element to remove */
X{
X    if ((ep == lp->l_cache) || ((ep != lp->l_first) && (ep != lp->l_last)))
X        lp->l_cache = NULL;
X    if (ep->e_next)
X        ep->e_next->e_prev = ep->e_prev;
X    if (ep->e_prev)
X        ep->e_prev->e_next = ep->e_next;
X    if (ep == lp->l_first) {
X        lp->l_first = ep->e_next;
X        lp->l_cacheindex--;
X    }
X    if (ep == lp->l_last)
X        lp->l_last = ep->e_prev;
X    lp->l_count--;
X    elemfree(ep);
X}
X
X
X/*
X * Search a list for the specified value starting at the specified index.
X * Returns the element number (zero based) of the found value, or -1 if
X * the value was not found.
X */
Xlong
Xlistsearch(lp, vp, index)
X    LIST *lp;
X    VALUE *vp;
X    long index;
X{
X    register LISTELEM *ep;
X
X    if (index < 0)
X        index = 0;
X    ep = listelement(lp, index);
X    while (ep) {
X        if (!comparevalue(&ep->e_value, vp)) {
X            lp->l_cache = ep;
X            lp->l_cacheindex = index;
X            return index;
X        }
X        ep = ep->e_next;
X        index++;
X    }
X    return -1;
X}
X
X
X/*
X * Search a list backwards for the specified value starting at the
X * specified index.  Returns the element number (zero based) of the
X * found value, or -1 if the value was not found.
X */
Xlong
Xlistrsearch(lp, vp, index)
X    LIST *lp;
X    VALUE *vp;
X    long index;
X{
X    register LISTELEM *ep;
X
X    if (index >= lp->l_count)
X        index = lp->l_count - 1;
X    ep = listelement(lp, index);
X    while (ep) {
X        if (!comparevalue(&ep->e_value, vp)) {
X            lp->l_cache = ep;
X            lp->l_cacheindex = index;
X            return index;
X        }
X        ep = ep->e_prev;
X        index--;
X    }
X    return -1;
X}
X
X
X/*
X * Index into a list and return the address for the value corresponding
X * to that index.  Returns NULL if the element does not exist.
X */
XVALUE *
Xlistindex(lp, index)
X    LIST *lp;        /* list to index into */
X    long index;        /* index of desired element */
X{
X    LISTELEM *ep;
X
X    ep = listelement(lp, index);
X    if (ep == NULL)
X        return NULL;
X    return &ep->e_value;
X}
X
X
X/*
X * Return the element at a specified index number of a list.
X * The list is indexed starting at zero, and negative indices
X * indicate to index from the end of the list.  This routine finds
X * the element by chaining through the list from the closest one
X * of the first, last, and cached elements.  Returns NULL if the
X * element does not exist.
X */
Xstatic LISTELEM *
Xlistelement(lp, index)
X    register LIST *lp;    /* list to index into */
X    long index;        /* index of desired element */
X{
X    register LISTELEM *ep;    /* current list element */
X    long dist;        /* distance to element */
X    long temp;        /* temporary distance */
X    BOOL forward;        /* TRUE if need to walk forwards */
X
X    if (index < 0)
X        index += lp->l_count;
X    if ((index < 0) || (index >= lp->l_count))
X        return NULL;
X    /*
X     * Check quick special cases first.
X     */
X    if (index == 0)
X        return lp->l_first;
X    if (index == 1)
X        return lp->l_first->e_next;
X    if (index == lp->l_count - 1)
X        return lp->l_last;
X    if ((index == lp->l_cacheindex) && lp->l_cache)
X        return lp->l_cache;
X    /*
X     * Calculate whether it is better to go forwards from
X     * the first element or backwards from the last element.
X     */
X    forward = ((index * 2) <= lp->l_count);
X    if (forward) {
X        dist = index;
X        ep = lp->l_first;
X    } else {
X        dist = (lp->l_count - 1) - index;
X        ep = lp->l_last;
X    }
X    /*
X     * Now see if we have a cached element and if so, whether or
X     * not the distance from it is better than the above distance.
X     */
X    if (lp->l_cache) {
X        temp = index - lp->l_cacheindex;
X        if ((temp >= 0) && (temp < dist)) {
X            dist = temp;
X            ep = lp->l_cache;
X            forward = TRUE;
X        }
X        if ((temp < 0) && (-temp < dist)) {
X            dist = -temp;
X            ep = lp->l_cache;
X            forward = FALSE;
X        }
X    }
X    /*
X     * Now walk forwards or backwards from the selected element
X     * until we reach the correct element.  Cache the location of
X     * the found element for future use.
X     */
X    if (forward) {
X        while (dist-- > 0)
X            ep = ep->e_next;
X    } else {
X        while (dist-- > 0)
X            ep = ep->e_prev;
X    }
X    lp->l_cache = ep;
X    lp->l_cacheindex = index;
X    return ep;
X}
X
X
X/*
X * Compare two lists to see if they are identical.
X * Returns TRUE if they are different.
X */
XBOOL
Xlistcmp(lp1, lp2)
X    LIST *lp1, *lp2;
X{
X    LISTELEM *e1, *e2;
X    long count;
X
X    if (lp1 == lp2)
X        return FALSE;
X    if (lp1->l_count != lp2->l_count)
X        return TRUE;
X    e1 = lp1->l_first;
X    e2 = lp2->l_first;
X    count = lp1->l_count;
X    while (count-- > 0) {
X        if (comparevalue(&e1->e_value, &e2->e_value))
X            return TRUE;
X        e1 = e1->e_next;
X        e2 = e2->e_next;
X    }
X    return FALSE;
X}
X
X
X/*
X * Copy a list
X */
XLIST *
Xlistcopy(oldlp)
X    LIST *oldlp;
X{
X    LIST *lp;
X    LISTELEM *oldep;
X
X    lp = listalloc();
X    oldep = oldlp->l_first;
X    while (oldep) {
X        insertlistlast(lp, &oldep->e_value);
X        oldep = oldep->e_next;
X    }
X    return lp;
X}
X
X
X/*
X * Allocate an element for a list.
X */
Xstatic LISTELEM *
Xelemalloc()
X{
X    LISTELEM *ep;
X
X    ep = (LISTELEM *) allocitem(&elementfreelist);
X    if (ep == NULL)
X        error("Cannot allocate list element");
X    ep->e_next = NULL;
X    ep->e_prev = NULL;
X    ep->e_value.v_type = V_NULL;
X    return ep;
X}
X
X
X/*
X * Free a list element, along with any contained value.
X */
Xstatic void
Xelemfree(ep)
X    LISTELEM *ep;
X{
X    if (ep->e_value.v_type != V_NULL)
X        freevalue(&ep->e_value);
X    freeitem(&elementfreelist, (FREEITEM *) ep);
X}
X
X
X/*
X * Allocate a new list header.
X */
XLIST *
Xlistalloc()
X{
X    register LIST *lp;
X
X    lp = (LIST *) allocitem(&headerfreelist);
X    if (lp == NULL)
X        error("Cannot allocate list header");
X    lp->l_first = NULL;
X    lp->l_last = NULL;
X    lp->l_cache = NULL;
X    lp->l_cacheindex = 0;
X    lp->l_count = 0;
X    return lp;
X}
X
X
X/*
X * Free a list header, along with all of its list elements.
X */
Xvoid
Xlistfree(lp)
X    register LIST *lp;
X{
X    register LISTELEM *ep;
X
X    while (lp->l_count-- > 0) {
X        ep = lp->l_first;
X        lp->l_first = ep->e_next;
X        elemfree(ep);
X    }
X    freeitem(&headerfreelist, (FREEITEM *) lp);
X}
X
X
X/*
X * Print out a list along with the specified number of its elements.
X * The elements are printed out in shortened form.
X */
Xvoid
Xlistprint(lp, maxprint)
X    LIST *lp;
X    long maxprint;
X{
X    long count;
X    long index;
X    LISTELEM *ep;
X
X    if (maxprint > lp->l_count)
X        maxprint = lp->l_count;
X    count = 0;
X    ep = lp->l_first;
X    index = lp->l_count;
X    while (index-- > 0) {
X        if ((ep->e_value.v_type != V_NUM) ||
X            (!qiszero(ep->e_value.v_num)))
X                count++;
X        ep = ep->e_next;
X    }
X    if (maxprint > 0)
X        math_str("\n");
X    math_fmt("list (%ld element%s, %ld nonzero)", lp->l_count,
X        ((lp->l_count == 1) ? "" : "s"), count);
X    if (maxprint <= 0)
X        return;
X
X    /*
X     * Walk through the first few list elements, printing their
X     * value in short and unambiguous format.
X     */
X    math_str(":\n");
X    ep = lp->l_first;
X    for (index = 0; index < maxprint; index++) {
X        math_fmt("  [[%ld]] = ", index);
X        printvalue(&ep->e_value, PRINT_SHORT | PRINT_UNAMBIG);
X        math_str("\n");
X        ep = ep->e_next;
X    }
X    if (maxprint < lp->l_count)
X        math_str("  ...\n");
X}
X
X/* END CODE */
END_OF_FILE
if test 11116 -ne `wc -c <'listfunc.c'`; then
    echo shar: \"'listfunc.c'\" unpacked with wrong size!
fi
# end of 'listfunc.c'
fi
if test -f 'qmod.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'qmod.c'\"
else
echo shar: Extracting \"'qmod.c'\" \(10845 characters\)
sed "s/^X//" >'qmod.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 * Modular arithmetic routines for normal numbers, and also using
X * the faster REDC algorithm.
X */
X
X#include <stdio.h>
X#include "math.h"
X
X
X/*
X * Structure used for caching REDC information.
X */
Xtypedef struct    {
X    NUMBER    *num;        /* modulus being cached */
X    REDC    *redc;        /* REDC information for modulus */
X    long    age;        /* age counter for reallocation */
X} REDC_CACHE;
X
X
Xstatic long redc_age;            /* current age counter */
Xstatic REDC_CACHE redc_cache[MAXREDC];    /* cached REDC info */
X
X
Xstatic REDC *qfindredc();
X
X
X/*
X * Return the remainder when one number is divided by another.
X * The second argument cannot be negative.  The result is normalized
X * to lie in the range 0 to q2, and works for fractional values as:
X *    qmod(q1, q2) = q1 - (int(q1 / q2) * q2).
X */
XNUMBER *
Xqmod(q1, q2)
X    register NUMBER *q1, *q2;
X{
X    ZVALUE res;
X    NUMBER *q, *tmp;
X
X    if (qisneg(q2) || qiszero(q2))
X        error("Non-positive modulus");
X    if (qisint(q1) && qisint(q2)) {        /* easy case */
X        zmod(q1->num, q2->num, &res);
X        if (iszero(res)) {
X            freeh(res.v);
X            return qlink(&_qzero_);
X        }
X        if (isone(res)) {
X            freeh(res.v);
X            return qlink(&_qone_);
X        }
X        q = qalloc();
X        q->num = res;
X        return q;
X    }
X    q = qquo(q1, q2);
X    tmp = qmul(q, q2);
X    qfree(q);
X    q = qsub(q1, tmp);
X    qfree(tmp);
X    if (qisneg(q)) {
X        tmp = qadd(q2, q);
X        qfree(q);
X        q = tmp;
X    }
X    return q;
X}
X
X
X/*
X * Given two numbers (a and b), calculate the quotient (c) and remainder (d)
X * when a is divided by b.  This is defined so 0 <= d < b, and c is integral,
X * and a = b * c + d.  The results are returned indirectly through pointers.
X * This works for integers or fractions.  Returns whether or not there is a
X * remainder.  Examples:
X *    qquomod(11, 4, &x, &y) sets x to 2, y to 3, and returns TRUE.
X *    qquomod(-7, 3, &x, &y) sets x to -3, y to 2, and returns TRUE.
X */
XBOOL
Xqquomod(q1, q2, retqdiv, retqmod)
X    NUMBER *q1, *q2;        /* numbers to do quotient with */
X    NUMBER **retqdiv;        /* returned quotient */
X    NUMBER **retqmod;        /* returned modulo */
X{
X    NUMBER *qq, *qm, *tmp;
X
X    if (qisneg(q2) || qiszero(q2))
X        error("Non-positive modulus");
X
X    if (qisint(q1) && qisint(q2)) {        /* integer case */
X        qq = qalloc();
X        qm = qalloc();
X        zdiv(q1->num, q2->num, &qq->num, &qm->num);
X        if (!qisneg(q1) || qiszero(qm)) {
X            *retqdiv = qq;
X            *retqmod = qm;
X            return !qiszero(qm);
X        }
X
X        /*
X         * Need to fix up negative results.
X         */
X        tmp = qdec(qq);
X        qfree(qq);
X        qq = tmp;
X        tmp = qsub(q2, qm);
X        qfree(qm);
X        qm = tmp;
X        *retqdiv = qq;
X        *retqmod = qm;
X        return TRUE;
X    }
X
X    /*
X     * Fractional case.
X     */
X    qq = qquo(q1, q2);
X    tmp = qmul(qq, q2);
X    qm = qsub(q1, tmp);
X    qfree(tmp);
X    if (qisneg(qm)) {
X        tmp = qadd(qm, q2);
X        qfree(qm);
X        qm = tmp;
X        tmp = qdec(qq);
X        qfree(qq);
X        qq = tmp;
X    }
X    *retqdiv = qq;
X    *retqmod = qm;
X    return !qiszero(qm);
X}
X
X
X#if 0
X/*
X * Return the product of two integers modulo a third integer.
X * The result is in the range 0 to q3 - 1 inclusive.
X *    q4 = (q1 * q2) mod q3.
X */
XNUMBER *
Xqmulmod(q1, q2, q3)
X    NUMBER *q1, *q2, *q3;
X{
X    NUMBER *q;
X
X    if (qisneg(q3) || qiszero(q3))
X        error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        error("Non-integers for qmulmod");
X    if (qiszero(q1) || qiszero(q2) || qisunit(q3))
X        return qlink(&_qzero_);
X    q = qalloc();
X    zmulmod(q1->num, q2->num, q3->num, &q->num);
X    return q;
X}
X
X
X/*
X * Return the square of an integer modulo another integer.
X * The result is in the range 0 to q2 - 1 inclusive.
X *    q2 = (q1^2) mod q2.
X */
XNUMBER *
Xqsquaremod(q1, q2)
X    NUMBER *q1, *q2;
X{
X    NUMBER *q;
X
X    if (qisneg(q2) || qiszero(q2))
X        error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for qsquaremod");
X    if (qiszero(q1) || qisunit(q2))
X        return qlink(&_qzero_);
X    if (qisunit(q1))
X        return qlink(&_qone_);
X    q = qalloc();
X    zsquaremod(q1->num, q2->num, &q->num);
X    return q;
X}
X
X
X/*
X * Return the sum of two integers modulo a third integer.
X * The result is in the range 0 to q3 - 1 inclusive.
X *    q4 = (q1 + q2) mod q3.
X */
XNUMBER *
Xqaddmod(q1, q2, q3)
X    NUMBER *q1, *q2, *q3;
X{
X    NUMBER *q;
X
X    if (qisneg(q3) || qiszero(q3))
X        error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        error("Non-integers for qaddmod");
X    q = qalloc();
X    zaddmod(q1->num, q2->num, q3->num, &q->num);
X    return q;
X}
X
X
X/*
X * Return the difference of two integers modulo a third integer.
X * The result is in the range 0 to q3 - 1 inclusive.
X *    q4 = (q1 - q2) mod q3.
X */
XNUMBER *
Xqsubmod(q1, q2, q3)
X    NUMBER *q1, *q2, *q3;
X{
X    NUMBER *q;
X
X    if (qisneg(q3) || qiszero(q3))
X        error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        error("Non-integers for qsubmod");
X    if (q1 == q2)
X        return qlink(&_qzero_);
X    q = qalloc();
X    zsubmod(q1->num, q2->num, q3->num, &q->num);
X    return q;
X}
X
X
X/*
X * Return the negative of an integer modulo another integer.
X * The result is in the range 0 to q2 - 1 inclusive.
X *    q2 = (-q1) mod q2.
X */
XNUMBER *
Xqnegmod(q1, q2)
X    NUMBER *q1, *q2;
X{
X    NUMBER *q;
X
X    if (qisneg(q2) || qiszero(q2))
X        error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for qnegmod");
X    if (qiszero(q1) || qisunit(q2))
X        return qlink(&_qzero_);
X    q = qalloc();
X    znegmod(q1->num, q2->num, &q->num);
X    return q;
X}
X#endif
X
X
X/*
X * Return the integer congruent to an integer whose absolute value is smallest.
X * This is a unique integer in the range int((q2-1)/2 to int(q2/2), inclusive.
X * For example, for a modulus of 7, returned values are [-3, 3], and for a
X * modulus of 8, returned values are [-3, 4].
X */
XNUMBER *
Xqminmod(q1, q2)
X    NUMBER *q1, *q2;
X{
X    NUMBER *q;
X
X    if (qisneg(q2) || qiszero(q2))
X        error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2))
X        error("Non-integers for qminmod");
X    if (qiszero(q1) || (q1->num.len < q2->num.len - 1))
X        return qlink(q1);
X    q = qalloc();
X    zminmod(q1->num, q2->num, &q->num);
X    return q;
X}
X
X
X/*
X * Return whether or not two integers are congruent modulo a third integer.
X * Returns TRUE if the numbers are not congruent, and FALSE if they are.
X */
XBOOL
Xqcmpmod(q1, q2, q3)
X    NUMBER *q1, *q2, *q3;
X{
X    if (qisneg(q3) || qiszero(q3))
X        error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        error("Non-integers for qcmpmod");
X    if (q1 == q2)
X        return FALSE;
X    return zcmpmod(q1->num, q2->num, q3->num);
X}
X
X
X/*
X * Convert an integer into REDC format for use in faster modular arithmetic.
X * The number can be negative or out of modulus range.
X */
XNUMBER *
Xqredcin(q1, q2)
X    NUMBER *q1;        /* number to convert into REDC format */
X    NUMBER *q2;        /* modulus */
X{
X    REDC *rp;        /* REDC information */
X    NUMBER *r;        /* result */
X
X    if (qisfrac(q1))
X        error("Non-integer for qredcin");
X    rp = qfindredc(q2);
X    if (qiszero(q1))
X        return qlink(&_qzero_);
X    r = qalloc();
X    zredcencode(rp, q1->num, &r->num);
X    return r;
X}
X
X
X/*
X * Convert a REDC format number back into a normal integer.
X * The resulting number is in the range 0 to the modulus - 1.
X */
XNUMBER *
Xqredcout(q1, q2)
X    NUMBER *q1;        /* number to convert out of REDC format */
X    NUMBER *q2;        /* modulus */
X{
X    REDC *rp;        /* REDC information */
X    NUMBER *r;        /* result */
X
X    if (qisfrac(q1) || qisneg(q1))
X        error("Non-positive integer required for qredcout");
X    rp = qfindredc(q2);
X    if (qiszero(q1))
X        return qlink(&_qzero_);
X    r = qalloc();
X    zredcdecode(rp, q1->num, &r->num);
X    if (isunit(r->num)) {
X        qfree(r);
X        r = qlink(&_qone_);
X    }
X    return r;
X}
X
X
X/*
X * Multiply two REDC format numbers together producing a REDC format result.
X * This multiplication is done modulo the specified modulus.
X */
XNUMBER *
Xqredcmul(q1, q2, q3)
X    NUMBER *q1, *q2;    /* REDC numbers to be multiplied */
X    NUMBER *q3;        /* modulus */
X{
X    REDC *rp;        /* REDC information */
X    NUMBER *r;        /* result */
X
X    if (qisfrac(q1) || qisneg(q1) || qisfrac(q2) || qisneg(q2))
X        error("Non-positive integers required for qredcmul");
X    rp = qfindredc(q3);
X    if (qiszero(q1) || qiszero(q2))
X        return qlink(&_qzero_);
X    r = qalloc();
X    zredcmul(rp, q1->num, q2->num, &r->num);
X    return r;
X}
X
X
X/*
X * Square a REDC format number to produce a REDC format result.
X * This squaring is done modulo the specified modulus.
X */
XNUMBER *
Xqredcsquare(q1, q2)
X    NUMBER *q1;        /* REDC number to be squared */
X    NUMBER *q2;        /* modulus */
X{
X    REDC *rp;        /* REDC information */
X    NUMBER *r;        /* result */
X
X    if (qisfrac(q1) || qisneg(q1))
X        error("Non-positive integer required for qredcsquare");
X    rp = qfindredc(q2);
X    if (qiszero(q1))
X        return qlink(&_qzero_);
X    r = qalloc();
X    zredcsquare(rp, q1->num, &r->num);
X    return r;
X}
X
X
X/*
X * Raise a REDC format number to the indicated power producing a REDC
X * format result.  This is done modulo the specified modulus.  The
X * power to be raised to is a normal number.
X */
XNUMBER *
Xqredcpower(q1, q2, q3)
X    NUMBER *q1;        /* REDC number to be raised */
X    NUMBER *q2;        /* power to be raised to */
X    NUMBER *q3;        /* modulus */
X{
X    REDC *rp;        /* REDC information */
X    NUMBER *r;        /* result */
X
X    if (qisfrac(q1) || qisneg(q1) || qisfrac(q2) || qisneg(q2))
X        error("Non-positive integers required for qredcpower");
X    rp = qfindredc(q3);
X    if (qiszero(q1) || qisunit(q3))
X        return qlink(&_qzero_);
X    r = qalloc();
X    zredcpower(rp, q1->num, q2->num, &r->num);
X    return r;
X}
X
X
X/*
X * Search for and return the REDC information for the specified number.
X * The information is cached into a local table so that future calls
X * for this information will be quick.  If the table fills up, then
X * the oldest cached entry is reused.
X */
Xstatic REDC *
Xqfindredc(q)
X    NUMBER *q;        /* modulus to find REDC information of */
X{
X    register REDC_CACHE *rcp;
X    REDC_CACHE *bestrcp;
X
X    /*
X     * First try for an exact pointer match in the table.
X     */
X    for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
X        if (q == rcp->num) {
X            rcp->age = ++redc_age;
X            return rcp->redc;
X        }
X    }
X
X    /*
X     * Search the table again looking for a value which matches.
X     */
X    for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
X        if (rcp->age && (qcmp(q, rcp->num) == 0)) {
X            rcp->age = ++redc_age;
X            return rcp->redc;
X        }
X    }
X
X    /*
X     * Must invalidate an existing entry in the table.
X     * Find the oldest (or first unused) entry.
X     * But first make sure the modulus will be reasonable.
X     */
X    if (qisfrac(q) || qiseven(q) || qisneg(q))
X        error("REDC modulus must be positive odd integer");
X
X    bestrcp = NULL;
X    for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
X        if ((bestrcp == NULL) || (rcp->age < bestrcp->age))
X            bestrcp = rcp;
X    }
X
X    /*
X     * Found the best entry.
X     * Free the old information for the entry if necessary,
X     * then initialize it.
X     */
X    rcp = bestrcp;
X    if (rcp->age) {
X        rcp->age = 0;
X        qfree(rcp->num);
X        zredcfree(rcp->redc);
X    }
X
X    rcp->redc = zredcalloc(q->num);
X    rcp->num = qlink(q);
X    rcp->age = ++redc_age;
X    return rcp->redc;
X}
X
X/* END CODE */
END_OF_FILE
if test 10845 -ne `wc -c <'qmod.c'`; then
    echo shar: \"'qmod.c'\" unpacked with wrong size!
fi
# end of 'qmod.c'
fi
echo shar: End of archive 5 \(of 21\).
cp /dev/null ark5isdone
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