Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/unix/volume27/calc-2.9.0/part11

Text File | 1993-12-07 | 58.1 KB | 2,448 lines

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

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

#!/bin/sh
# this is part 11 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/value.c continued
#
CurArch=11
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/value.c"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/value.c
X * Raise one value to another value's power, within the specified error.
X * Result is placed in the indicated location.
X */
Xvoid
Xpowervalue(v1, v2, v3, vres)
X    VALUE *v1, *v2, *v3, *vres;
X{
X    NUMBER *epsilon;
X    COMPLEX *c, ctmp;
X
X    vres->v_type = V_NULL;
X    if (v3->v_type != V_NUM)
X        math_error("Non-real epsilon value for power");
X    epsilon = v3->v_num;
X    if (qisneg(epsilon) || qiszero(epsilon))
X        math_error("Non-positive epsilon value for power");
X    switch (TWOVAL(v1->v_type, v2->v_type)) {
X        case TWOVAL(V_NUM, V_NUM):
X            vres->v_num = qpower(v1->v_num, v2->v_num, epsilon);
X            vres->v_type = V_NUM;
X            return;
X        case TWOVAL(V_NUM, V_COM):
X            ctmp.real = v1->v_num;
X            ctmp.imag = &_qzero_;
X            ctmp.links = 1;
X            vres->v_com = cpower(&ctmp, v2->v_com, epsilon);
X            break;
X        case TWOVAL(V_COM, V_NUM):
X            ctmp.real = v2->v_num;
X            ctmp.imag = &_qzero_;
X            ctmp.links = 1;
X            vres->v_com = cpower(v1->v_com, &ctmp, epsilon);
X            break;
X        case TWOVAL(V_COM, V_COM):
X            vres->v_com = cpower(v1->v_com, v2->v_com, epsilon);
X            break;
X        default:
X            math_error("Illegal value for raising to power");
X    }
X    /*
X     * Here for any complex result.
X     */
X    vres->v_type = V_COM;
X    c = vres->v_com;
X    if (!cisreal(c))
X        return;
X    vres->v_num = qlink(c->real);
X    vres->v_type = V_NUM;
X    comfree(c);
X}
X
X
X/*
X * Divide one arbitrary value by another one.
X * Result is placed in the indicated location.
X */
Xvoid
Xdivvalue(v1, v2, vres)
X    VALUE *v1, *v2, *vres;
X{
X    COMPLEX *c;
X    COMPLEX ctmp;
X    VALUE tmpval;
X
X    vres->v_type = V_NULL;
X    switch (TWOVAL(v1->v_type, v2->v_type)) {
X        case TWOVAL(V_NUM, V_NUM):
X            vres->v_num = qdiv(v1->v_num, v2->v_num);
X            vres->v_type = V_NUM;
X            return;
X        case TWOVAL(V_COM, V_NUM):
X            vres->v_com = cdivq(v1->v_com, v2->v_num);
X            vres->v_type = V_COM;
X            return;
X        case TWOVAL(V_NUM, V_COM):
X            if (qiszero(v1->v_num)) {
X                vres->v_num = qlink(&_qzero_);
X                vres->v_type = V_NUM;
X                return;
X            }
X            ctmp.real = v1->v_num;
X            ctmp.imag = &_qzero_;
X            ctmp.links = 1;
X            vres->v_com = cdiv(&ctmp, v2->v_com);
X            vres->v_type = V_COM;
X            return;
X        case TWOVAL(V_COM, V_COM):
X            vres->v_com = cdiv(v1->v_com, v2->v_com);
X            vres->v_type = V_COM;
X            c = vres->v_com;
X            if (cisreal(c)) {
X                vres->v_num = qlink(c->real);
X                vres->v_type = V_NUM;
X                comfree(c);
X            }
X            return;
X        case TWOVAL(V_MAT, V_NUM):
X        case TWOVAL(V_MAT, V_COM):
X            invertvalue(v2, &tmpval);
X            vres->v_mat = matmulval(v1->v_mat, &tmpval);
X            vres->v_type = V_MAT;
X            freevalue(&tmpval);
X            return;
X        default:
X            if ((v1->v_type != V_OBJ) && (v2->v_type != V_OBJ))
X                math_error("Non-compatible values for divide");
X            *vres = objcall(OBJ_DIV, v1, v2, NULL_VALUE);
X            return;
X    }
X}
X
X
X/*
X * Divide one arbitrary value by another one keeping only the integer part.
X * Result is placed in the indicated location.
X */
Xvoid
Xquovalue(v1, v2, vres)
X    VALUE *v1, *v2, *vres;
X{
X    COMPLEX *c;
X
X    vres->v_type = V_NULL;
X    switch (TWOVAL(v1->v_type, v2->v_type)) {
X        case TWOVAL(V_NUM, V_NUM):
X            vres->v_num = qquo(v1->v_num, v2->v_num);
X            vres->v_type = V_NUM;
X            return;
X        case TWOVAL(V_COM, V_NUM):
X            vres->v_com = cquoq(v1->v_com, v2->v_num);
X            vres->v_type = V_COM;
X            c = vres->v_com;
X            if (cisreal(c)) {
X                vres->v_num = qlink(c->real);
X                vres->v_type = V_NUM;
X                comfree(c);
X            }
X            return;
X        case TWOVAL(V_MAT, V_NUM):
X        case TWOVAL(V_MAT, V_COM):
X            vres->v_mat = matquoval(v1->v_mat, v2);
X            vres->v_type = V_MAT;
X            return;
X        default:
X            if ((v1->v_type != V_OBJ) && (v2->v_type != V_OBJ))
X                math_error("Non-compatible values for quotient");
X            *vres = objcall(OBJ_QUO, v1, v2, NULL_VALUE);
X            return;
X    }
X}
X
X
X/*
X * Divide one arbitrary value by another one keeping only the remainder.
X * Result is placed in the indicated location.
X */
Xvoid
Xmodvalue(v1, v2, vres)
X    VALUE *v1, *v2, *vres;
X{
X    COMPLEX *c;
X
X    vres->v_type = V_NULL;
X    switch (TWOVAL(v1->v_type, v2->v_type)) {
X        case TWOVAL(V_NUM, V_NUM):
X            vres->v_num = qmod(v1->v_num, v2->v_num);
X            vres->v_type = V_NUM;
X            return;
X        case TWOVAL(V_COM, V_NUM):
X            vres->v_com = cmodq(v1->v_com, v2->v_num);
X            vres->v_type = V_COM;
X            c = vres->v_com;
X            if (cisreal(c)) {
X                vres->v_num = qlink(c->real);
X                vres->v_type = V_NUM;
X                comfree(c);
X            }
X            return;
X        case TWOVAL(V_MAT, V_NUM):
X        case TWOVAL(V_MAT, V_COM):
X            vres->v_mat = matmodval(v1->v_mat, v2);
X            vres->v_type = V_MAT;
X            return;
X        default:
X            if ((v1->v_type != V_OBJ) && (v2->v_type != V_OBJ))
X                math_error("Non-compatible values for mod");
X            *vres = objcall(OBJ_MOD, v1, v2, NULL_VALUE);
X            return;
X    }
X}
X
X
X/*
X * Test an arbitrary value to see if it is equal to "zero".
X * The definition of zero varies depending on the value type.  For example,
X * the null string is "zero", and a matrix with zero values is "zero".
X * Returns TRUE if value is not equal to zero.
X */
XBOOL
Xtestvalue(vp)
X    VALUE *vp;
X{
X    VALUE val;
X
X    switch (vp->v_type) {
X        case V_NUM:
X            return !qiszero(vp->v_num);
X        case V_COM:
X            return !ciszero(vp->v_com);
X        case V_STR:
X            return (vp->v_str[0] != '\0');
X        case V_MAT:
X            return mattest(vp->v_mat);
X        case V_LIST:
X            return (vp->v_list->l_count != 0);
X        case V_ASSOC:
X            return (vp->v_assoc->a_count != 0);
X        case V_FILE:
X            return validid(vp->v_file);
X        case V_NULL:
X            return FALSE;
X        case V_OBJ:
X            val = objcall(OBJ_TEST, vp, NULL_VALUE, NULL_VALUE);
X            return (val.v_int != 0);
X        default:
X            return TRUE;
X    }
X}
X
X
X/*
X * Compare two values for equality.
X * Returns TRUE if the two values differ.
X */
XBOOL
Xcomparevalue(v1, v2)
X    VALUE *v1, *v2;
X{
X    int r;
X    VALUE val;
X
X    if ((v1->v_type == V_OBJ) || (v2->v_type == V_OBJ)) {
X        val = objcall(OBJ_CMP, v1, v2, NULL_VALUE);
X        return (val.v_int != 0);
X    }
X    if (v1 == v2)
X        return FALSE;
X    if (v1->v_type != v2->v_type)
X        return TRUE;
X    switch (v1->v_type) {
X        case V_NUM:
X            r = qcmp(v1->v_num, v2->v_num);
X            break;
X        case V_COM:
X            r = ccmp(v1->v_com, v2->v_com);
X            break;
X        case V_STR:
X            r = ((v1->v_str != v2->v_str) &&
X                ((v1->v_str[0] - v2->v_str[0]) ||
X                strcmp(v1->v_str, v2->v_str)));
X            break;
X        case V_MAT:
X            r = matcmp(v1->v_mat, v2->v_mat);
X            break;
X        case V_LIST:
X            r = listcmp(v1->v_list, v2->v_list);
X            break;
X        case V_ASSOC:
X            r = assoccmp(v1->v_assoc, v2->v_assoc);
X            break;
X        case V_NULL:
X            r = FALSE;
X            break;
X        case V_FILE:
X            r = (v1->v_file != v2->v_file);
X            break;
X        default:
X            math_error("Illegal values for comparevalue");
X    }
X    return (r != 0);
X}
X
X
X/*
X * Compare two values for their relative values.
X * Returns minus one if the first value is less than the second one,
X * one if the first value is greater than the second one, and
X * zero if they are equal.
X */
XFLAG
Xrelvalue(v1, v2)
X    VALUE *v1, *v2;
X{
X    int r;
X    VALUE val;
X
X    if ((v1->v_type == V_OBJ) || (v2->v_type == V_OBJ)) {
X        val = objcall(OBJ_REL, v1, v2, NULL_VALUE);
X        return val.v_int;
X    }
X    if (v1 == v2)
X        return 0;
X    if (v1->v_type != v2->v_type)
X        math_error("Relative comparison of differing types");
X    switch (v1->v_type) {
X        case V_NUM:
X            r = qrel(v1->v_num, v2->v_num);
X            break;
X        case V_STR:
X            r = strcmp(v1->v_str, v2->v_str);
X            break;
X        case V_NULL:
X            r = 0;
X            break;
X        default:
X            math_error("Illegal value for relative comparison");
X    }
X    if (r < 0)
X        return -1;
X    return (r != 0);
X}
X
X
X/*
X * Calculate a hash value for a value.
X * The hash does not have to be a perfect one, it is only used for
X * making associations faster.
X */
XHASH
Xhashvalue(vp)
X    VALUE *vp;
X{
X    switch (vp->v_type) {
X        case V_INT:
X            return ((long) vp->v_int);
X            break;
X        case V_NUM:
X            return qhash(vp->v_num);
X            break;
X        case V_COM:
X            return chash(vp->v_com);
X            break;
X        case V_STR:
X            return hashstr(vp->v_str);
X            break;
X        case V_NULL:
X            return 0;
X            break;
X        case V_OBJ:
X            return objhash(vp->v_obj);
X            break;
X        case V_LIST:
X            return listhash(vp->v_list);
X            break;
X        case V_ASSOC:
X            return assochash(vp->v_assoc);
X            break;
X        case V_MAT:
X            return mathash(vp->v_mat);
X            break;
X        case V_FILE:
X            return ((long) vp->v_file);
X            break;
X        default:
X            math_error("Hashing unknown value");
X    }
X    return 0;
X}
X
X
X/*
X * Print the value of a descriptor in one of several formats.
X * If flags contains PRINT_SHORT, then elements of arrays and lists
X * will not be printed.  If flags contains PRINT_UNAMBIG, then quotes
X * are placed around strings and the null value is explicitly printed.
X */
Xvoid
Xprintvalue(vp, flags)
X    VALUE *vp;
X{
X    switch (vp->v_type) {
X        case V_NUM:
X            qprintnum(vp->v_num, MODE_DEFAULT);
X            break;
X        case V_COM:
X            comprint(vp->v_com);
X            break;
X        case V_STR:
X            if (flags & PRINT_UNAMBIG)
X                math_chr('\"');
X            math_str(vp->v_str);
X            if (flags & PRINT_UNAMBIG)
X                math_chr('\"');
X            break;
X        case V_NULL:
X            if (flags & PRINT_UNAMBIG)
X                math_str("NULL");
X            break;
X        case V_OBJ:
X            (void) objcall(OBJ_PRINT, vp, NULL_VALUE, NULL_VALUE);
X            break;
X        case V_LIST:
X            listprint(vp->v_list,
X                ((flags & PRINT_SHORT) ? 0L : maxprint));
X            break;
X        case V_ASSOC:
X            assocprint(vp->v_assoc,
X                ((flags & PRINT_SHORT) ? 0L : maxprint));
X            break;
X        case V_MAT:
X            matprint(vp->v_mat,
X                ((flags & PRINT_SHORT) ? 0L : maxprint));
X            break;
X        case V_FILE:
X            printid(vp->v_file, flags);
X            break;
X        default:
X            math_error("Printing unknown value");
X    }
X}
X
X/* END CODE */
SHAR_EOF
echo "File calc2.9.0/value.c is complete"
chmod 0644 calc2.9.0/value.c || echo "restore of calc2.9.0/value.c fails"
set `wc -c calc2.9.0/value.c`;Sum=$1
if test "$Sum" != "29836"
then echo original size 29836, current size $Sum;fi
echo "x - extracting calc2.9.0/value.h (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/value.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 * Definitions of general values and related routines used by the calculator.
X */
X
X#ifndef    VALUE_H
X#define    VALUE_H
X
X#include "cmath.h"
X
X
X#define MAXDIM        4    /* maximum number of dimensions in matrices */
X#define USUAL_ELEMENTS    4    /* usual number of elements for objects */
X
X
X/*
X * Flags to modify results from the printvalue routine.
X * These flags are OR'd together.
X */
X#define    PRINT_NORMAL    0x00    /* print in normal manner */
X#define    PRINT_SHORT    0x01    /* print in short format (no elements) */
X#define    PRINT_UNAMBIG    0x02    /* print in non-ambiguous manner */
X
X
X/*
X * Definition of values of various types.
X */
Xtypedef struct value VALUE;
Xtypedef struct object OBJECT;
Xtypedef struct matrix MATRIX;
Xtypedef struct list LIST;
Xtypedef    struct assoc ASSOC;
Xtypedef    long FILEID;
X
X
Xstruct value {
X    short v_type;            /* type of value */
X    short v_subtype;        /* other data related to some types */
X    union {
X        long vv_int;        /* small integer value */
X        FILEID vv_file;        /* id of opened file */
X        NUMBER *vv_num;        /* arbitrary sized numeric value */
X        COMPLEX *vv_com;    /* complex number */
X        VALUE *vv_addr;        /* address of variable value */
X        MATRIX *vv_mat;        /* address of matrix */
X        LIST *vv_list;        /* address of list */
X        ASSOC *vv_assoc;    /* address of association */
X        OBJECT *vv_obj;        /* address of object */
X        char *vv_str;        /* string value */
X    } v_union;
X};
X
X
X/*
X * For ease in referencing
X */
X#define v_int    v_union.vv_int
X#define    v_file    v_union.vv_file
X#define v_num    v_union.vv_num
X#define v_com    v_union.vv_com
X#define v_addr    v_union.vv_addr
X#define v_str    v_union.vv_str
X#define v_mat    v_union.vv_mat
X#define    v_list    v_union.vv_list
X#define    v_assoc    v_union.vv_assoc
X#define v_obj    v_union.vv_obj
X#define    v_valid    v_union.vv_int
X
X
X/*
X * Value types.
X */
X#define V_NULL    0    /* null value */
X#define V_INT    1    /* normal integer */
X#define V_NUM    2    /* number */
X#define V_COM    3    /* complex number */
X#define V_ADDR    4    /* address of variable value */
X#define V_STR    5    /* address of string */
X#define V_MAT    6    /* address of matrix structure */
X#define    V_LIST    7    /* address of list structure */
X#define    V_ASSOC    8    /* address of association structure */
X#define V_OBJ    9    /* address of object structure */
X#define    V_FILE    10    /* opened file id */
X#define V_MAX    10    /* highest legal value */
X
X#define V_STRLITERAL    0    /* string subtype for literal str */
X#define V_STRALLOC    1    /* string subtype for allocated str */
X
X#define TWOVAL(a,b) ((a) * (V_MAX+1) + (b))    /* for switch of two values */
X
X#define    NULL_VALUE    ((VALUE *) 0)
X
X
Xextern void freevalue MATH_PROTO((VALUE *vp));
Xextern void copyvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void negvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void addvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void subvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void mulvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void squarevalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void invertvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void roundvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void broundvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void intvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void fracvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void incvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void decvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void conjvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void sqrtvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void rootvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *v3,
X    VALUE *vres));
Xextern void absvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void normvalue MATH_PROTO((VALUE *vp, VALUE *vres));
Xextern void shiftvalue MATH_PROTO((VALUE *v1, VALUE *v2, BOOL rightshift,
X    VALUE *vres));
Xextern void scalevalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void powivalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void powervalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *v3,
X    VALUE *vres));
Xextern void divvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void quovalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern void modvalue MATH_PROTO((VALUE *v1, VALUE *v2, VALUE *vres));
Xextern BOOL testvalue MATH_PROTO((VALUE *vp));
Xextern BOOL comparevalue MATH_PROTO((VALUE *v1, VALUE *v2));
Xextern FLAG relvalue MATH_PROTO((VALUE *v1, VALUE *v2));
Xextern HASH hashvalue MATH_PROTO((VALUE *vp));
Xextern void printvalue MATH_PROTO((VALUE *vp, int flags));
X
X
X
X/*
X * Structure of a matrix.
X */
Xstruct matrix {
X    long m_dim;        /* dimension of matrix */
X    long m_size;        /* total number of elements */
X    long m_min[MAXDIM];    /* minimum bound for indices */
X    long m_max[MAXDIM];    /* maximum bound for indices */
X    VALUE m_table[1];    /* actually varying length table */
X};
X
X#define matsize(n) (sizeof(MATRIX) - sizeof(VALUE) + ((n) * sizeof(VALUE)))
X
X
Xextern MATRIX *matadd MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern MATRIX *matsub MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern MATRIX *matmul MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern MATRIX *matneg MATH_PROTO((MATRIX *m));
Xextern MATRIX *matalloc MATH_PROTO((long size));
Xextern MATRIX *matcopy MATH_PROTO((MATRIX *m));
Xextern MATRIX *matsquare MATH_PROTO((MATRIX *m));
Xextern MATRIX *matinv MATH_PROTO((MATRIX *m));
Xextern MATRIX *matscale MATH_PROTO((MATRIX *m, long n));
Xextern MATRIX *matshift MATH_PROTO((MATRIX *m, long n));
Xextern MATRIX *matmulval MATH_PROTO((MATRIX *m, VALUE *vp));
Xextern MATRIX *matpowi MATH_PROTO((MATRIX *m, NUMBER *q));
Xextern MATRIX *matconj MATH_PROTO((MATRIX *m));
Xextern MATRIX *matquoval MATH_PROTO((MATRIX *m, VALUE *vp));
Xextern MATRIX *matmodval MATH_PROTO((MATRIX *m, VALUE *vp));
Xextern MATRIX *matint MATH_PROTO((MATRIX *m));
Xextern MATRIX *matfrac MATH_PROTO((MATRIX *m));
Xextern MATRIX *matround MATH_PROTO((MATRIX *m, long places));
Xextern MATRIX *matbround MATH_PROTO((MATRIX *m, long places));
Xextern MATRIX *mattrans MATH_PROTO((MATRIX *m));
Xextern MATRIX *matcross MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern BOOL mattest MATH_PROTO((MATRIX *m));
Xextern BOOL matcmp MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern long matsearch MATH_PROTO((MATRIX *m, VALUE *vp, long index));
Xextern long matrsearch MATH_PROTO((MATRIX *m, VALUE *vp, long index));
Xextern HASH mathash MATH_PROTO((MATRIX *m));
Xextern VALUE matdet MATH_PROTO((MATRIX *m));
Xextern VALUE matdot MATH_PROTO((MATRIX *m1, MATRIX *m2));
Xextern void matfill MATH_PROTO((MATRIX *m, VALUE *v1, VALUE *v2));
Xextern void matfree MATH_PROTO((MATRIX *m));
Xextern void matprint MATH_PROTO((MATRIX *m, long max_print));
Xextern VALUE *matindex MATH_PROTO((MATRIX *mp, BOOL create, long dim,
X    VALUE *indices));
X
X
X#if 0
Xextern BOOL matisident MATH_PROTO((MATRIX *m));
X#endif
X
X
X
X/*
X * List definitions.
X * An individual list element.
X */
Xtypedef struct listelem LISTELEM;
Xstruct listelem {
X    LISTELEM *e_next;    /* next element in list (or NULL) */
X    LISTELEM *e_prev;    /* previous element in list (or NULL) */
X    VALUE e_value;        /* value of this element */
X};
X
X
X/*
X * Structure for a list of elements.
X */
Xstruct list {
X    LISTELEM *l_first;    /* first list element (or NULL) */
X    LISTELEM *l_last;    /* last list element (or NULL) */
X    LISTELEM *l_cache;    /* cached list element (or NULL) */
X    long l_cacheindex;    /* index of cached element (or undefined) */
X    long l_count;        /* total number of elements in the list */
X};
X
X
Xextern void insertlistfirst MATH_PROTO((LIST *lp, VALUE *vp));
Xextern void insertlistlast MATH_PROTO((LIST *lp, VALUE *vp));
Xextern void insertlistmiddle MATH_PROTO((LIST *lp, long index, VALUE *vp));
Xextern void removelistfirst MATH_PROTO((LIST *lp, VALUE *vp));
Xextern void removelistlast MATH_PROTO((LIST *lp, VALUE *vp));
Xextern void removelistmiddle MATH_PROTO((LIST *lp, long index, VALUE *vp));
Xextern void listfree MATH_PROTO((LIST *lp));
Xextern void listprint MATH_PROTO((LIST *lp, long max_print));
Xextern long listsearch MATH_PROTO((LIST *lp, VALUE *vp, long index));
Xextern long listrsearch MATH_PROTO((LIST *lp, VALUE *vp, long index));
Xextern HASH listhash MATH_PROTO((LIST *lp));
Xextern BOOL listcmp MATH_PROTO((LIST *lp1, LIST *lp2));
Xextern VALUE *listfindex MATH_PROTO((LIST *lp, long index));
Xextern LIST *listalloc MATH_PROTO((void));
Xextern LIST *listcopy MATH_PROTO((LIST *lp));
X
X
X
X/*
X * Structures for associations.
X * Associations are "indexed" by one or more arbitrary values, and are
X * stored in a hash table with their hash values for quick indexing.
X */
Xtypedef    struct assocelem ASSOCELEM;
Xstruct assocelem {
X    ASSOCELEM *e_next;    /* next element in list (or NULL) */
X    long e_dim;        /* dimension of indexing for this element */
X    HASH e_hash;        /* hash value for this element */
X    VALUE e_value;        /* value of association */
X    VALUE e_indices[1];    /* index values (variable length) */
X};
X
X
Xstruct assoc {
X    long a_count;        /* number of elements in the association */
X    long a_size;        /* current size of association hash table */
X    ASSOCELEM **a_table;    /* current hash table for elements */
X};
X
X
Xextern ASSOC *assocalloc MATH_PROTO((long initsize));
Xextern ASSOC *assoccopy MATH_PROTO((ASSOC *ap));
Xextern void assocfree MATH_PROTO((ASSOC *ap));
Xextern void assocprint MATH_PROTO((ASSOC *ap, long max_print));
Xextern long assocsearch MATH_PROTO((ASSOC *ap, VALUE *vp, long index));
Xextern long assocrsearch MATH_PROTO((ASSOC *ap, VALUE *vp, long index));
Xextern HASH assochash MATH_PROTO((ASSOC *ap));
Xextern BOOL assoccmp MATH_PROTO((ASSOC *ap1, ASSOC *ap2));
Xextern VALUE *assocfindex MATH_PROTO((ASSOC *ap, long index));
Xextern VALUE *associndex MATH_PROTO((ASSOC *ap, BOOL create, long dim,
X    VALUE *indices));
X
X
X/*
X * Object actions.
X */
X#define OBJ_PRINT    0    /* print the value */
X#define OBJ_ONE        1    /* create the multiplicative identity */
X#define OBJ_TEST    2    /* test a value for "zero" */
X#define OBJ_ADD        3    /* add two values */
X#define OBJ_SUB        4    /* subtrace one value from another */
X#define OBJ_NEG        5    /* negate a value */
X#define OBJ_MUL        6    /* multiply two values */
X#define OBJ_DIV        7    /* divide one value by another */
X#define OBJ_INV        8    /* invert a value */
X#define OBJ_ABS        9    /* take absolute value of value */
X#define OBJ_NORM    10    /* take the norm of a value */
X#define OBJ_CONJ    11    /* take the conjugate of a value */
X#define OBJ_POW        12    /* take the power function */
X#define OBJ_SGN        13    /* return the sign of a value */
X#define OBJ_CMP        14    /* compare two values for equality */
X#define OBJ_REL        15    /* compare two values for inequality */
X#define OBJ_QUO        16    /* integer quotient of values */
X#define OBJ_MOD        17    /* remainder of division of values */
X#define OBJ_INT        18    /* integer part of */
X#define OBJ_FRAC    19    /* fractional part of */
X#define OBJ_INC        20    /* increment by one */
X#define OBJ_DEC        21    /* decrement by one */
X#define OBJ_SQUARE    22    /* square value */
X#define OBJ_SCALE    23    /* scale by power of two */
X#define OBJ_SHIFT    24    /* shift left (or right) by number of bits */
X#define OBJ_ROUND    25    /* round to specified decimal places */
X#define OBJ_BROUND    26    /* round to specified binary places */
X#define OBJ_ROOT    27    /* take nth root of value */
X#define OBJ_SQRT    28    /* take square root of value */
X#define OBJ_MAXFUNC    28    /* highest function */
X
X
X/*
X * Definition of an object type.
X * This is actually a varying sized structure.
X */
Xtypedef struct {
X    char *name;            /* name of object */
X    int count;            /* number of elements defined */
X    int actions[OBJ_MAXFUNC+1];    /* function indices for actions */
X    int elements[1];        /* element indexes (MUST BE LAST) */
X} OBJECTACTIONS;
X
X#define objectactionsize(elements) \
X    (sizeof(OBJECTACTIONS) + ((elements) - 1) * sizeof(int))
X
X
X/*
X * Structure of an object.
X * This is actually a varying sized structure.
X * However, there are always at least USUAL_ELEMENTS values in the object.
X */
Xstruct object {
X    OBJECTACTIONS *o_actions;    /* action table for this object */
X    VALUE o_table[USUAL_ELEMENTS];    /* object values (MUST BE LAST) */
X};
X
X#define objectsize(elements) \
X    (sizeof(OBJECT) + ((elements) - USUAL_ELEMENTS) * sizeof(VALUE))
X
X
Xextern OBJECT *objcopy MATH_PROTO((OBJECT *op));
Xextern OBJECT *objalloc MATH_PROTO((long index));
Xextern VALUE objcall MATH_PROTO((int action, VALUE *v1, VALUE *v2, VALUE *v3));
Xextern void objfree MATH_PROTO((OBJECT *op));
Xextern void objuncache MATH_PROTO((void));
Xextern int addelement MATH_PROTO((char *name));
Xextern void defineobject MATH_PROTO((char *name, int indices[], int count));
Xextern int checkobject MATH_PROTO((char *name));
Xextern void showobjfuncs MATH_PROTO((void));
Xextern int findelement MATH_PROTO((char *name));
Xextern int objoffset MATH_PROTO((OBJECT *op, long index));
Xextern HASH objhash MATH_PROTO((OBJECT *op));
X
X#endif
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/value.h || echo "restore of calc2.9.0/value.h fails"
set `wc -c calc2.9.0/value.h`;Sum=$1
if test "$Sum" != "12825"
then echo original size 12825, current size $Sum;fi
echo "x - extracting calc2.9.0/version.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/version.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 * version - determine the version of calc
X */
X
X#include "calc.h"
X
X#define MAJOR_VER    2    /* major version */
X#define MINOR_VER    9    /* minor version */
X#define PATCH_LEVEL    0    /* patch level */
X
X
Xvoid
Xversion(stream)
X    FILE *stream;    /* stream to write version on */
X{
X    fprintf(stream,
X        "C-style arbitrary precision calculator (version %d.%d.%d)\n",
X        MAJOR_VER, MINOR_VER, PATCH_LEVEL);
X}
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/version.c || echo "restore of calc2.9.0/version.c fails"
set `wc -c calc2.9.0/version.c`;Sum=$1
if test "$Sum" != "564"
then echo original size 564, current size $Sum;fi
echo "x - extracting calc2.9.0/zfunc.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/zfunc.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Extended precision integral arithmetic non-primitive routines
X */
X
X#include "zmath.h"
X
Xstatic ZVALUE primeprod;        /* product of primes under 100 */
XZVALUE _tenpowers_[32];            /* table of 10^2^n */
X
X
X/*
X * Compute the factorial of a number.
X */
Xvoid
Xzfact(z, dest)
X    ZVALUE z, *dest;
X{
X    long ptwo;        /* count of powers of two */
X    long n;            /* current multiplication value */
X    long m;            /* reduced multiplication value */
X    long mul;        /* collected value to multiply by */
X    ZVALUE res, temp;
X
X    if (zisneg(z))
X        math_error("Negative argument for factorial");
X    if (zisbig(z))
X        math_error("Very large factorial");
X    n = (zistiny(z) ? z1tol(z) : z2tol(z));
X    ptwo = 0;
X    mul = 1;
X    res = _one_;
X    /*
X     * Multiply numbers together, but squeeze out all powers of two.
X     * We will put them back in at the end.  Also collect multiple
X     * numbers together until there is a risk of overflow.
X     */
X    for (; n > 1; n--) {
X        for (m = n; ((m & 0x1) == 0); m >>= 1)
X            ptwo++;
X        mul *= m;
X        if (mul < BASE1/2)
X            continue;
X        zmuli(res, mul, &temp);
X        zfree(res);
X        res = temp;
X        mul = 1;
X    }
X    /*
X     * Multiply by the remaining value, then scale result by
X     * the proper power of two.
X     */
X    if (mul > 1) {
X        zmuli(res, mul, &temp);
X        zfree(res);
X        res = temp;
X    }
X    zshift(res, ptwo, &temp);
X    zfree(res);
X    *dest = temp;
X}
X
X
X/*
X * Compute the product of the primes up to the specified number.
X */
Xvoid
Xzpfact(z, dest)
X    ZVALUE z, *dest;
X{
X    long n;            /* limiting number to multiply by */
X    long p;            /* current prime */
X    long i;            /* test value */
X    long mul;        /* collected value to multiply by */
X    ZVALUE res, temp;
X
X    if (zisneg(z))
X        math_error("Negative argument for factorial");
X    if (zisbig(z))
X        math_error("Very large factorial");
X    n = (zistiny(z) ? z1tol(z) : z2tol(z));
X    /*
X     * Multiply by the primes in order, collecting multiple numbers
X     * together until there is a chance of overflow.
X     */
X    mul = 1 + (n > 1);
X    res = _one_;
X    for (p = 3; p <= n; p += 2) {
X        for (i = 3; (i * i) <= p; i += 2) {
X            if ((p % i) == 0)
X                goto next;
X        }
X        mul *= p;
X        if (mul < BASE1/2)
X            continue;
X        zmuli(res, mul, &temp);
X        zfree(res);
X        res = temp;
X        mul = 1;
Xnext: ;
X    }
X    /*
X     * Multiply by the final value if any.
X     */
X    if (mul > 1) {
X        zmuli(res, mul, &temp);
X        zfree(res);
X        res = temp;
X    }
X    *dest = res;
X}
X
X
X/*
X * Compute the least common multiple of all the numbers up to the
X * specified number.
X */
Xvoid
Xzlcmfact(z, dest)
X    ZVALUE z, *dest;
X{
X    long n;            /* limiting number to multiply by */
X    long p;            /* current prime */
X    long pp;        /* power of prime */
X    long i;            /* test value */
X    ZVALUE res, temp;
X
X    if (zisneg(z) || ziszero(z))
X        math_error("Non-positive argument for lcmfact");
X    if (zisbig(z))
X        math_error("Very large lcmfact");
X    n = (zistiny(z) ? z1tol(z) : z2tol(z));
X    /*
X     * Multiply by powers of the necessary odd primes in order.
X     * The power for each prime is the highest one which is not
X     * more than the specified number.
X     */
X    res = _one_;
X    for (p = 3; p <= n; p += 2) {
X        for (i = 3; (i * i) <= p; i += 2) {
X            if ((p % i) == 0)
X                goto next;
X        }
X        i = p;
X        while (i <= n) {
X            pp = i;
X            i *= p;
X        }
X        zmuli(res, pp, &temp);
X        zfree(res);
X        res = temp;
Xnext: ;
X    }
X    /*
X     * Finish by scaling by the necessary power of two.
X     */
X    zshift(res, zhighbit(z), dest);
X    zfree(res);
X}
X
X
X/*
X * Compute the permutation function  M! / (M - N)!.
X */
Xvoid
Xzperm(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    long count;
X    ZVALUE cur, tmp, ans;
X
X    if (zisneg(z1) || zisneg(z2))
X        math_error("Negative argument for permutation");
X    if (zrel(z1, z2) < 0)
X        math_error("Second arg larger than first in permutation");
X    if (zisbig(z2))
X        math_error("Very large permutation");
X    count = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
X    zcopy(z1, &ans);
X    zsub(z1, _one_, &cur);
X    while (--count > 0) {
X        zmul(ans, cur, &tmp);
X        zfree(ans);
X        ans = tmp;
X        zsub(cur, _one_, &tmp);
X        zfree(cur);
X        cur = tmp;
X    }
X    zfree(cur);
X    *res = ans;
X}
X
X
X/*
X * Compute the combinatorial function  M! / ( N! * (M - N)! ).
X */
Xvoid
Xzcomb(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE ans;
X    ZVALUE mul, div, temp;
X    FULL count, i;
X    HALF dh[2];
X
X    if (zisneg(z1) || zisneg(z2))
X        math_error("Negative argument for combinatorial");
X    zsub(z1, z2, &temp);
X    if (zisneg(temp)) {
X        zfree(temp);
X        math_error("Second arg larger than first for combinatorial");
X    }
X    if (zisbig(z2) && zisbig(temp)) {
X        zfree(temp);
X        math_error("Very large combinatorial");
X    }
X    count = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
X    i = (zistiny(temp) ? z1tol(temp) : z2tol(temp));
X    if (zisbig(z2) || (!zisbig(temp) && (i < count)))
X        count = i;
X    zfree(temp);
X    mul = z1;
X    div.sign = 0;
X    div.v = dh;
X    ans = _one_;
X    for (i = 1; i <= count; i++) {
X        dh[0] = i & BASE1;
X        dh[1] = i / BASE;
X        div.len = 1 + (dh[1] != 0);
X        zmul(ans, mul, &temp);
X        zfree(ans);
X        zquo(temp, div, &ans);
X        zfree(temp);
X        zsub(mul, _one_, &temp);
X        if (mul.v != z1.v)
X            zfree(mul);
X        mul = temp;
X    }
X    if (mul.v != z1.v)
X        zfree(mul);
X    *res = ans;
X}
X
X
X/*
X * Perform a probabilistic primality test (algorithm P in Knuth).
X * Returns FALSE if definitely not prime, or TRUE if probably prime.
X * Count determines how many times to check for primality.
X * The chance of a non-prime passing this test is less than (1/4)^count.
X * For example, a count of 100 fails for only 1 in 10^60 numbers.
X */
XBOOL
Xzprimetest(z, count)
X    ZVALUE z;        /* number to test for primeness */
X    long count;
X{
X    long ij, ik, ix;
X    ZVALUE zm1, z1, z2, z3, ztmp;
X    HALF val[2];
X
X    z.sign = 0;
X    if (ziseven(z))        /* if even, not prime if not 2 */
X        return (zistwo(z) != 0);
X    /*
X     * See if the number is small, and is either a small prime,
X     * or is divisible by a small prime.
X     */
X    if (zistiny(z) && (*z.v <= (HALF)(101*101-1))) {
X        ix = *z.v;
X        for (ik = 3; (ik <= 97) && ((ik * ik) <= ix); ik += 2)
X            if ((ix % ik) == 0)
X                return FALSE;
X        return TRUE;
X    }
X    /*
X     * See if the number is divisible by one of the primes 3, 5,
X     * 7, 11, or 13.  This is a very easy check.
X     */
X    ij = zmodi(z, 15015L);
X    if (!(ij % 3) || !(ij % 5) || !(ij % 7) || !(ij % 11) || !(ij % 13))
X        return FALSE;
X    /*
X     * Check the gcd of the number and the product of more of the first
X     * few odd primes.  We must build the prime product on the first call.
X     */
X    ztmp.sign = 0;
X    ztmp.len = 1;
X    ztmp.v = val;
X    if (primeprod.len == 0) {
X        val[0] = 101;
X        zpfact(ztmp, &primeprod);
X    }
X    zgcd(z, primeprod, &z1);
X    if (!zisunit(z1)) {
X        zfree(z1);
X        return FALSE;
X    }
X    zfree(z1);
X    /*
X     * Not divisible by a small prime, so onward with the real test.
X     * Make sure the count is limited by the number of odd numbers between
X     * three and the number being tested.
X     */
X    ix = ((zistiny(z) ? z1tol(z) : z2tol(z) - 3) / 2);
X    if (count > ix) count = ix;
X    zsub(z, _one_, &zm1);
X    ik = zlowbit(zm1);
X    zshift(zm1, -ik, &z1);
X    /*
X     * Loop over various "random" numbers, testing each one.
X     * These numbers are the odd numbers starting from three.
X     */
X    for (ix = 0; ix < count; ix++) {
X        val[0] = (ix * 2) + 3;
X        ij = 0;
X        zpowermod(ztmp, z1, z, &z3);
X        for (;;) {
X            if (zisone(z3)) {
X                if (ij)    /* number is definitely not prime */
X                    goto notprime;
X                break;
X            }
X            if (zcmp(z3, zm1) == 0)
X                break;
X            if (++ij >= ik)
X                goto notprime;    /* number is definitely not prime */
X            zsquare(z3, &z2);
X            zfree(z3);
X            zmod(z2, z, &z3);
X            zfree(z2);
X        }
X        zfree(z3);
X    }
X    zfree(zm1);
X    zfree(z1);
X    return TRUE;    /* number might be prime */
X
Xnotprime:
X    zfree(z3);
X    zfree(zm1);
X    zfree(z1);
X    return FALSE;
X}
X
X
X/*
X * Compute the Jacobi function (p / q) for odd q.
X * If q is prime then the result is:
X *    1 if p == x^2 (mod q) for some x.
X *    -1 otherwise.
X * If q is not prime, then the result is not meaningful if it is 1.
X * This function returns 0 if q is even or q < 0.
X */
XFLAG
Xzjacobi(z1, z2)
X    ZVALUE z1, z2;
X{
X    ZVALUE p, q, tmp;
X    long lowbit;
X    int val;
X
X    if (ziseven(z2) || zisneg(z2))
X        return 0;
X    val = 1;
X    if (ziszero(z1) || zisone(z1))
X        return val;
X    if (zisunit(z1)) {
X        if ((*z2.v - 1) & 0x2)
X            val = -val;
X        return val;
X    }
X    zcopy(z1, &p);
X    zcopy(z2, &q);
X    for (;;) {
X        zmod(p, q, &tmp);
X        zfree(p);
X        p = tmp;
X        if (ziszero(p)) {
X            zfree(p);
X            p = _one_;
X        }
X        if (ziseven(p)) {
X            lowbit = zlowbit(p);
X            zshift(p, -lowbit, &tmp);
X            zfree(p);
X            p = tmp;
X            if ((lowbit & 1) && (((*q.v & 0x7) == 3) || ((*q.v & 0x7) == 5)))
X                val = -val;
X        }
X        if (zisunit(p)) {
X            zfree(p);
X            zfree(q);
X            return val;
X        }
X        if ((*p.v & *q.v & 0x3) == 3)
X            val = -val;
X        tmp = q;
X        q = p;
X        p = tmp;
X    }
X}
X
X
X/*
X * Return the Fibonacci number F(n).
X * This is evaluated by recursively using the formulas:
X *    F(2N+1) = F(N+1)^2 + F(N)^2
X * and
X *    F(2N) = F(N+1)^2 - F(N-1)^2
X */
Xvoid
Xzfib(z, res)
X    ZVALUE z, *res;
X{
X    unsigned long i;
X    long n;
X    int sign;
X    ZVALUE fnm1, fn, fnp1;        /* consecutive fibonacci values */
X    ZVALUE t1, t2, t3;
X
X    if (zisbig(z))
X        math_error("Very large Fibonacci number");
X    n = (zistiny(z) ? z1tol(z) : z2tol(z));
X    if (n == 0) {
X        *res = _zero_;
X        return;
X    }
X    sign = z.sign && ((n & 0x1) == 0);
X    if (n <= 2) {
X        *res = _one_;
X        res->sign = (BOOL)sign;
X        return;
X    }
X    i = TOPFULL;
X    while ((i & n) == 0)
X        i >>= 1L;
X    i >>= 1L;
X    fnm1 = _zero_;
X    fn = _one_;
X    fnp1 = _one_;
X    while (i) {
X        zsquare(fnm1, &t1);
X        zsquare(fn, &t2);
X        zsquare(fnp1, &t3);
X        zfree(fnm1);
X        zfree(fn);
X        zfree(fnp1);
X        zadd(t2, t3, &fnp1);
X        zsub(t3, t1, &fn);
X        zfree(t1);
X        zfree(t2);
X        zfree(t3);
X        if (i & n) {
X            fnm1 = fn;
X            fn = fnp1;
X            zadd(fnm1, fn, &fnp1);
X        } else
X            zsub(fnp1, fn, &fnm1);
X        i >>= 1L;
X    }
X    zfree(fnm1);
X    zfree(fnp1);
X    *res = fn;
X    res->sign = (BOOL)sign;
X}
X
X
X/*
X * Compute the result of raising one number to the power of another
X * The second number is assumed to be non-negative.
X * It cannot be too large except for trivial cases.
X */
Xvoid
Xzpowi(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    int sign;        /* final sign of number */
X    unsigned long power;    /* power to raise to */
X    unsigned long bit;    /* current bit value */
X    long twos;        /* count of times 2 is in result */
X    ZVALUE ans, temp;
X
X    sign = (z1.sign && zisodd(z2));
X    z1.sign = 0;
X    z2.sign = 0;
X    if (ziszero(z2)) {    /* number raised to power 0 */
X        if (ziszero(z1))
X            math_error("Zero raised to zero power");
X        *res = _one_;
X        return;
X    }
X    if (zisleone(z1)) {    /* 0, 1, or -1 raised to a power */
X        ans = _one_;
X        ans.sign = (BOOL)sign;
X        if (*z1.v == 0)
X            ans = _zero_;
X        *res = ans;
X        return;
X    }
X    if (zisbig(z2))
X        math_error("Raising to very large power");
X    power = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
X    if (zistwo(z1)) {    /* two raised to a power */
X        zbitvalue((long) power, res);
X        return;
X    }
X    /*
X     * See if this is a power of ten
X     */
X    if (zistiny(z1) && (*z1.v == 10)) {
X        ztenpow((long) power, res);
X        res->sign = (BOOL)sign;
X        return;
X    }
X    /*
X     * Handle low powers specially
X     */
X    if (power <= 4) {
X        switch ((int) power) {
X            case 1:
X                ans.len = z1.len;
X                ans.v = alloc(ans.len);
X                zcopyval(z1, ans);
X                ans.sign = (BOOL)sign;
X                *res = ans;
X                return;
X            case 2:
X                zsquare(z1, res);
X                return;
X            case 3:
X                zsquare(z1, &temp);
X                zmul(z1, temp, res);
X                zfree(temp);
X                res->sign = (BOOL)sign;
X                return;
X            case 4:
X                zsquare(z1, &temp);
X                zsquare(temp, res);
X                zfree(temp);
X                return;
X        }
X    }
X    /*
X     * Shift out all powers of twos so the multiplies are smaller.
X     * We will shift back the right amount when done.
X     */
X    twos = 0;
X    if (ziseven(z1)) {
X        twos = zlowbit(z1);
X        ans.v = alloc(z1.len);
X        ans.len = z1.len;
X        zcopyval(z1, ans);
X        zshiftr(ans, twos);
X        ztrim(&ans);
X        z1 = ans;
X        twos *= power;
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    zsquare(z1, &ans);
X    if (bit & power) {
X        zmul(ans, z1, &temp);
X        zfree(ans);
X        ans = temp;
X    }
X    bit >>= 1L;
X    while (bit) {
X        zsquare(ans, &temp);
X        zfree(ans);
X        ans = temp;
X        if (bit & power) {
X            zmul(ans, z1, &temp);
X            zfree(ans);
X            ans = temp;
X        }
X        bit >>= 1L;
X    }
X    /*
X     * Scale back up by proper power of two
X     */
X    if (twos) {
X        zshift(ans, twos, &temp);
X        zfree(ans);
X        ans = temp;
X        zfree(z1);
X    }
X    ans.sign = (BOOL)sign;
X    *res = ans;
X}
X
X
X/*
X * Compute ten to the specified power
X * This saves some work since the squares of ten are saved.
X */
Xvoid
Xztenpow(power, res)
X    long power;
X    ZVALUE *res;
X{
X    long i;
X    ZVALUE ans;
X    ZVALUE temp;
X
X    if (power <= 0) {
X        *res = _one_;
X        return;
X    }
X    ans = _one_;
X    _tenpowers_[0] = _ten_;
X    for (i = 0; power; i++) {
X        if (_tenpowers_[i].len == 0)
X            zsquare(_tenpowers_[i-1], &_tenpowers_[i]);
X        if (power & 0x1) {
X            zmul(ans, _tenpowers_[i], &temp);
X            zfree(ans);
X            ans = temp;
X        }
X        power /= 2;
X    }
X    *res = ans;
X}
X
X
X/*
X * Calculate modular inverse suppressing unnecessary divisions.
X * This is based on the Euclidian algorithm for large numbers.
X * (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
X * Returns TRUE if there is no solution because the numbers
X * are not relatively prime.
X */
XBOOL
Xzmodinv(u, v, res)
X    ZVALUE u, v;
X    ZVALUE *res;
X{
X    FULL    q1, q2, ui3, vi3, uh, vh, A, B, C, D, T;
X    ZVALUE    u2, u3, v2, v3, qz, tmp1, tmp2, tmp3;
X
X    if (zisneg(u) || zisneg(v) || (zrel(u, v) >= 0))
X        zmod(u, v, &v3);
X    else
X        zcopy(u, &v3);
X    zcopy(v, &u3);
X    u2 = _zero_;
X    v2 = _one_;
X
X    /*
X     * Loop here while the size of the numbers remain above
X     * the size of a FULL.  Throughout this loop u3 >= v3.
X     */
X    while ((u3.len > 1) && !ziszero(v3)) {
X        uh = (((FULL) u3.v[u3.len - 1]) << BASEB) + u3.v[u3.len - 2];
X        vh = 0;
X        if ((v3.len + 1) >= u3.len)
X            vh = v3.v[v3.len - 1];
X        if (v3.len == u3.len)
X            vh = (vh << BASEB) + v3.v[v3.len - 2];
X        A = 1;
X        B = 0;
X        C = 0;
X        D = 1;
X
X        /*
X         * Calculate successive quotients of the continued fraction
X         * expansion using only single precision arithmetic until
X         * greater precision is required.
X         */
X        while ((vh + C) && (vh + D)) {
X            q1 = (uh + A) / (vh + C);
X            q2 = (uh + B) / (vh + D);
X            if (q1 != q2)
X                break;
X            T = A - q1 * C;
X            A = C;
X            C = T;
X            T = B - q1 * D;
X            B = D;
X            D = T;
X            T = uh - q1 * vh;
X            uh = vh;
X            vh = T;
X        }
X
X        /*
X         * If B is zero, then we made no progress because
X         * the calculation requires a very large quotient.
X         * So we must do this step of the calculation in
X         * full precision
X         */
X        if (B == 0) {
X            zquo(u3, v3, &qz);
X            zmul(qz, v2, &tmp1);
X            zsub(u2, tmp1, &tmp2);
X            zfree(tmp1);
X            zfree(u2);
X            u2 = v2;
X            v2 = tmp2;
X            zmul(qz, v3, &tmp1);
X            zsub(u3, tmp1, &tmp2);
X            zfree(tmp1);
X            zfree(u3);
X            u3 = v3;
X            v3 = tmp2;
X            zfree(qz);
X            continue;
X        }
X        /*
X         * Apply the calculated A,B,C,D numbers to the current
X         * values to update them as if the full precision
X         * calculations had been carried out.
X         */
X        zmuli(u2, (long) A, &tmp1);
X        zmuli(v2, (long) B, &tmp2);
X        zadd(tmp1, tmp2, &tmp3);
X        zfree(tmp1);
X        zfree(tmp2);
X        zmuli(u2, (long) C, &tmp1);
X        zmuli(v2, (long) D, &tmp2);
X        zfree(u2);
X        zfree(v2);
X        u2 = tmp3;
X        zadd(tmp1, tmp2, &v2);
X        zfree(tmp1);
X        zfree(tmp2);
X        zmuli(u3, (long) A, &tmp1);
X        zmuli(v3, (long) B, &tmp2);
X        zadd(tmp1, tmp2, &tmp3);
X        zfree(tmp1);
X        zfree(tmp2);
X        zmuli(u3, (long) C, &tmp1);
X        zmuli(v3, (long) D, &tmp2);
X        zfree(u3);
X        zfree(v3);
X        u3 = tmp3;
X        zadd(tmp1, tmp2, &v3);
X        zfree(tmp1);
X        zfree(tmp2);
X    }
X
X    /*
X     * Here when the remaining numbers become single precision in size.
X     * Finish the procedure using single precision calculations.
X     */
X    if (ziszero(v3) && !zisone(u3)) {
X        zfree(u3);
X        zfree(v3);
X        zfree(u2);
X        zfree(v2);
X        return TRUE;
X    }
X    ui3 = (zistiny(u3) ? z1tol(u3) : z2tol(u3));
X    vi3 = (zistiny(v3) ? z1tol(v3) : z2tol(v3));
X    zfree(u3);
X    zfree(v3);
X    while (vi3) {
X        q1 = ui3 / vi3;
X        zmuli(v2, (long) q1, &tmp1);
X        zsub(u2, tmp1, &tmp2);
X        zfree(tmp1);
X        zfree(u2);
X        u2 = v2;
X        v2 = tmp2;
X        q2 = ui3 - q1 * vi3;
X        ui3 = vi3;
X        vi3 = q2;
X    }
X    zfree(v2);
X    if (ui3 != 1) {
X        zfree(u2);
X        return TRUE;
X    }
X    if (zisneg(u2)) {
X        zadd(v, u2, res);
X        zfree(u2);
X        return FALSE;
X    }
X    *res = u2;
X    return FALSE;
X}
X
X
X#if 0
X/*
X * Approximate the quotient of two integers by another set of smaller
X * integers.  This uses continued fractions to determine the smaller set.
X */
Xvoid
Xzapprox(z1, z2, res1, res2)
X    ZVALUE z1, z2, *res1, *res2;
X{
X    int sign;
X    ZVALUE u1, v1, u3, v3, q, t1, t2, t3;
X
X    sign = ((z1.sign != 0) ^ (z2.sign != 0));
X    z1.sign = 0;
X    z2.sign = 0;
X    v3 = z2;
X    u3 = z1;
X    u1 = _one_;
X    v1 = _zero_;
X    while (!ziszero(v3)) {
X        zdiv(u3, v3, &q, &t1);
X        zmul(v1, q, &t2);
X        zsub(u1, t2, &t3);
X        zfree(q);
X        zfree(t2);
X        zfree(u1);
X        if ((u3.v != z1.v) && (u3.v != z2.v))
X            zfree(u3);
X        u1 = v1;
X        u3 = v3;
X        v1 = t3;
X        v3 = t1;
X    }
X    if (!zisunit(u3))
X        math_error("Non-relativly prime numbers for approx");
X    if ((u3.v != z1.v) && (u3.v != z2.v))
X        zfree(u3);
X    if ((v3.v != z1.v) && (v3.v != z2.v))
X        zfree(v3);
X    zfree(v1);
X    zmul(u1, z1, &t1);
X    zsub(t1, _one_, &t2);
X    zfree(t1);
X    zquo(t2, z2, &t1);
X    zfree(t2);
X    u1.sign = (BOOL)sign;
X    t1.sign = 0;
X    *res1 = t1;
X    *res2 = u1;
X}
X#endif
X
X
X/*
X * Binary gcd algorithm
X * This algorithm taken from Knuth
X */
Xvoid
Xzgcd(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE u, v, t;
X    register long j, k, olen, mask;
X    register HALF h;
X    HALF *oldv1, *oldv2;
X
X    z1.sign = 0;
X    z2.sign = 0;
X    if (ziszero(z1)) {
X        zcopy(z2, res);
X        return;
X    }
X    if (ziszero(z2)) {
X        zcopy(z1, res);
X        return;
X    }
X    if (zisunit(z1) || zisunit(z2)) {
X        *res = _one_;
X        return;
X    }
X    /*
X     * First see if one number is very much larger than the other.
X     * If so, then divide as necessary to get the numbers near each other.
X     */
X    oldv1 = z1.v;
X    oldv2 = z2.v;
X    if (z1.len < z2.len) {
X        t = z1;
X        z1 = z2;
X        z2 = t;
X    }
X    while ((z1.len > (z2.len + 5)) && !ziszero(z2)) {
X        zmod(z1, z2, &t);
X        if ((z1.v != oldv1) && (z1.v != oldv2))
X            zfree(z1);
X        z1 = z2;
X        z2 = t;
X    }
X    /*
X     * Ok, now do the binary method proper
X     */
X    u.len = z1.len;
X    v.len = z2.len;
X    u.sign = 0;
X    v.sign = 0;
X    if (!ztest(z1)) {
X        v.v = alloc(v.len);
X        zcopyval(z2, v);
X        *res = v;
X        goto done;
X    }
X    if (!ztest(z2)) {
X        u.v = alloc(u.len);
X        zcopyval(z1, u);
X        *res = u;
X        goto done;
X    }
X    u.v = alloc(u.len);
X    v.v = alloc(v.len);
X    zcopyval(z1, u);
X    zcopyval(z2, v);
X    k = 0;
X    while (u.v[k] == 0 && v.v[k] == 0)
X        ++k;
X    mask = 01;
X    h = u.v[k] | v.v[k];
X    k *= BASEB;
X    while (!(h & mask)) {
X        mask <<= 1;
X        k++;
X    }
X    zshiftr(u, k);
X    zshiftr(v, k);
X    ztrim(&u);
X    ztrim(&v);
X    if (zisodd(u)) {
X        t.v = alloc(v.len);
X        t.len = v.len;
X        zcopyval(v, t);
X        t.sign = !v.sign;
X    } else {
X        t.v = alloc(u.len);
X        t.len = u.len;
X        zcopyval(u, t);
X        t.sign = u.sign;
X    }
X    while (ztest(t)) {
X        j = 0;
X        while (!t.v[j])
X            ++j;
X        mask = 01;
X        h = t.v[j];
X        j *= BASEB;
X        while (!(h & mask)) {
X            mask <<= 1;
X            j++;
X        }
X        zshiftr(t, j);
X        ztrim(&t);
X        if (ztest(t) > 0) {
X            zfree(u);
X            u = t;
X        } else {
X            zfree(v);
X            v = t;
X            v.sign = !t.sign;
X        }
X        zsub(u, v, &t);
X    }
X    zfree(t);
X    zfree(v);
X    if (k) {
X        olen = u.len;
X        u.len += k / BASEB + 1;
X        u.v = (HALF *) realloc(u.v, u.len * sizeof(HALF));
X        while (olen != u.len)
X            u.v[olen++] = 0;
X        zshiftl(u, k);
X    }
X    ztrim(&u);
X    *res = u;
X
Xdone:
X    if ((z1.v != oldv1) && (z1.v != oldv2))
X        zfree(z1);
X    if ((z2.v != oldv1) && (z2.v != oldv2))
X        zfree(z2);
X}
X
X
X/*
X * Compute the lcm of two integers (least common multiple).
X * This is done using the formula:  gcd(a,b) * lcm(a,b) = a * b.
X */
Xvoid
Xzlcm(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE temp1, temp2;
X
X    zgcd(z1, z2, &temp1);
X    zquo(z1, temp1, &temp2);
X    zfree(temp1);
X    zmul(temp2, z2, res);
X    zfree(temp2);
X}
X
X
X/*
X * Return whether or not two numbers are relatively prime to each other.
X */
XBOOL
Xzrelprime(z1, z2)
X    ZVALUE z1, z2;            /* numbers to be tested */
X{
X    FULL rem1, rem2;        /* remainders */
X    ZVALUE rem;
X    BOOL result;
X
X    z1.sign = 0;
X    z2.sign = 0;
X    if (ziseven(z1) && ziseven(z2))    /* false if both even */
X        return FALSE;
X    if (zisunit(z1) || zisunit(z2))    /* true if either is a unit */
X        return TRUE;
X    if (ziszero(z1) || ziszero(z2))    /* false if either is zero */
X        return FALSE;
X    if (zistwo(z1) || zistwo(z2))    /* true if either is two */
X        return TRUE;
X    /*
X     * Try reducing each number by the product of the first few odd primes
X     * to see if any of them are a common factor.
X     */
X    rem1 = zmodi(z1, 3L * 5 * 7 * 11 * 13);
X    rem2 = zmodi(z2, 3L * 5 * 7 * 11 * 13);
X    if (((rem1 % 3) == 0) && ((rem2 % 3) == 0))
X        return FALSE;
X    if (((rem1 % 5) == 0) && ((rem2 % 5) == 0))
X        return FALSE;
X    if (((rem1 % 7) == 0) && ((rem2 % 7) == 0))
X        return FALSE;
X    if (((rem1 % 11) == 0) && ((rem2 % 11) == 0))
X        return FALSE;
X    if (((rem1 % 13) == 0) && ((rem2 % 13) == 0))
X        return FALSE;
X    /*
X     * Try a new batch of primes now
X     */
X    rem1 = zmodi(z1, 17L * 19 * 23);
X    rem2 = zmodi(z2, 17L * 19 * 23);
X    if (((rem1 % 17) == 0) && ((rem2 % 17) == 0))
X        return FALSE;
X    if (((rem1 % 19) == 0) && ((rem2 % 19) == 0))
X        return FALSE;
X    if (((rem1 % 23) == 0) && ((rem2 % 23) == 0))
X        return FALSE;
X    /*
X     * Yuk, we must actually compute the gcd to know the answer
X     */
X    zgcd(z1, z2, &rem);
X    result = zisunit(rem);
X    zfree(rem);
X    return result;
X}
X
X
X/*
X * Compute the log of one number base another, to the closest integer.
X * This is the largest integer which when the second number is raised to it,
X * the resulting value is less than or equal to the first number.
X * Example:  zlog(123456, 10) = 5.
X */
Xlong
Xzlog(z1, z2)
X    ZVALUE z1, z2;
X{
X    register ZVALUE *zp;        /* current square */
X    long power;            /* current power */
X    long worth;            /* worth of current square */
X    ZVALUE val;            /* current value of power */
X    ZVALUE temp;            /* temporary */
X    ZVALUE squares[32];        /* table of squares of base */
X
X    /*
X     * Make sure that the numbers are nonzero and the base is greater than one.
X     */
X    if (zisneg(z1) || ziszero(z1) || zisneg(z2) || zisleone(z2))
X        math_error("Bad arguments for log");
X    /*
X     * Reject trivial cases.
X     */
X    if (z1.len < z2.len)
X        return 0;
X    if ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))
X        return 0;
X    power = zrel(z1, z2);
X    if (power <= 0)
X        return (power + 1);
X    /*
X     * Handle any power of two special.
X     */
X    if (zisonebit(z2))
X        return (zhighbit(z1) / zlowbit(z2));
X    /*
X     * Handle base 10 special
X     */
X    if ((z2.len == 1) && (*z2.v == 10))
X        return zlog10(z1);
X    /*
X     * Now loop by squaring the base each time, and see whether or
X     * not each successive square is still smaller than the number.
X     */
X    worth = 1;
X    zp = &squares[0];
X    *zp = z2;
X    while (((zp->len * 2) - 1) <= z1.len) {    /* while square not too large */
X        zsquare(*zp, zp + 1);
X        zp++;
X        worth *= 2;
X    }
X    /*
X     * Now back down the squares, and multiply them together to see
X     * exactly how many times the base can be raised by.
X     */
X    val = _one_;
X    power = 0;
X    for (; zp >= squares; zp--, worth /= 2) {
X        if ((val.len + zp->len - 1) <= z1.len) {
X            zmul(val, *zp, &temp);
X            if (zrel(z1, temp) >= 0) {
X                zfree(val);
X                val = temp;
X                power += worth;
X            } else
X                zfree(temp);
X        }
X        if (zp != squares)
X            zfree(*zp);
X    }
X    return power;
X}
X
X
X/*
X * Return the integral log base 10 of a number.
X */
Xlong
Xzlog10(z)
X    ZVALUE z;
X{
X    register ZVALUE *zp;        /* current square */
X    long power;            /* current power */
X    long worth;            /* worth of current square */
X    ZVALUE val;            /* current value of power */
X    ZVALUE temp;            /* temporary */
X
X    if (!zispos(z))
X        math_error("Non-positive number for log10");
X    /*
X     * Loop by squaring the base each time, and see whether or
X     * not each successive square is still smaller than the number.
X     */
X    worth = 1;
X    zp = &_tenpowers_[0];
X    *zp = _ten_;
X    while (((zp->len * 2) - 1) <= z.len) {    /* while square not too large */
X        if (zp[1].len == 0)
X            zsquare(*zp, zp + 1);
X        zp++;
X        worth *= 2;
X    }
X    /*
X     * Now back down the squares, and multiply them together to see
X     * exactly how many times the base can be raised by.
X     */
X    val = _one_;
X    power = 0;
X    for (; zp >= _tenpowers_; zp--, worth /= 2) {
X        if ((val.len + zp->len - 1) <= z.len) {
X            zmul(val, *zp, &temp);
X            if (zrel(z, temp) >= 0) {
X                zfree(val);
X                val = temp;
X                power += worth;
X            } else
X                zfree(temp);
X        }
X    }
X    return power;
X}
X
X
X/*
X * Return the number of times that one number will divide another.
X * This works similarly to zlog, except that divisions must be exact.
X * For example, zdivcount(540, 3) = 3, since 3^3 divides 540 but 3^4 won't.
X */
Xlong
Xzdivcount(z1, z2)
X    ZVALUE z1, z2;
X{
X    long count;        /* number of factors removed */
X    ZVALUE tmp;        /* ignored return value */
X
X    count = zfacrem(z1, z2, &tmp);
X    zfree(tmp);
X    return count;
X}
X
X
X/*
X * Remove all occurances of the specified factor from a number.
X * Also returns the number of factors removed as a function return value.
X * Example:  zfacrem(540, 3, &x) returns 3 and sets x to 20.
X */
Xlong
Xzfacrem(z1, z2, rem)
X    ZVALUE z1, z2, *rem;
X{
X    register ZVALUE *zp;        /* current square */
X    long count;            /* total count of divisions */
X    long worth;            /* worth of current square */
X    ZVALUE temp1, temp2, temp3;    /* temporaries */
X    ZVALUE squares[32];        /* table of squares of factor */
X
X    z1.sign = 0;
X    z2.sign = 0;
X    /*
X     * Make sure factor isn't 0 or 1.
X     */
X    if (zisleone(z2))
X        math_error("Bad argument for facrem");
X    /*
X     * Reject trivial cases.
X     */
X    if ((z1.len < z2.len) || (zisodd(z1) && ziseven(z2)) ||
X        ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))) {
X        rem->v = alloc(z1.len);
X        rem->len = z1.len;
X        rem->sign = 0;
X        zcopyval(z1, *rem);
X        return 0;
X    }
X    /*
X     * Handle any power of two special.
X     */
X    if (zisonebit(z2)) {
X        count = zlowbit(z1) / zlowbit(z2);
X        rem->v = alloc(z1.len);
X        rem->len = z1.len;
X        rem->sign = 0;
X        zcopyval(z1, *rem);
X        zshiftr(*rem, count);
X        ztrim(rem);
X        return count;
X    }
X    /*
X     * See if the factor goes in even once.
X     */
X    zdiv(z1, z2, &temp1, &temp2);
X    if (!ziszero(temp2)) {
X        zfree(temp1);
X        zfree(temp2);
X        rem->v = alloc(z1.len);
X        rem->len = z1.len;
X        rem->sign = 0;
X        zcopyval(z1, *rem);
X        return 0;
X    }
X    zfree(temp2);
X    z1 = temp1;
X    /*
X     * Now loop by squaring the factor each time, and see whether
X     * or not each successive square will still divide the number.
X     */
X    count = 1;
X    worth = 1;
X    zp = &squares[0];
X    *zp = z2;
X    while (((zp->len * 2) - 1) <= z1.len) {    /* while square not too large */
X        zsquare(*zp, &temp1);
X        zdiv(z1, temp1, &temp2, &temp3);
X        if (!ziszero(temp3)) {
X            zfree(temp1);
X            zfree(temp2);
X            zfree(temp3);
X            break;
X        }
X        zfree(temp3);
X        zfree(z1);
X        z1 = temp2;
X        *++zp = temp1;
X        worth *= 2;
X        count += worth;
X    }
X    /*
X     * Now back down the list of squares, and see if the lower powers
X     * will divide any more times.
X     */
X    for (; zp >= squares; zp--, worth /= 2) {
X        if (zp->len <= z1.len) {
X            zdiv(z1, *zp, &temp1, &temp2);
X            if (ziszero(temp2)) {
X                temp3 = z1;
X                z1 = temp1;
X                temp1 = temp3;
X                count += worth;
X            }
X            zfree(temp1);
X            zfree(temp2);
X        }
X        if (zp != squares)
X            zfree(*zp);
X    }
X    *rem = z1;
X    return count;
X}
X
X
X/*
X * Keep dividing a number by the gcd of it with another number until the
X * result is relatively prime to the second number.
X */
Xvoid
Xzgcdrem(z1, z2, res)
X    ZVALUE z1, z2, *res;
X{
X    ZVALUE tmp1, tmp2;
X
X    /*
X     * Begin by taking the gcd for the first time.
X     * If the number is already relatively prime, then we are done.
X     */
X    zgcd(z1, z2, &tmp1);
X    if (zisunit(tmp1) || ziszero(tmp1)) {
X        res->len = z1.len;
X        res->v = alloc(z1.len);
X        res->sign = z1.sign;
X        zcopyval(z1, *res);
X        return;
X    }
X    zquo(z1, tmp1, &tmp2);
X    z1 = tmp2;
X    z2 = tmp1;
X    /*
X     * Now keep alternately taking the gcd and removing factors until
X     * the gcd becomes one.
X     */
X    while (!zisunit(z2)) {
X        (void) zfacrem(z1, z2, &tmp1);
X        zfree(z1);
X        z1 = tmp1;
X        zgcd(z1, z2, &tmp1);
X        zfree(z2);
X        z2 = tmp1;
X    }
X    *res = z1;
X}
X
X
X/*
X * Find the lowest prime factor of a number if one can be found.
X * Search is conducted for the first count primes.
X * Returns 1 if no factor was found.
X */
Xlong
Xzlowfactor(z, count)
X    ZVALUE z;
X    long count;
X{
X    FULL p, d;
X    ZVALUE div, tmp;
X    HALF divval[2];
X
X    if ((--count < 0) || ziszero(z))
X        return 1;
X    if (ziseven(z))
X        return 2;
X    div.sign = 0;
X    div.v = divval;
X    for (p = 3; (count > 0); p += 2) {
X        for (d = 3; (d * d) <= p; d += 2)
X            if ((p % d) == 0)
X                goto next;
X        divval[0] = (p & BASE1);
X        divval[1] = (p / BASE);
X        div.len = 1 + (p >= BASE);
X        zmod(z, div, &tmp);
X        if (ziszero(tmp))
X            return p;
X        zfree(tmp);
X        count--;
Xnext:;
X    }
X    return 1;
X}
X
X
X/*
X * Return the number of digits (base 10) in a number, ignoring the sign.
X */
Xlong
Xzdigits(z1)
X    ZVALUE z1;
X{
X    long count, val;
X
X    z1.sign = 0;
X    if (zistiny(z1)) {    /* do small numbers ourself */
X        count = 1;
X        val = 10;
X        while (*z1.v >= (HALF)val) {
X            count++;
X            val *= 10;
X        }
X        return count;
X    }
X    return (zlog10(z1) + 1);
X}
X
X
X/*
X * Return the single digit at the specified decimal place of a number,
X * where 0 means the rightmost digit.  Example:  zdigit(1234, 1) = 3.
X */
XFLAG
Xzdigit(z1, n)
X    ZVALUE z1;
X    long n;
X{
X    ZVALUE tmp1, tmp2;
X    FLAG res;
X
X    z1.sign = 0;
X    if (ziszero(z1) || (n < 0) || (n / BASEDIG >= z1.len))
X        return 0;
X    if (n == 0)
X        return zmodi(z1, 10L);
X    if (n == 1)
X        return zmodi(z1, 100L) / 10;
X    if (n == 2)
X        return zmodi(z1, 1000L) / 100;
X    if (n == 3)
X        return zmodi(z1, 10000L) / 1000;
X    ztenpow(n, &tmp1);
X    zquo(z1, tmp1, &tmp2);
X    res = zmodi(tmp2, 10L);
X    zfree(tmp1);
X    zfree(tmp2);
X    return res;
X}
X
X
X/*
X * Find the square root of a number.  This is the greatest integer whose
X * square is less than or equal to the number. Returns TRUE if the
X * square root is exact.
X */
XBOOL
Xzsqrt(z1, dest)
X    ZVALUE z1, *dest;
X{
X    ZVALUE try, quo, rem, old, temp;
X    FULL iquo, val;
X    long i,j;
X
X    if (z1.sign)
X        math_error("Square root of negative number");
X    if (ziszero(z1)) {
X        *dest = _zero_;
X        return TRUE;
X    }
X    if ((*z1.v < 4) && zistiny(z1)) {
X        *dest = _one_;
X        return (*z1.v == 1);
X    }
X    /*
X     * Pick the square root of the leading one or two digits as a first guess.
X     */
X    val = z1.v[z1.len-1];
X    if ((z1.len & 0x1) == 0)
X        val = (val * BASE) + z1.v[z1.len-2];
X
X    /*
X     * Find the largest power of 2 that when squared
X     * is <= val > 0.  Avoid multiply overflow by doing
X     * a careful check at the BASE boundary.
X     */
X    j = 1<<(BASEB+BASEB-2);
X    if (val > j) {
X        iquo = BASE;
X    } else {
X        i = BASEB-1;
X        while (j > val) {
X            --i;
X            j >>= 2;
X        }
X        iquo = bitmask[i];
X    }
X
X    for (i = 8; i > 0; i--)
X        iquo = (iquo + (val / iquo)) / 2;
X    if (iquo > BASE1)
X        iquo = BASE1;
X    /*
X     * Allocate the numbers to use for the main loop.
X     * The size and high bits of the final result are correctly set here.
X     * Notice that the remainder of the test value is rubbish, but this
X     * is unimportant.
X     */
X    try.sign = 0;
X    try.len = (z1.len + 1) / 2;
X    try.v = alloc(try.len);
X    try.v[try.len-1] = (HALF)iquo;
X    old.sign = 0;
X    old.v = alloc(try.len);
X    old.len = 1;
X    *old.v = 0;
X    /*
X     * Main divide and average loop
X     */
X    for (;;) {
X        zdiv(z1, try, &quo, &rem);
X        i = zrel(try, quo);
X        if ((i == 0) && ziszero(rem)) {    /* exact square root */
X            zfree(rem);
X            zfree(quo);
X            zfree(old);
X            *dest = try;
X            return TRUE;
X        }
X        zfree(rem);
X        if (i <= 0) {
X            /*
X            * Current try is less than or equal to the square root since
X            * it is less than the quotient.  If the quotient is equal to
X            * the try, we are all done.  Also, if the try is equal to the
X            * old try value, we are done since no improvement occurred.
X            * If not, save the improved value and loop some more.
X            */
X            if ((i == 0) || (zcmp(old, try) == 0)) {
X                zfree(quo);
X                zfree(old);
X                *dest = try;
X                return FALSE;
X            }
X            old.len = try.len;
X            zcopyval(try, old);
X        }
X        /* average current try and quotent for the new try */
X        zadd(try, quo, &temp);
X        zfree(quo);
X        zfree(try);
X        try = temp;
X        zshiftr(try, 1L);
X        zquicktrim(try);
X    }
X}
X
X
X/*
X * Take an arbitrary root of a number (to the greatest integer).
X * This uses the following iteration to get the Kth root of N:
X *    x = ((K-1) * x + N / x^(K-1)) / K
X */
Xvoid
Xzroot(z1, z2, dest)
X    ZVALUE z1, z2, *dest;
X{
X    ZVALUE try, quo, old, temp, temp2;
X    ZVALUE k1;            /* holds k - 1 */
X    int sign;
X    long i, k, highbit;
X    SIUNION sival;
X
X    sign = z1.sign;
X    if (sign && ziseven(z2))
X        math_error("Even root of negative number");
X    if (ziszero(z2) || zisneg(z2))
X        math_error("Non-positive root");
X    if (ziszero(z1)) {    /* root of zero */
X        *dest = _zero_;
X        return;
X    }
X    if (zisunit(z2)) {    /* first root */
X        zcopy(z1, dest);
X        return;
X    }
X    if (zisbig(z2)) {    /* humongous root */
X        *dest = _one_;
X        dest->sign = (HALF)sign;
X        return;
X    }
X    k = (zistiny(z2) ? z1tol(z2) : z2tol(z2));
X    highbit = zhighbit(z1);
X    if (highbit < k) {    /* too high a root */
X        *dest = _one_;
X        dest->sign = (HALF)sign;
X        return;
X    }
X    sival.ivalue = k - 1;
X    k1.v = &sival.silow;
X    k1.len = 1 + (sival.sihigh != 0);
X    k1.sign = 0;
X    z1.sign = 0;
X    /*
X     * Allocate the numbers to use for the main loop.
X     * The size and high bits of the final result are correctly set here.
X     * Notice that the remainder of the test value is rubbish, but this
SHAR_EOF
echo "End of part 11"
echo "File calc2.9.0/zfunc.c is continued in part 12"
echo "12" > s2_seq_.tmp
exit 0