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

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

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

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

#!/bin/sh
# this is part 9 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/qmath.c continued
#
CurArch=9
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
     exit 1; fi
( read Scheck
  if test "$Scheck" != $CurArch
  then echo "Please unpack part $Scheck next!"
       exit 1;
  else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file calc2.9.0/qmath.c"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/qmath.c
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 (!zisunit(num))
X        zsquare(num, &q->num);
X    if (!zisunit(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        math_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        math_error("Very large scale value");
X    numshift = zisodd(q->num) ? 0 : zlowbit(q->num);
X    denshift = zisodd(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 (q1 == q2)
X        return qlink(q1);
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 (q1 == q2)
X        return qlink(q1);
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        math_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        math_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 (ziszero(res)) {
X        zfree(res);
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        math_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 (ziszero(res)) {
X        zfree(res);
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 (zisunit(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 * Determine whether or not one number exactly divides another one.
X * Returns TRUE if the first number is an integer multiple of the second one.
X */
XBOOL
Xqdivides(q1, q2)
X    NUMBER *q1, *q2;
X{
X    if (qiszero(q1))
X        return TRUE;
X    if (qisint(q1) && qisint(q2)) {
X        if (qisunit(q2))
X            return TRUE;
X        return zdivides(q1->num, q2->num);
X    }
X    return zdivides(q1->num, q2->num) && zdivides(q2->den, q1->den);
X}
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 (zisunit(q2->den))
X        z1 = q1->num;
X    else if (zisone(q1->num))
X        z1 = q2->den;
X    else {
X        z1f = 1;
X        zmul(q1->num, q2->den, &z1);
X    }
X    if (zisunit(q1->den))
X        z2 = q2->num;
X    else if (zisone(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        zfree(z1);
X    if (z2f)
X        zfree(z2);
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 (zisunit(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            math_error("Not enough memory");
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    zfree(q->num);
X    zfree(q->den);
X    a = (union allocNode *) q;
X    a->link = freeNum;
X    freeNum = a;
X}
X
X/* END CODE */
SHAR_EOF
echo "File calc2.9.0/qmath.c is complete"
chmod 0644 calc2.9.0/qmath.c || echo "restore of calc2.9.0/qmath.c fails"
set `wc -c calc2.9.0/qmath.c`;Sum=$1
if test "$Sum" != "20680"
then echo original size 20680, current size $Sum;fi
echo "x - extracting calc2.9.0/qmath.h (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/qmath.h &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Data structure declarations for extended precision rational arithmetic.
X */
X
X#ifndef    QMATH_H
X#define    QMATH_H
X
X#include "zmath.h"
X
X
X/*
X * Rational arithmetic definitions.
X */
Xtypedef struct {
X    ZVALUE num;        /* numerator (containing sign) */
X    ZVALUE den;        /* denominator (always positive) */
X    long links;        /* number of links to this value */
X} NUMBER;
X
X
X/*
X * Input. output, allocation, and conversion routines.
X */
Xextern NUMBER *qalloc MATH_PROTO((void));
Xextern NUMBER *qcopy MATH_PROTO((NUMBER *q));
Xextern NUMBER *iitoq MATH_PROTO((long i1, long i2));
Xextern NUMBER *atoq MATH_PROTO((char *str));
Xextern NUMBER *itoq MATH_PROTO((long i));
Xextern long qtoi MATH_PROTO((NUMBER *q));
Xextern long qparse MATH_PROTO((char *str, int flags));
Xextern void qfreenum MATH_PROTO((NUMBER *q));
Xextern void qprintnum MATH_PROTO((NUMBER *q, int mode));
Xextern void qprintff MATH_PROTO((NUMBER *q, long width, long precision));
Xextern void qprintfe MATH_PROTO((NUMBER *q, long width, long precision));
Xextern void qprintfr MATH_PROTO((NUMBER *q, long width, BOOL force));
Xextern void qprintfd MATH_PROTO((NUMBER *q, long width));
Xextern void qprintfx MATH_PROTO((NUMBER *q, long width));
Xextern void qprintfb MATH_PROTO((NUMBER *q, long width));
Xextern void qprintfo MATH_PROTO((NUMBER *q, long width));
X
X
X
X/*
X * Basic numeric routines.
X */
Xextern NUMBER *qaddi MATH_PROTO((NUMBER *q, long i));
Xextern NUMBER *qmuli MATH_PROTO((NUMBER *q, long i));
Xextern NUMBER *qdivi MATH_PROTO((NUMBER *q, long i));
Xextern NUMBER *qadd MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qsub MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qmul MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qdiv MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qquo MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qmod MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qmin MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qmax MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qand MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qor MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qxor MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qpowermod MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *q3));
Xextern NUMBER *qpowi MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qsquare MATH_PROTO((NUMBER *q));
Xextern NUMBER *qneg MATH_PROTO((NUMBER *q));
Xextern NUMBER *qsign MATH_PROTO((NUMBER *q));
Xextern NUMBER *qint MATH_PROTO((NUMBER *q));
Xextern NUMBER *qfrac MATH_PROTO((NUMBER *q));
Xextern NUMBER *qnum MATH_PROTO((NUMBER *q));
Xextern NUMBER *qden MATH_PROTO((NUMBER *q));
Xextern NUMBER *qinv MATH_PROTO((NUMBER *q));
Xextern NUMBER *qabs MATH_PROTO((NUMBER *q));
Xextern NUMBER *qinc MATH_PROTO((NUMBER *q));
Xextern NUMBER *qdec MATH_PROTO((NUMBER *q));
Xextern NUMBER *qshift MATH_PROTO((NUMBER *q, long n));
Xextern NUMBER *qtrunc MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qround MATH_PROTO((NUMBER *q, long places));
Xextern NUMBER *qbtrunc MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qbround MATH_PROTO((NUMBER *q, long places));
Xextern NUMBER *qscale MATH_PROTO((NUMBER *q, long i));
Xextern BOOL qdivides MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern BOOL qcmp MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern BOOL qcmpi MATH_PROTO((NUMBER *q, long i));
Xextern FLAG qrel MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern FLAG qreli MATH_PROTO((NUMBER *q, long i));
Xextern BOOL qisset MATH_PROTO((NUMBER *q, long i));
X
X
X/*
X * More complicated numeric functions.
X */
Xextern NUMBER *qcomb MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qgcd MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qlcm MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qfact MATH_PROTO((NUMBER *q));
Xextern NUMBER *qpfact MATH_PROTO((NUMBER *q));
Xextern NUMBER *qminv MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qfacrem MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qperm MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qgcdrem MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qlowfactor MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qfib MATH_PROTO((NUMBER *q));
Xextern NUMBER *qcfappr MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qisqrt MATH_PROTO((NUMBER *q));
Xextern NUMBER *qjacobi MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qiroot MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qbappr MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qlcmfact MATH_PROTO((NUMBER *q));
Xextern NUMBER *qminmod MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qredcin MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qredcout MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qredcmul MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *q3));
Xextern NUMBER *qredcsquare MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qredcpower MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *q3));
Xextern BOOL qprimetest MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern BOOL qissquare MATH_PROTO((NUMBER *q));
Xextern long qilog2 MATH_PROTO((NUMBER *q));
Xextern long qilog10 MATH_PROTO((NUMBER *q));
Xextern BOOL qcmpmod MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *q3));
Xextern BOOL qquomod MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER **retdiv,
X    NUMBER **retmod));
Xextern FLAG qnear MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *epsilon));
Xextern FLAG qdigit MATH_PROTO((NUMBER *q, long i));
Xextern long qprecision MATH_PROTO((NUMBER *q));
Xextern long qplaces MATH_PROTO((NUMBER *q));
Xextern long qdigits MATH_PROTO((NUMBER *q));
Xextern HASH qhash MATH_PROTO((NUMBER *q));
Xextern void setepsilon MATH_PROTO((NUMBER *q));
X
X#if 0
Xextern NUMBER *qbitvalue MATH_PROTO((long i));
Xextern NUMBER *qmulmod MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *q3));
Xextern NUMBER *qsquaremod MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern NUMBER *qaddmod MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *q3));
Xextern NUMBER *qsubmod MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *q3));
Xextern NUMBER *qreadval MATH_PROTO((FILE *fp));
Xextern NUMBER *qnegmod MATH_PROTO((NUMBER *q1, NUMBER *q2));
Xextern BOOL qbittest MATH_PROTO((NUMBER *q, long i));
Xextern FLAG qtest MATH_PROTO((NUMBER *q));
X#endif
X
X
X/*
X * Transcendental functions.  These all take an epsilon argument to
X * specify the required accuracy of the calculation.
X */
Xextern NUMBER *qsqrt MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qpower MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *epsilon));
Xextern NUMBER *qroot MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *epsilon));
Xextern NUMBER *qcos MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qsin MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qexp MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qln MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qtan MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qacos MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qasin MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qatan MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qatan2 MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *epsilon));
Xextern NUMBER *qhypot MATH_PROTO((NUMBER *q1, NUMBER *q2, NUMBER *epsilon));
Xextern NUMBER *qcosh MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qsinh MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qtanh MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qacosh MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qasinh MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qatanh MATH_PROTO((NUMBER *q, NUMBER *epsilon));
Xextern NUMBER *qlegtoleg MATH_PROTO((NUMBER *q, NUMBER *epsilon, BOOL wantneg));
Xextern NUMBER *qpi MATH_PROTO((NUMBER *epsilon));
X
X
X/*
X * macro expansions to speed this thing up
X */
X#define qiszero(q)    (ziszero((q)->num))
X#define qisneg(q)    (zisneg((q)->num))
X#define qispos(q)    (zispos((q)->num))
X#define qisint(q)    (zisunit((q)->den))
X#define qisfrac(q)    (!zisunit((q)->den))
X#define qisunit(q)    (zisunit((q)->num) && zisunit((q)->den))
X#define qisone(q)    (zisone((q)->num) && zisunit((q)->den))
X#define qisnegone(q)    (zisnegone((q)->num) && zisunit((q)->den))
X#define qistwo(q)    (zistwo((q)->num) && zisunit((q)->den))
X#define qiseven(q)    (zisunit((q)->den) && ziseven((q)->num))
X#define qisodd(q)    (zisunit((q)->den) && zisodd((q)->num))
X#define qistwopower(q)    (zisunit((q)->den) && zistwopower((q)->num))
X
X#define qhighbit(q)    (zhighbit((q)->num))
X#define qlowbit(q)    (zlowbit((q)->num))
X#define qdivcount(q1, q2)    (zdivcount((q1)->num, (q2)->num))
X#define qilog(q1, q2)    (zlog((q1)->num, (q2)->num))
X#define qlink(q)    ((q)->links++, (q))
X
X#define qfree(q)    {if (--((q)->links) <= 0) qfreenum(q);}
X
X
X/*
X * Flags for qparse calls
X */
X#define QPF_SLASH    0x1    /* allow slash for fractional number */
X#define QPF_IMAG    0x2    /* allow trailing 'i' for imaginary number */
X
X
X#ifdef VARARGS
Xextern void qprintf();
X#else
Xextern void qprintf MATH_PROTO((char *, ...));
X#endif
X
X
X/*
X * constants used often by the arithmetic routines
X */
Xextern NUMBER _qzero_, _qone_, _qnegone_, _qonehalf_;
Xextern BOOL _sinisneg_;        /* whether sin(x) < 0 (set by cos(x)) */
Xextern long _epsilonprec_;    /* binary precision of epsilon */
Xextern NUMBER *_epsilon_;    /* default error for real functions */
Xextern long _outdigits_;    /* current output digits for float or exp */
Xextern int _outmode_;        /* current output mode */
X
X#endif
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/qmath.h || echo "restore of calc2.9.0/qmath.h fails"
set `wc -c calc2.9.0/qmath.h`;Sum=$1
if test "$Sum" != "9457"
then echo original size 9457, current size $Sum;fi
echo "x - extracting calc2.9.0/qmod.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/qmod.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Modular arithmetic routines for normal numbers, and also using
X * the faster REDC algorithm.
X */
X
X#include "qmath.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 MATH_PROTO((NUMBER *q));
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        math_error("Non-positive modulus");
X    if (qisint(q1) && qisint(q2)) {        /* easy case */
X        zmod(q1->num, q2->num, &res);
X        if (ziszero(res)) {
X            zfree(res);
X            return qlink(&_qzero_);
X        }
X        if (zisone(res)) {
X            zfree(res);
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        math_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        math_error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        math_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        math_error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2))
X        math_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        math_error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        math_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        math_error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        math_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        math_error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2))
X        math_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        math_error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2))
X        math_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        math_error("Non-positive modulus");
X    if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X        math_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        math_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        math_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 (zisunit(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        math_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        math_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        math_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        math_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 */
SHAR_EOF
chmod 0644 calc2.9.0/qmod.c || echo "restore of calc2.9.0/qmod.c fails"
set `wc -c calc2.9.0/qmod.c`;Sum=$1
if test "$Sum" != "10958"
then echo original size 10958, current size $Sum;fi
echo "x - extracting calc2.9.0/qtrans.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/qtrans.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Transcendental functions for real numbers.
X * These are sin, cos, exp, ln, power, cosh, sinh.
X */
X
X#include "qmath.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        math_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        math_error("Illegal epsilon value for sine");
X    if (qiszero(q))
X        return qlink(q);
X    epsilon2 = qsquare(epsilon);
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        math_error("Illegal epsilon value for tangent");
X    if (qiszero(q))
X        return qlink(q);
X    epsilon2 = qsquare(epsilon);
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        math_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        math_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, *tmp3, *epsilon2;
X
X    if (qisneg(epsilon) || qiszero(epsilon))
X        math_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        math_error("Argument too large for acos");
X    /*
X     * Calculate the result using the formula:
X     *    acos(x) = asin(sqrt(1 - x^2)).
X     * The formula is only good for positive x, so we must fix up the
X     * result for negative values.
X     */
X    epsilon2 = qscale(epsilon, -8L);
X    tmp1 = qlegtoleg(q, epsilon2, FALSE);
X    qfree(epsilon2);
X    tmp2 = qasin(tmp1, epsilon);
X    qfree(tmp1);
X    if (!qisneg(q))
X        return tmp2;
X    /*
X     * For negative values, we need to subtract the asin from pi.
X     */
X    tmp1 = qpi(epsilon);
X    tmp3 = qsub(tmp1, tmp2);
X    qfree(tmp1);
X    qfree(tmp2);
X    tmp1 = qbround(tmp3, qprecision(epsilon) + 1);
X    qfree(tmp3);
X    return tmp1;
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        math_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        math_error("Illegal epsilon value for atan2");
X    if (qiszero(qy) && qiszero(qx))
X        math_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        math_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        zfree(comb);
X        zmuli(tmp1, t * (2 * N - 1), &comb);
X        zfree(tmp1);
X        zsquare(comb, &tmp1);
X        zmul(comb, tmp1, &tmp2);
X        zfree(tmp1);
X        zmuli(tmp2, 42 * N + 5, &tmp1);
X        zfree(tmp2);
X        zshift(sum, 12L, &tmp2);
X        zfree(sum);
X        zadd(tmp1, tmp2, &sum);
X        t = zhighbit(tmp1);
X        zfree(tmp1);
X        zfree(tmp2);
X        shift += 12;
X    } while ((shift - t) < bits);
X    qtmp.num = _one_;
X    qtmp.den = sum;
X    t1 = qscale(&qtmp, shift);
X    zfree(sum);
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        math_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            math_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, *maxr;
X    long shift, bits, bits2;
X    int j, k;
X    FULL n;
X    BOOL neg;
X
X    if (qisneg(q) || qiszero(q))
X        math_error("log of non-positive number");
X    if (qisneg(epsilon) || qiszero(epsilon))
X        math_error("Illegal epsilon for ln");
X    if (qisone(q))
X        return qlink(&_qzero_);
X    /*
X     * If the number is less than one, invert it and remember that
X     * the result is to be negative.
X     */
X    neg = FALSE;
X    if (zrel(q->num, q->den) < 0) {
X        neg = TRUE;
X        q = qinv(q);
X    } else
X        q = qlink(q);
X    j = 16;
X    k = zhighbit(q->num) - zhighbit(q->den) + 1;
X    while (k >>= 1)
X        j++;
X    epsilon2 = qscale(epsilon, -j);
X    bits = qprecision(epsilon) + 1;
X    bits2 = qprecision(epsilon2) + 5;
X    /*
X     * By repeated square-roots scale number down to a value close
X     * to 1 so that Taylor series to be used will converge rapidly.
X     * The effect of scaling will be reversed by a later shift.
X     */
X    maxr = iitoq(65537L, 65536L);
X    shift = 1;
X    while (qrel(q, maxr) > 0) {
X        tmp1 = qsqrt(q, epsilon2);
X        qfree(q);
X        q = tmp1;
X        shift++;
X    }
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    /*
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    /*
X     * Calculate the final result by multiplying by the proper power
X     * of two to undo the square roots done at the top, and possibly
X     * negating the result.
X     */
X    tmp1 = qscale(sum, shift);
X    qfree(sum);
X    sum = qbround(tmp1, bits);
X    qfree(tmp1);
X    if (neg) {
X        tmp1 = qneg(sum);
X        qfree(sum);
X        sum = tmp1;
X    }
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        math_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 (ziseven(q2->num))
X            math_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        math_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            math_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        math_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        math_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        math_error("Illegal epsilon value for acosh");
X    if (qisone(q))
X        return qlink(&_qzero_);
X    if (qreli(q, 1L) < 0)
X        math_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        math_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        math_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        math_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 */
SHAR_EOF
chmod 0644 calc2.9.0/qtrans.c || echo "restore of calc2.9.0/qtrans.c fails"
set `wc -c calc2.9.0/qtrans.c`;Sum=$1
if test "$Sum" != "21232"
then echo original size 21232, current size $Sum;fi
echo "x - extracting calc2.9.0/stdarg.h (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/stdarg.h &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X */
X
X#include "args.h"
X
X#ifndef STDARG_H
X#define STDARG_H
X
X#ifdef VARARGS
X
X#include <varargs.h>
X
X#else /*VARARG*/
X
X#ifdef STDARG
X
X#if defined(__STDC__)
X#include <stdarg.h>
X#else
X%;%;%  BOGUS DEFINE!  Must use ANSI C when STDARG is defined  %;%;%
X#endif
X
X#else /*STDARG*/
X
X/*
X * SIMULATE_STDARG
X *
X * WARNING: This type of stdarg makes assumptions about the stack
X *         that may not be true on your system.  You may want to
X *        define STDARG (if using ANSI C) or VARARGS.
X */
X
Xtypedef char *va_list;
X#define va_start(ap,parmn) (void)((ap) = (char*)(&(parmn) + 1))
X#define va_end(ap) (void)((ap) = 0)
X#define va_arg(ap, type) \
X    (((type*)((ap) = ((ap) + sizeof(type))))[-1])
X
X#endif /*STDARG*/
X#endif /*VARARG*/
X
X#endif
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/stdarg.h || echo "restore of calc2.9.0/stdarg.h fails"
set `wc -c calc2.9.0/stdarg.h`;Sum=$1
if test "$Sum" != "904"
then echo original size 904, current size $Sum;fi
echo "x - extracting calc2.9.0/string.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/string.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * String list routines.
X */
X
X#include "calc.h"
X#include "string.h"
X
X#define STR_TABLECHUNK    100    /* how often to reallocate string table */
X#define STR_CHUNK    2000    /* size of string storage allocation */
X#define STR_UNIQUE    100    /* size of string to allocate separately */
X
X
Xstatic char *chartable;        /* single character string table */
X
Xstatic struct {
X    long l_count;        /* count of strings in table */
X    long l_maxcount;    /* maximum strings storable in table */
X    long l_avail;        /* characters available in current string */
X    char *l_alloc;        /* next available string storage */
X    char **l_table;        /* current string table */
X} literals;
X
X
X/*
X * Initialize or reinitialize a string header for use.
X */
Xvoid
Xinitstr(hp)
X    register STRINGHEAD *hp;    /* structure to be inited */
X{
X    if (hp->h_list == NULL) {
X        hp->h_list = (char *)malloc(2000);
X        hp->h_avail = 2000;
X        hp->h_used = 0;
X    }
X    hp->h_avail += hp->h_used;
X    hp->h_used = 0;
X    hp->h_count = 0;
X    hp->h_list[0] = '\0';
X    hp->h_list[1] = '\0';
X}
X
X
X/*
X * Copy a string to the end of a list of strings, and return the address
X * of the copied string.  Returns NULL if the string could not be copied.
X * No checks are made to see if the string is already in the list.
X * The string cannot be null or have imbedded nulls.
X */
Xchar *
Xaddstr(hp, str)
X    register STRINGHEAD *hp;    /* header of string storage */
X    char *str;        /* string to be added */
X{
X    char *retstr;        /* returned string pointer */
X    char *list;        /* string list */
X    long newsize;        /* new size of string list */
X    long len;        /* length of current string */
X
X    if ((str == NULL) || (*str == '\0'))
X        return NULL;
X    len = strlen(str) + 1;
X    if (hp->h_avail <= len) {
X        newsize = len + 2000 + hp->h_used + hp->h_avail;
X        list = (char *)realloc(hp->h_list, newsize);
X        if (list == NULL)
X            return NULL;
X        hp->h_list = list;
X        hp->h_avail = newsize - hp->h_used;
X    }
X    retstr = hp->h_list + hp->h_used;
X    hp->h_used += len;
X    hp->h_avail -= len;
X    hp->h_count++;
X    strcpy(retstr, str);
X    retstr[len] = '\0';
X    return retstr;
X}
X
X
X/*
X * Return a null-terminated string which consists of a single character.
X * The table is initialized on the first call.
X */
Xchar *
Xcharstr(ch)
X{
X    char *cp;
X    int i;
X
X    if (chartable == NULL) {
X        cp = (char *)malloc(512);
X        if (cp == NULL)
X            math_error("Cannot allocate character table");
X        for (i = 0; i < 256; i++) {
X            *cp++ = (char)i;
X            *cp++ = '\0';
X        }
X        chartable = cp - 512;
X    }
X    return &chartable[(ch & 0xff) * 2];
X}
X
X
X/*
X * Find a string with the specified name and return its number in the
X * string list.  The first string is numbered zero.  Minus one is returned
X * if the string is not found.
X */
Xlong
Xfindstr(hp, str)
X    STRINGHEAD *hp;        /* header of string storage */
X    register char *str;    /* string to be added */
X{
X    register char *test;    /* string being tested */
X    long len;        /* length of string being found */
X    long testlen;        /* length of test string */
X    long index;        /* index of string */
X
X    if ((hp->h_count <= 0) || (str == NULL))
X        return -1;
X    len = strlen(str);
X    test = hp->h_list;
X    index = 0;
X    while (*test) {
X        testlen = strlen(test);
X        if ((testlen == len) && (*test == *str) && (strcmp(test, str) == 0))
X            return index;
X        test += (testlen + 1);
X        index++;
X    }
X    return -1;
X}
X
X
X/*
X * Return the name of a string with the given index.
X * If the index is illegal, a pointer to an empty string is returned.
X */
Xchar *
Xnamestr(hp, n)
X    STRINGHEAD *hp;        /* header of string storage */
X    long n;
X{
X    register char *str;    /* current string */
X
X    if ((unsigned long)n >= hp->h_count)
X        return "";
X    str = hp->h_list;
X    while (*str) {
X        if (--n < 0)
X            return str;
X        str += (strlen(str) + 1);
X    }
X    return "";
X}
X
X
X/*
X * Useful routine to return the index of one string within another one
X * which has the format:  "str1\0str2\0str3\0...strn\0\0".  Index starts
X * at one for the first string.  Returns zero if the string being checked
X * is not contained in the formatted string.
X */
Xlong
Xstringindex(format, test)
X    register char *format;    /* string formatted into substrings */
X    char *test;        /* string to be found in formatted string */
X{
X    long index;        /* found index */
X    long len;        /* length of current piece of string */
X    long testlen;        /* length of test string */
SHAR_EOF
echo "End of part 9"
echo "File calc2.9.0/string.c is continued in part 10"
echo "10" > s2_seq_.tmp
exit 0