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

Text File | 1993-12-07 | 59.1 KB | 2,521 lines

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

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

#!/bin/sh
# this is part 6 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/lint.sed continued
#
CurArch=6
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/lint.sed"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/lint.sed
X/: warning: possible pointer alignment problem$/d
X/^Lint pass[0-9][0-9]*:$/d
X/^[a-zA-Z][a-zA-Z0-9_-]*\.c:[     ]*$/d
X/^addglobal, arg\. 2 used inconsistently[     ]/d
X/^addopptr, arg\. 2 used inconsistently[     ]/d
X/^codegen\.c([0-9]*):getassignment returns value which is sometimes ignored$/d
X/^errno used([     ]*func\.c([0-9]*)[     ]*), but not defined$/d
X/^exit value declared inconsistently[     ]/d
X/^fclose returns value which is sometimes ignored$/d
X/^fflush returns value which is always ignored$/d
X/^fprintf returns value which is always ignored$/d
X/^fputc returns value which is always ignored$/d
X/^fputs returns value which is always ignored$/d
X/^free, arg\. 1 used inconsistently[     ]/d
X/^geteuid value declared inconsistently[     ]/d
X/^geteuid value used inconsistently[     ]/d
X/^getpwuid, arg\. 1 used inconsistently[     ]/d
X/^malloc, arg\. 1 used inconsistently[     ]/d
X/^math_setdigits returns value which is always ignored$/d
X/^math_setmode returns value which is sometimes ignored$/d
X/^memcpy returns value which is always ignored$/d
X/^memcpy value declared inconsistently[     ]/d
X/^memcpy, arg\. [1-3] used inconsistently[     ]/d
X/^memset value declared inconsistently[     ]/d
X/^printf returns value which is always ignored$/d
X/^putc returns value which is always ignored$/d
X/^qcfappr, arg\. 2 used inconsistently[     ]/d
X/^realloc, arg\. [1-2] used inconsistently[     ]/d
X/^sprintf returns value which is always ignored/d
X/^strcat returns value which is always ignored/d
X/^strcpy returns value which is always ignored/d
X/^strncpy returns value which is always ignored/d
X/^strncpy, arg\. [1-3] used inconsistently[     ]/d
X/^system returns value which is always ignored/d
X/^times returns value which is always ignored/d
X/^vsprintf returns value which is always ignored/d
SHAR_EOF
echo "File calc2.9.0/lint.sed is complete"
chmod 0644 calc2.9.0/lint.sed || echo "restore of calc2.9.0/lint.sed fails"
set `wc -c calc2.9.0/lint.sed`;Sum=$1
if test "$Sum" != "1807"
then echo original size 1807, current size $Sum;fi
echo "x - extracting calc2.9.0/listfunc.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/listfunc.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 * List handling routines.
X * Lists can be composed of any types of values, mixed if desired.
X * Lists are doubly linked so that elements can be inserted or
X * deleted efficiently at any point in the list.  A pointer is
X * kept to the most recently indexed element so that sequential
X * accesses are fast.
X */
X
X#include "value.h"
X
X
Xstatic LISTELEM *elemalloc MATH_PROTO((void));
Xstatic LISTELEM *listelement MATH_PROTO((LIST *lp, long index));
Xstatic void elemfree MATH_PROTO((LISTELEM *ep));
Xstatic void removelistelement MATH_PROTO((LIST *lp, LISTELEM *ep));
X
X
X/*
X * Free lists for list headers and list elements.
X */
Xstatic FREELIST    headerfreelist = {
X    sizeof(LIST),        /* size of list header */
X    20            /* number of free headers to keep */
X};
X
Xstatic FREELIST elementfreelist = {
X    sizeof(LISTELEM),    /* size of list element */
X    1000            /* number of free list elements to keep */
X};
X
X
X/*
X * Insert an element before the first element of a list.
X */
Xvoid
Xinsertlistfirst(lp, vp)
X    LIST *lp;        /* list to put element onto */
X    VALUE *vp;        /* value to be inserted */
X{
X    LISTELEM *ep;        /* list element */
X
X    ep = elemalloc();
X    copyvalue(vp, &ep->e_value);
X    if (lp->l_count == 0)
X        lp->l_last = ep;
X    else {
X        lp->l_cacheindex++;
X        lp->l_first->e_prev = ep;
X        ep->e_next = lp->l_first;
X    }
X    lp->l_first = ep;
X    lp->l_count++;
X}
X
X
X/*
X * Insert an element after the last element of a list.
X */
Xvoid
Xinsertlistlast(lp, vp)
X    LIST *lp;        /* list to put element onto */
X    VALUE *vp;        /* value to be inserted */
X{
X    LISTELEM *ep;        /* list element */
X
X    ep = elemalloc();
X    copyvalue(vp, &ep->e_value);
X    if (lp->l_count == 0)
X        lp->l_first = ep;
X    else {
X        lp->l_last->e_next = ep;
X        ep->e_prev = lp->l_last;
X    }
X    lp->l_last = ep;
X    lp->l_count++;
X}
X
X
X/*
X * Insert an element into the middle of list at the given index (zero based).
X * The specified index will select the new element, so existing elements
X * at or beyond the index will be shifted down one position.  It is legal
X * to specify an index which is right at the end of the list, in which
X * case the element is appended to the list.
X */
Xvoid
Xinsertlistmiddle(lp, index, vp)
X    LIST *lp;        /* list to put element onto */
X    long index;        /* element number to insert in front of */
X    VALUE *vp;        /* value to be inserted */
X{
X    LISTELEM *ep;        /* list element */
X    LISTELEM *oldep;    /* old list element at desired index */
X
X    if (index == 0) {
X        insertlistfirst(lp, vp);
X        return;
X    }
X    if (index == lp->l_count) {
X        insertlistlast(lp, vp);
X        return;
X    }
X    oldep = NULL;
X    if ((index >= 0) && (index < lp->l_count))
X        oldep = listelement(lp, index);
X    if (oldep == NULL)
X        math_error("Index out of bounds for list insertion");
X    ep = elemalloc();
X    copyvalue(vp, &ep->e_value);
X    ep->e_next = oldep;
X    ep->e_prev = oldep->e_prev;
X    ep->e_prev->e_next = ep;
X    oldep->e_prev = ep;
X    lp->l_cache = ep;
X    lp->l_cacheindex = index;
X    lp->l_count++;
X}
X
X
X/*
X * Remove the first element from a list, returning its value.
X * Returns the null value if no more elements exist.
X */
Xvoid
Xremovelistfirst(lp, vp)
X    LIST *lp;        /* list to have element removed */
X    VALUE *vp;        /* location of the value */
X{
X    if (lp->l_count == 0) {
X        vp->v_type = V_NULL;
X        return;
X    }
X    *vp = lp->l_first->e_value;
X    lp->l_first->e_value.v_type = V_NULL;
X    removelistelement(lp, lp->l_first);
X}
X
X
X/*
X * Remove the last element from a list, returning its value.
X * Returns the null value if no more elements exist.
X */
Xvoid
Xremovelistlast(lp, vp)
X    LIST *lp;        /* list to have element removed */
X    VALUE *vp;        /* location of the value */
X{
X    if (lp->l_count == 0) {
X        vp->v_type = V_NULL;
X        return;
X    }
X    *vp = lp->l_last->e_value;
X    lp->l_last->e_value.v_type = V_NULL;
X    removelistelement(lp, lp->l_last);
X}
X
X
X/*
X * Remove the element with the given index from a list, returning its value.
X */
Xvoid
Xremovelistmiddle(lp, index, vp)
X    LIST *lp;        /* list to have element removed */
X    long index;        /* list element to be removed */
X    VALUE *vp;        /* location of the value */
X{
X    LISTELEM *ep;        /* element being removed */
X
X    ep = NULL;
X    if ((index >= 0) && (index < lp->l_count))
X        ep = listelement(lp, index);
X    if (ep == NULL)
X        math_error("Index out of bounds for list deletion");
X    *vp = ep->e_value;
X    ep->e_value.v_type = V_NULL;
X    removelistelement(lp, ep);
X}
X
X
X/*
X * Remove an arbitrary element from a list.
X * The value contained in the element is freed.
X */
Xstatic void
Xremovelistelement(lp, ep)
X    register LIST *lp;        /* list header */
X    register LISTELEM *ep;        /* list element to remove */
X{
X    if ((ep == lp->l_cache) || ((ep != lp->l_first) && (ep != lp->l_last)))
X        lp->l_cache = NULL;
X    if (ep->e_next)
X        ep->e_next->e_prev = ep->e_prev;
X    if (ep->e_prev)
X        ep->e_prev->e_next = ep->e_next;
X    if (ep == lp->l_first) {
X        lp->l_first = ep->e_next;
X        lp->l_cacheindex--;
X    }
X    if (ep == lp->l_last)
X        lp->l_last = ep->e_prev;
X    lp->l_count--;
X    elemfree(ep);
X}
X
X
X/*
X * Search a list for the specified value starting at the specified index.
X * Returns the element number (zero based) of the found value, or -1 if
X * the value was not found.
X */
Xlong
Xlistsearch(lp, vp, index)
X    LIST *lp;
X    VALUE *vp;
X    long index;
X{
X    register LISTELEM *ep;
X
X    if (index < 0)
X        index = 0;
X    ep = listelement(lp, index);
X    while (ep) {
X        if (!comparevalue(&ep->e_value, vp)) {
X            lp->l_cache = ep;
X            lp->l_cacheindex = index;
X            return index;
X        }
X        ep = ep->e_next;
X        index++;
X    }
X    return -1;
X}
X
X
X/*
X * Search a list backwards for the specified value starting at the
X * specified index.  Returns the element number (zero based) of the
X * found value, or -1 if the value was not found.
X */
Xlong
Xlistrsearch(lp, vp, index)
X    LIST *lp;
X    VALUE *vp;
X    long index;
X{
X    register LISTELEM *ep;
X
X    if (index >= lp->l_count)
X        index = lp->l_count - 1;
X    ep = listelement(lp, index);
X    while (ep) {
X        if (!comparevalue(&ep->e_value, vp)) {
X            lp->l_cache = ep;
X            lp->l_cacheindex = index;
X            return index;
X        }
X        ep = ep->e_prev;
X        index--;
X    }
X    return -1;
X}
X
X
X/*
X * Index into a list and return the address for the value corresponding
X * to that index.  Returns NULL if the element does not exist.
X */
XVALUE *
Xlistfindex(lp, index)
X    LIST *lp;        /* list to index into */
X    long index;        /* index of desired element */
X{
X    LISTELEM *ep;
X
X    ep = listelement(lp, index);
X    if (ep == NULL)
X        return NULL;
X    return &ep->e_value;
X}
X
X
X/*
X * Return the element at a specified index number of a list.
X * The list is indexed starting at zero, and negative indices
X * indicate to index from the end of the list.  This routine finds
X * the element by chaining through the list from the closest one
X * of the first, last, and cached elements.  Returns NULL if the
X * element does not exist.
X */
Xstatic LISTELEM *
Xlistelement(lp, index)
X    register LIST *lp;    /* list to index into */
X    long index;        /* index of desired element */
X{
X    register LISTELEM *ep;    /* current list element */
X    long dist;        /* distance to element */
X    long temp;        /* temporary distance */
X    BOOL forward;        /* TRUE if need to walk forwards */
X
X    if (index < 0)
X        index += lp->l_count;
X    if ((index < 0) || (index >= lp->l_count))
X        return NULL;
X    /*
X     * Check quick special cases first.
X     */
X    if (index == 0)
X        return lp->l_first;
X    if (index == 1)
X        return lp->l_first->e_next;
X    if (index == lp->l_count - 1)
X        return lp->l_last;
X    if ((index == lp->l_cacheindex) && lp->l_cache)
X        return lp->l_cache;
X    /*
X     * Calculate whether it is better to go forwards from
X     * the first element or backwards from the last element.
X     */
X    forward = ((index * 2) <= lp->l_count);
X    if (forward) {
X        dist = index;
X        ep = lp->l_first;
X    } else {
X        dist = (lp->l_count - 1) - index;
X        ep = lp->l_last;
X    }
X    /*
X     * Now see if we have a cached element and if so, whether or
X     * not the distance from it is better than the above distance.
X     */
X    if (lp->l_cache) {
X        temp = index - lp->l_cacheindex;
X        if ((temp >= 0) && (temp < dist)) {
X            dist = temp;
X            ep = lp->l_cache;
X            forward = TRUE;
X        }
X        if ((temp < 0) && (-temp < dist)) {
X            dist = -temp;
X            ep = lp->l_cache;
X            forward = FALSE;
X        }
X    }
X    /*
X     * Now walk forwards or backwards from the selected element
X     * until we reach the correct element.  Cache the location of
X     * the found element for future use.
X     */
X    if (forward) {
X        while (dist-- > 0)
X            ep = ep->e_next;
X    } else {
X        while (dist-- > 0)
X            ep = ep->e_prev;
X    }
X    lp->l_cache = ep;
X    lp->l_cacheindex = index;
X    return ep;
X}
X
X
X/*
X * Compare two lists to see if they are identical.
X * Returns TRUE if they are different.
X */
XBOOL
Xlistcmp(lp1, lp2)
X    LIST *lp1, *lp2;
X{
X    LISTELEM *e1, *e2;
X    long count;
X
X    if (lp1 == lp2)
X        return FALSE;
X    if (lp1->l_count != lp2->l_count)
X        return TRUE;
X    e1 = lp1->l_first;
X    e2 = lp2->l_first;
X    count = lp1->l_count;
X    while (count-- > 0) {
X        if (comparevalue(&e1->e_value, &e2->e_value))
X            return TRUE;
X        e1 = e1->e_next;
X        e2 = e2->e_next;
X    }
X    return FALSE;
X}
X
X
X/*
X * Copy a list
X */
XLIST *
Xlistcopy(oldlp)
X    LIST *oldlp;
X{
X    LIST *lp;
X    LISTELEM *oldep;
X
X    lp = listalloc();
X    oldep = oldlp->l_first;
X    while (oldep) {
X        insertlistlast(lp, &oldep->e_value);
X        oldep = oldep->e_next;
X    }
X    return lp;
X}
X
X
X/*
X * Allocate an element for a list.
X */
Xstatic LISTELEM *
Xelemalloc()
X{
X    LISTELEM *ep;
X
X    ep = (LISTELEM *) allocitem(&elementfreelist);
X    if (ep == NULL)
X        math_error("Cannot allocate list element");
X    ep->e_next = NULL;
X    ep->e_prev = NULL;
X    ep->e_value.v_type = V_NULL;
X    return ep;
X}
X
X
X/*
X * Free a list element, along with any contained value.
X */
Xstatic void
Xelemfree(ep)
X    LISTELEM *ep;
X{
X    if (ep->e_value.v_type != V_NULL)
X        freevalue(&ep->e_value);
X    freeitem(&elementfreelist, (FREEITEM *) ep);
X}
X
X
X/*
X * Allocate a new list header.
X */
XLIST *
Xlistalloc()
X{
X    register LIST *lp;
X
X    lp = (LIST *) allocitem(&headerfreelist);
X    if (lp == NULL)
X        math_error("Cannot allocate list header");
X    lp->l_first = NULL;
X    lp->l_last = NULL;
X    lp->l_cache = NULL;
X    lp->l_cacheindex = 0;
X    lp->l_count = 0;
X    return lp;
X}
X
X
X/*
X * Free a list header, along with all of its list elements.
X */
Xvoid
Xlistfree(lp)
X    register LIST *lp;
X{
X    register LISTELEM *ep;
X
X    while (lp->l_count-- > 0) {
X        ep = lp->l_first;
X        lp->l_first = ep->e_next;
X        elemfree(ep);
X    }
X    freeitem(&headerfreelist, (FREEITEM *) lp);
X}
X
X
X/*
X * Print out a list along with the specified number of its elements.
X * The elements are printed out in shortened form.
X */
Xvoid
Xlistprint(lp, max_print)
X    LIST *lp;
X    long max_print;
X{
X    long count;
X    long index;
X    LISTELEM *ep;
X
X    if (max_print > lp->l_count)
X        max_print = lp->l_count;
X    count = 0;
X    ep = lp->l_first;
X    index = lp->l_count;
X    while (index-- > 0) {
X        if ((ep->e_value.v_type != V_NUM) ||
X            (!qiszero(ep->e_value.v_num)))
X                count++;
X        ep = ep->e_next;
X    }
X    if (max_print > 0)
X        math_str("\n");
X    math_fmt("list (%ld element%s, %ld nonzero)", lp->l_count,
X        ((lp->l_count == 1) ? "" : "s"), count);
X    if (max_print <= 0)
X        return;
X
X    /*
X     * Walk through the first few list elements, printing their
X     * value in short and unambiguous format.
X     */
X    math_str(":\n");
X    ep = lp->l_first;
X    for (index = 0; index < max_print; index++) {
X        math_fmt("  [[%ld]] = ", index);
X        printvalue(&ep->e_value, PRINT_SHORT | PRINT_UNAMBIG);
X        math_str("\n");
X        ep = ep->e_next;
X    }
X    if (max_print < lp->l_count)
X        math_str("  ...\n");
X}
X
X
X/*
X * Return a trivial hash value for a list.
X */
XHASH
Xlisthash(lp)
X    LIST *lp;
X{
X    HASH hash;
X
X    hash = lp->l_count * 600011;
X    if (lp->l_count > 0)
X        hash = hash * 600043 + hashvalue(&lp->l_first->e_value);
X    if (lp->l_count > 1)
X        hash = hash * 600053 + hashvalue(&lp->l_last->e_value);
X    return hash;
X}
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/listfunc.c || echo "restore of calc2.9.0/listfunc.c fails"
set `wc -c calc2.9.0/listfunc.c`;Sum=$1
if test "$Sum" != "11559"
then echo original size 11559, current size $Sum;fi
echo "x - extracting calc2.9.0/matfunc.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/matfunc.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 rational arithmetic matrix functions.
X * Matrices can contain arbitrary types of elements.
X */
X
X#include "value.h"
X
X
Xstatic void matswaprow MATH_PROTO((MATRIX *m, long r1, long r2));
Xstatic void matsubrow MATH_PROTO((MATRIX *m, long oprow, long baserow,
X    VALUE *mulval));
Xstatic void matmulrow MATH_PROTO((MATRIX *m, long row, VALUE *mulval));
Xstatic MATRIX *matident MATH_PROTO((MATRIX *m));
X
X
X
X/*
X * Add two compatible matrices.
X */
XMATRIX *
Xmatadd(m1, m2)
X    MATRIX *m1, *m2;
X{
X    int dim;
X
X    long min1, min2, max1, max2, index;
X    VALUE *v1, *v2, *vres;
X    MATRIX *res;
X    MATRIX tmp;
X
X    if (m1->m_dim != m2->m_dim)
X        math_error("Incompatible matrix dimensions for add");
X    tmp.m_dim = m1->m_dim;
X    tmp.m_size = m1->m_size;
X    for (dim = 0; dim < m1->m_dim; dim++) {
X        min1 = m1->m_min[dim];
X        max1 = m1->m_max[dim];
X        min2 = m2->m_min[dim];
X        max2 = m2->m_max[dim];
X        if ((min1 && min2 && (min1 != min2)) || ((max1-min1) != (max2-min2)))
X            math_error("Incompatible matrix bounds for add");
X        tmp.m_min[dim] = (min1 ? min1 : min2);
X        tmp.m_max[dim] = tmp.m_min[dim] + (max1 - min1);
X    }
X    res = matalloc(m1->m_size);
X    *res = tmp;
X    v1 = m1->m_table;
X    v2 = m2->m_table;
X    vres = res->m_table;
X    for (index = m1->m_size; index > 0; index--)
X        addvalue(v1++, v2++, vres++);
X    return res;
X}
X
X
X/*
X * Subtract two compatible matrices.
X */
XMATRIX *
Xmatsub(m1, m2)
X    MATRIX *m1, *m2;
X{
X    int dim;
X    long min1, min2, max1, max2, index;
X    VALUE *v1, *v2, *vres;
X    MATRIX *res;
X    MATRIX tmp;
X
X    if (m1->m_dim != m2->m_dim)
X        math_error("Incompatible matrix dimensions for sub");
X    tmp.m_dim = m1->m_dim;
X    tmp.m_size = m1->m_size;
X    for (dim = 0; dim < m1->m_dim; dim++) {
X        min1 = m1->m_min[dim];
X        max1 = m1->m_max[dim];
X        min2 = m2->m_min[dim];
X        max2 = m2->m_max[dim];
X        if ((min1 && min2 && (min1 != min2)) || ((max1-min1) != (max2-min2)))
X            math_error("Incompatible matrix bounds for sub");
X        tmp.m_min[dim] = (min1 ? min1 : min2);
X        tmp.m_max[dim] = tmp.m_min[dim] + (max1 - min1);
X    }
X    res = matalloc(m1->m_size);
X    *res = tmp;
X    v1 = m1->m_table;
X    v2 = m2->m_table;
X    vres = res->m_table;
X    for (index = m1->m_size; index > 0; index--)
X        subvalue(v1++, v2++, vres++);
X    return res;
X}
X
X
X/*
X * Produce the negative of a matrix.
X */
XMATRIX *
Xmatneg(m)
X    MATRIX *m;
X{
X    register VALUE *val, *vres;
X    long index;
X    MATRIX *res;
X
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        negvalue(val++, vres++);
X    return res;
X}
X
X
X/*
X * Multiply two compatible matrices.
X */
XMATRIX *
Xmatmul(m1, m2)
X    MATRIX *m1, *m2;
X{
X    register MATRIX *res;
X    long i1, i2, max1, max2, index, maxindex;
X    VALUE *v1, *v2;
X    VALUE sum, tmp1, tmp2;
X
X    if ((m1->m_dim != 2) || (m2->m_dim != 2))
X        math_error("Matrix dimension must be two for mul");
X    if ((m1->m_max[1] - m1->m_min[1]) != (m2->m_max[0] - m2->m_min[0]))
X        math_error("Incompatible bounds for matrix mul");
X    max1 = (m1->m_max[0] - m1->m_min[0] + 1);
X    max2 = (m2->m_max[1] - m2->m_min[1] + 1);
X    maxindex = (m1->m_max[1] - m1->m_min[1] + 1);
X    res = matalloc(max1 * max2);
X    res->m_dim = 2;
X    res->m_min[0] = m1->m_min[0];
X    res->m_max[0] = m1->m_max[0];
X    res->m_min[1] = m2->m_min[1];
X    res->m_max[1] = m2->m_max[1];
X    for (i1 = 0; i1 < max1; i1++) {
X        for (i2 = 0; i2 < max2; i2++) {
X            sum.v_num = qlink(&_qzero_);
X            sum.v_type = V_NUM;
X            v1 = &m1->m_table[i1 * maxindex];
X            v2 = &m2->m_table[i2];
X            for (index = 0; index < maxindex; index++) {
X                mulvalue(v1, v2, &tmp1);
X                addvalue(&sum, &tmp1, &tmp2);
X                freevalue(&tmp1);
X                freevalue(&sum);
X                sum = tmp2;
X                v1++;
X                v2 += max2;
X            }
X            index = (i1 * max2) + i2;
X            res->m_table[index] = sum;
X        }
X    }
X    return res;
X}
X
X
X/*
X * Square a matrix.
X */
XMATRIX *
Xmatsquare(m)
X    MATRIX *m;
X{
X    register MATRIX *res;
X    long i1, i2, max, index;
X    VALUE *v1, *v2;
X    VALUE sum, tmp1, tmp2;
X
X    if (m->m_dim != 2)
X        math_error("Matrix dimension must be two for square");
X    if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
X        math_error("Squaring non-square matrix");
X    max = (m->m_max[0] - m->m_min[0] + 1);
X    res = matalloc(max * max);
X    res->m_dim = 2;
X    res->m_min[0] = m->m_min[0];
X    res->m_max[0] = m->m_max[0];
X    res->m_min[1] = m->m_min[1];
X    res->m_max[1] = m->m_max[1];
X    for (i1 = 0; i1 < max; i1++) {
X        for (i2 = 0; i2 < max; i2++) {
X            sum.v_num = qlink(&_qzero_);
X            sum.v_type = V_NUM;
X            v1 = &m->m_table[i1 * max];
X            v2 = &m->m_table[i2];
X            for (index = 0; index < max; index++) {
X                mulvalue(v1, v2, &tmp1);
X                addvalue(&sum, &tmp1, &tmp2);
X                freevalue(&tmp1);
X                freevalue(&sum);
X                sum = tmp2;
X                v1++;
X                v2 += max;
X            }
X            index = (i1 * max) + i2;
X            res->m_table[index] = sum;
X        }
X    }
X    return res;
X}
X
X
X/*
X * Compute the result of raising a square matrix to an integer power.
X * Negative powers mean the positive power of the inverse.
X * Note: This calculation could someday be improved for large powers
X * by using the characteristic polynomial of the matrix.
X */
XMATRIX *
Xmatpowi(m, q)
X    MATRIX *m;        /* matrix to be raised */
X    NUMBER *q;        /* power to raise it to */
X{
X    MATRIX *res, *tmp;
X    long power;        /* power to raise to */
X    unsigned long bit;    /* current bit value */
X
X    if (m->m_dim != 2)
X        math_error("Matrix dimension must be two for power");
X    if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
X        math_error("Raising non-square matrix to a power");
X    if (qisfrac(q))
X        math_error("Raising matrix to non-integral power");
X    if (zisbig(q->num))
X        math_error("Raising matrix to very large power");
X    power = (zistiny(q->num) ? z1tol(q->num) : z2tol(q->num));
X    if (qisneg(q))
X        power = -power;
X    /*
X     * Handle some low powers specially
X     */
X    if ((power <= 4) && (power >= -2)) {
X        switch ((int) power) {
X            case 0:
X                return matident(m);
X            case 1:
X                return matcopy(m);
X            case -1:
X                return matinv(m);
X            case 2:
X                return matsquare(m);
X            case -2:
X                tmp = matinv(m);
X                res = matsquare(tmp);
X                matfree(tmp);
X                return res;
X            case 3:
X                tmp = matsquare(m);
X                res = matmul(m, tmp);
X                matfree(tmp);
X                return res;
X            case 4:
X                tmp = matsquare(m);
X                res = matsquare(tmp);
X                matfree(tmp);
X                return res;
X        }
X    }
X    if (power < 0) {
X        m = matinv(m);
X        power = -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    res = matsquare(m);
X    if (bit & power) {
X        tmp = matmul(res, m);
X        matfree(res);
X        res = tmp;
X    }
X    bit >>= 1L;
X    while (bit) {
X        tmp = matsquare(res);
X        matfree(res);
X        res = tmp;
X        if (bit & power) {
X            tmp = matmul(res, m);
X            matfree(res);
X            res = tmp;
X        }
X        bit >>= 1L;
X    }
X    if (qisneg(q))
X        matfree(m);
X    return res;
X}
X
X
X/*
X * Calculate the cross product of two one dimensional matrices each
X * with three components.
X *    m3 = matcross(m1, m2);
X */
XMATRIX *
Xmatcross(m1, m2)
X    MATRIX *m1, *m2;
X{
X    MATRIX *res;
X    VALUE *v1, *v2, *vr;
X    VALUE tmp1, tmp2;
X
X    if ((m1->m_dim != 1) || (m2->m_dim != 1))
X        math_error("Matrix not 1d for cross product");
X    if ((m1->m_size != 3) || (m2->m_size != 3))
X        math_error("Matrix not size 3 for cross product");
X    res = matalloc(3L);
X    res->m_dim = 1;
X    res->m_min[0] = 0;
X    res->m_max[0] = 2;
X    v1 = m1->m_table;
X    v2 = m2->m_table;
X    vr = res->m_table;
X    mulvalue(v1 + 1, v2 + 2, &tmp1);
X    mulvalue(v1 + 2, v2 + 1, &tmp2);
X    subvalue(&tmp1, &tmp2, vr + 0);
X    freevalue(&tmp1);
X    freevalue(&tmp2);
X    mulvalue(v1 + 2, v2 + 0, &tmp1);
X    mulvalue(v1 + 0, v2 + 2, &tmp2);
X    subvalue(&tmp1, &tmp2, vr + 1);
X    freevalue(&tmp1);
X    freevalue(&tmp2);
X    mulvalue(v1 + 0, v2 + 1, &tmp1);
X    mulvalue(v1 + 1, v2 + 0, &tmp2);
X    subvalue(&tmp1, &tmp2, vr + 2);
X    freevalue(&tmp1);
X    freevalue(&tmp2);
X    return res;
X}
X
X
X/*
X * Return the dot product of two matrices.
X *    result = matdot(m1, m2);
X */
XVALUE
Xmatdot(m1, m2)
X    MATRIX *m1, *m2;
X{
X    VALUE *v1, *v2;
X    VALUE result, tmp1, tmp2;
X    long len;
X
X    if ((m1->m_dim != 1) || (m2->m_dim != 1))
X        math_error("Matrix not 1d for dot product");
X    if (m1->m_size != m2->m_size)
X        math_error("Incompatible matrix sizes for dot product");
X    v1 = m1->m_table;
X    v2 = m2->m_table;
X    mulvalue(v1, v2, &result);
X    len = m1->m_size;
X    while (--len > 0) {
X        mulvalue(++v1, ++v2, &tmp1);
X        addvalue(&result, &tmp1, &tmp2);
X        freevalue(&tmp1);
X        freevalue(&result);
X        result = tmp2;
X    }
X    return result;
X}
X
X
X/*
X * Scale the elements of a matrix by a specified power of two.
X */
XMATRIX *
Xmatscale(m, n)
X    MATRIX *m;        /* matrix to be scaled */
X    long n;
X{
X    register VALUE *val, *vres;
X    VALUE num;
X    long index;
X    MATRIX *res;        /* resulting matrix */
X
X    if (n == 0)
X        return matcopy(m);
X    num.v_type = V_NUM;
X    num.v_num = itoq(n);
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        scalevalue(val++, &num, vres++);
X    qfree(num.v_num);
X    return res;
X}
X
X
X/*
X * Shift the elements of a matrix by the specified number of bits.
X * Positive shift means leftwards, negative shift rightwards.
X */
XMATRIX *
Xmatshift(m, n)
X    MATRIX *m;        /* matrix to be scaled */
X    long n;
X{
X    register VALUE *val, *vres;
X    VALUE num;
X    long index;
X    MATRIX *res;        /* resulting matrix */
X
X    if (n == 0)
X        return matcopy(m);
X    num.v_type = V_NUM;
X    num.v_num = itoq(n);
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        shiftvalue(val++, &num, FALSE, vres++);
X    qfree(num.v_num);
X    return res;
X}
X
X
X/*
X * Multiply the elements of a matrix by a specified value.
X */
XMATRIX *
Xmatmulval(m, vp)
X    MATRIX *m;        /* matrix to be multiplied */
X    VALUE *vp;        /* value to multiply by */
X{
X    register VALUE *val, *vres;
X    long index;
X    MATRIX *res;
X
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        mulvalue(val++, vp, vres++);
X    return res;
X}
X
X
X/*
X * Divide the elements of a matrix by a specified value, keeping
X * only the integer quotient.
X */
XMATRIX *
Xmatquoval(m, vp)
X    MATRIX *m;        /* matrix to be divided */
X    VALUE *vp;        /* value to divide by */
X{
X    register VALUE *val, *vres;
X    long index;
X    MATRIX *res;
X
X    if ((vp->v_type == V_NUM) && qiszero(vp->v_num))
X        math_error("Division by zero");
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        quovalue(val++, vp, vres++);
X    return res;
X}
X
X
X/*
X * Divide the elements of a matrix by a specified value, keeping
X * only the remainder of the division.
X */
XMATRIX *
Xmatmodval(m, vp)
X    MATRIX *m;        /* matrix to be divided */
X    VALUE *vp;        /* value to divide by */
X{
X    register VALUE *val, *vres;
X    long index;
X    MATRIX *res;
X
X    if ((vp->v_type == V_NUM) && qiszero(vp->v_num))
X        math_error("Division by zero");
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        modvalue(val++, vp, vres++);
X    return res;
X}
X
X
XMATRIX *
Xmattrans(m)
X    MATRIX *m;            /* matrix to be transposed */
X{
X    register VALUE *v1, *v2;    /* current values */
X    long rows, cols;        /* rows and columns in new matrix */
X    long row, col;            /* current row and column */
X    MATRIX *res;
X
X    if (m->m_dim != 2)
X        math_error("Matrix dimension must be two for transpose");
X    res = matalloc(m->m_size);
X    res->m_dim = 2;
X    res->m_min[0] = m->m_min[1];
X    res->m_max[0] = m->m_max[1];
X    res->m_min[1] = m->m_min[0];
X    res->m_max[1] = m->m_max[0];
X    rows = (m->m_max[1] - m->m_min[1] + 1);
X    cols = (m->m_max[0] - m->m_min[0] + 1);
X    v1 = res->m_table;
X    for (row = 0; row < rows; row++) {
X        v2 = &m->m_table[row];
X        for (col = 0; col < cols; col++) {
X            copyvalue(v2, v1);
X            v1++;
X            v2 += rows;
X        }
X    }
X    return res;
X}
X
X
X/*
X * Produce a matrix with values all of which are conjugated.
X */
XMATRIX *
Xmatconj(m)
X    MATRIX *m;
X{
X    register VALUE *val, *vres;
X    long index;
X    MATRIX *res;
X
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        conjvalue(val++, vres++);
X    return res;
X}
X
X
X/*
X * Produce a matrix with values all of which have been rounded to the
X * specified number of decimal places.
X */
XMATRIX *
Xmatround(m, places)
X    MATRIX *m;
X    long places;
X{
X    register VALUE *val, *vres;
X    VALUE tmp;
X    long index;
X    MATRIX *res;
X
X    if (places < 0)
X        math_error("Negative number of places for matround");
X    res = matalloc(m->m_size);
X    *res = *m;
X    tmp.v_type = V_INT;
X    tmp.v_int = places;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        roundvalue(val++, &tmp, vres++);
X    return res;
X}
X
X
X/*
X * Produce a matrix with values all of which have been rounded to the
X * specified number of binary places.
X */
XMATRIX *
Xmatbround(m, places)
X    MATRIX *m;
X    long places;
X{
X    register VALUE *val, *vres;
X    VALUE tmp;
X    long index;
X    MATRIX *res;
X
X    if (places < 0)
X        math_error("Negative number of places for matbround");
X    res = matalloc(m->m_size);
X    *res = *m;
X    tmp.v_type = V_INT;
X    tmp.v_int = places;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        broundvalue(val++, &tmp, vres++);
X    return res;
X}
X
X
X/*
X * Produce a matrix with values all of which have been truncated to integers.
X */
XMATRIX *
Xmatint(m)
X    MATRIX *m;
X{
X    register VALUE *val, *vres;
X    long index;
X    MATRIX *res;
X
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        intvalue(val++, vres++);
X    return res;
X}
X
X
X/*
X * Produce a matrix with values all of which have only the fraction part left.
X */
XMATRIX *
Xmatfrac(m)
X    MATRIX *m;
X{
X    register VALUE *val, *vres;
X    long index;
X    MATRIX *res;
X
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = m->m_table;
X    vres = res->m_table;
X    for (index = m->m_size; index > 0; index--)
X        fracvalue(val++, vres++);
X    return res;
X}
X
X
X/*
X * Index a matrix normally by the specified set of index values.
X * Returns the address of the matrix element if it is valid, or generates
X * an error if the index values are out of range.  The create flag is TRUE
X * if the element is to be written, but this is ignored here.
X */
X/*ARGSUSED*/
XVALUE *
Xmatindex(mp, create, dim, indices)
X    MATRIX *mp;
X    BOOL create;
X    long dim;        /* dimension of the indexing */
X    VALUE *indices;        /* table of values being indexed by */
X{
X    NUMBER *q;        /* index value */
X    long index;        /* index value as an integer */
X    long offset;        /* current offset into array */
X    int i;            /* loop counter */
X
X    if ((dim <= 0) || (dim > MAXDIM))
X        math_error("Bad dimension %ld for matrix", dim);
X    if (mp->m_dim != dim)
X        math_error("Indexing a %ldd matrix as a %ldd matrix", mp->m_dim, dim);
X    offset = 0;
X    for (i = 0; i < dim; i++) {
X        if (indices->v_type != V_NUM)
X            math_error("Non-numeric index for matrix");
X        q = indices->v_num;
X        if (qisfrac(q))
X            math_error("Non-integral index for matrix");
X        index = qtoi(q);
X        if (zisbig(q->num) || (index < mp->m_min[i]) || (index > mp->m_max[i]))
X            math_error("Index out of bounds for matrix");
X        offset *= (mp->m_max[i] - mp->m_min[i] + 1);
X        offset += (index - mp->m_min[i]);
X        indices++;
X    }
X    return mp->m_table + offset;
X}
X
X
X/*
X * Search a matrix for the specified value, starting with the specified index.
X * Returns the index of the found value, or -1 if the value was not found.
X */
Xlong
Xmatsearch(m, vp, index)
X    MATRIX *m;
X    VALUE *vp;
X    long index;
X{
X    register VALUE *val;
X
X    if (index < 0)
X        index = 0;
X    val = &m->m_table[index];
X    while (index < m->m_size) {
X        if (!comparevalue(vp, val))
X            return index;
X        index++;
X        val++;
X    }
X    return -1;
X}
X
X
X/*
X * Search a matrix backwards for the specified value, starting with the
X * specified index.  Returns the index of the found value, or -1 if the
X * value was not found.
X */
Xlong
Xmatrsearch(m, vp, index)
X    MATRIX *m;
X    VALUE *vp;
X    long index;
X{
X    register VALUE *val;
X
X    if (index >= m->m_size)
X        index = m->m_size - 1;
X    val = &m->m_table[index];
X    while (index >= 0) {
X        if (!comparevalue(vp, val))
X            return index;
X        index--;
X        val--;
X    }
X    return -1;
X}
X
X
X/*
X * Fill all of the elements of a matrix with one of two specified values.
X * All entries are filled with the first specified value, except that if
X * the matrix is square and the second value pointer is non-NULL, then
X * all diagonal entries are filled with the second value.  This routine
X * affects the supplied matrix directly, and doesn't return a copy.
X */
Xvoid
Xmatfill(m, v1, v2)
X    MATRIX *m;        /* matrix to be filled */
X    VALUE *v1;        /* value to fill most of matrix with */
X    VALUE *v2;        /* value for diagonal entries (or NULL) */
X{
X    register VALUE *val;
X    long row, col;
X    long rows;
X    long index;
X
X    if (v2 && ((m->m_dim != 2) ||
X        ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))))
X            math_error("Filling diagonals of non-square matrix");
X    val = m->m_table;
X    for (index = m->m_size; index > 0; index--)
X        freevalue(val++);
X    val = m->m_table;
X    if (v2 == NULL) {
X        for (index = m->m_size; index > 0; index--)
X            copyvalue(v1, val++);
X        return;
X    }
X    rows = m->m_max[0] - m->m_min[0] + 1;
X    for (row = 0; row < rows; row++) {
X        for (col = 0; col < rows; col++) {
X            copyvalue(((row != col) ? v1 : v2), val++);
X        }
X    }
X}
X
X
X/*
X * Set a copy of a square matrix to the identity matrix.
X */
Xstatic MATRIX *
Xmatident(m)
X    MATRIX *m;
X{
X    register VALUE *val;    /* current value */
X    long row, col;        /* current row and column */
X    long rows;        /* number of rows */
X    MATRIX *res;        /* resulting matrix */
X
X    if (m->m_dim != 2)
X        math_error("Matrix dimension must be two for setting to identity");
X    if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
X        math_error("Matrix must be square for setting to identity");
X    res = matalloc(m->m_size);
X    *res = *m;
X    val = res->m_table;
X    rows = (res->m_max[0] - res->m_min[0] + 1);
X    for (row = 0; row < rows; row++) {
X        for (col = 0; col < rows; col++) {
X            val->v_type = V_NUM;
X            val->v_num = ((row == col) ? qlink(&_qone_) : qlink(&_qzero_));
X            val++;
X        }
X    }
X    return res;
X}
X
X
X/*
X * Calculate the inverse of a matrix if it exists.
X * This is done by using transformations on the supplied matrix to convert
X * it to the identity matrix, and simultaneously applying the same set of
X * transformations to the identity matrix.
X */
XMATRIX *
Xmatinv(m)
X    MATRIX *m;
X{
X    MATRIX *res;        /* matrix to become the inverse */
X    long rows;        /* number of rows */
X    long cur;        /* current row being worked on */
X    long row, col;        /* temp row and column values */
X    VALUE *val;        /* current value in matrix*/
X    VALUE mulval;        /* value to multiply rows by */
X    VALUE tmpval;        /* temporary value */
X
X    if (m->m_dim != 2)
X        math_error("Matrix dimension must be two for inverse");
X    if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
X        math_error("Inverting non-square matrix");
X    /*
X     * Begin by creating the identity matrix with the same attributes.
X     */
X    res = matalloc(m->m_size);
X    *res = *m;
X    rows = (m->m_max[0] - m->m_min[0] + 1);
X    val = res->m_table;
X    for (row = 0; row < rows; row++) {
X        for (col = 0; col < rows; col++) {
X            if (row == col)
X                val->v_num = qlink(&_qone_);
X            else
X                val->v_num = qlink(&_qzero_);
X            val->v_type = V_NUM;
X            val++;
X        }
X    }
X    /*
X     * Now loop over each row, and eliminate all entries in the
X     * corresponding column by using row operations.  Do the same
X     * operations on the resulting matrix.  Copy the original matrix
X     * so that we don't destroy it.
X     */
X    m = matcopy(m);
X    for (cur = 0; cur < rows; cur++) {
X        /*
X         * Find the first nonzero value in the rest of the column
X         * downwards from [cur,cur].  If there is no such value, then
X         * the matrix is not invertible.  If the first nonzero entry
X         * is not the current row, then swap the two rows to make the
X         * current one nonzero.
X         */
X        row = cur;
X        val = &m->m_table[(row * rows) + row];
X        while (testvalue(val) == 0) {
X            if (++row >= rows) {
X                matfree(m);
X                matfree(res);
X                math_error("Matrix is not invertible");
X            }
X            val += rows;
X        }
X        invertvalue(val, &mulval);
X        if (row != cur) {
X            matswaprow(m, row, cur);
X            matswaprow(res, row, cur);
X        }
X        /*
X         * Now for every other nonzero entry in the current column, subtract
X         * the appropriate multiple of the current row to force that entry
X         * to become zero.
X         */
X        val = &m->m_table[cur];
X        for (row = 0; row < rows; row++, val += rows) {
X            if ((row == cur) || (testvalue(val) == 0))
X                continue;
X            mulvalue(val, &mulval, &tmpval);
X            matsubrow(m, row, cur, &tmpval);
X            matsubrow(res, row, cur, &tmpval);
X            freevalue(&tmpval);
X        }
X        freevalue(&mulval);
X    }
X    /*
X     * Now the original matrix has nonzero entries only on its main diagonal.
X     * Scale the rows of the result matrix by the inverse of those entries.
X     */
X    val = m->m_table;
X    for (row = 0; row < rows; row++) {
X        if ((val->v_type != V_NUM) || !qisone(val->v_num)) {
X            invertvalue(val, &mulval);
X            matmulrow(res, row, &mulval);
X            freevalue(&mulval);
X        }
X        val += (rows + 1);
X    }
X    matfree(m);
X    return res;
X}
X
X
X/*
X * Calculate the determinant of a square matrix.
X * This is done using row operations to create an upper-diagonal matrix.
X */
XVALUE
Xmatdet(m)
X    MATRIX *m;
X{
X    long rows;        /* number of rows */
X    long cur;        /* current row being worked on */
X    long row;        /* temp row values */
X    int neg;        /* whether to negate determinant */
X    VALUE *val;        /* current value */
X    VALUE mulval, tmpval;    /* other values */
X
X    if (m->m_dim != 2)
X        math_error("Matrix dimension must be two for determinant");
X    if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
X        math_error("Non-square matrix for determinant");
X    /*
X     * Loop over each row, and eliminate all lower entries in the
X     * corresponding column by using row operations.  Copy the original
X     * matrix so that we don't destroy it.
X     */
X    neg = 0;
X    m = matcopy(m);
X    rows = (m->m_max[0] - m->m_min[0] + 1);
X    for (cur = 0; cur < rows; cur++) {
X        /*
X         * Find the first nonzero value in the rest of the column
X         * downwards from [cur,cur].  If there is no such value, then
X         * the determinant is zero.  If the first nonzero entry is not
X         * the current row, then swap the two rows to make the current
X         * one nonzero, and remember that the determinant changes sign.
X         */
X        row = cur;
X        val = &m->m_table[(row * rows) + row];
X        while (testvalue(val) == 0) {
X            if (++row >= rows) {
X                matfree(m);
X                mulval.v_type = V_NUM;
X                mulval.v_num = qlink(&_qzero_);
X                return mulval;
X            }
X            val += rows;
X        }
X        invertvalue(val, &mulval);
X        if (row != cur) {
X            matswaprow(m, row, cur);
X            neg = !neg;
X        }
X        /*
X         * Now for every other nonzero entry lower down in the current column,
X         * subtract the appropriate multiple of the current row to force that
X         * entry to become zero.
X         */
X        row = cur + 1;
X        val = &m->m_table[(row * rows) + cur];
X        for (; row < rows; row++, val += rows) {
X            if (testvalue(val) == 0)
X                continue;
X            mulvalue(val, &mulval, &tmpval);
X            matsubrow(m, row, cur, &tmpval);
X            freevalue(&tmpval);
X        }
X        freevalue(&mulval);
X    }
X    /*
X     * Now the matrix is upper-diagonal, and the determinant is the
X     * product of the main diagonal entries, and is possibly negated.
X     */
X    val = m->m_table;
X    mulval.v_type = V_NUM;
X    mulval.v_num = qlink(&_qone_);
X    for (row = 0; row < rows; row++) {
X        mulvalue(&mulval, val, &tmpval);
X        freevalue(&mulval);
X        mulval = tmpval;
X        val += (rows + 1);
X    }
X    matfree(m);
X    if (neg) {
X        negvalue(&mulval, &tmpval);
X        freevalue(&mulval);
X        return tmpval;
X    }
X    return mulval;
X}
X
X
X/*
X * Local utility routine to swap two rows of a square matrix.
X * No checks are made to verify the legality of the arguments.
X */
Xstatic void
Xmatswaprow(m, r1, r2)
X    MATRIX *m;
X    long r1, r2;
X{
X    register VALUE *v1, *v2;
X    register long rows;
X    VALUE tmp;
X
X    if (r1 == r2)
X        return;
X    rows = (m->m_max[0] - m->m_min[0] + 1);
X    v1 = &m->m_table[r1 * rows];
X    v2 = &m->m_table[r2 * rows];
X    while (rows-- > 0) {
X        tmp = *v1;
X        *v1 = *v2;
X        *v2 = tmp;
X        v1++;
X        v2++;
X    }
X}
X
X
X/*
X * Local utility routine to subtract a multiple of one row to another one.
X * The row to be changed is oprow, the row to be subtracted is baserow.
X * No checks are made to verify the legality of the arguments.
X */
Xstatic void
Xmatsubrow(m, oprow, baserow, mulval)
X    MATRIX *m;
X    long oprow, baserow;
X    VALUE *mulval;
X{
X    register VALUE *vop, *vbase;
X    register long entries;
X    VALUE tmp1, tmp2;
X
X    entries = (m->m_max[0] - m->m_min[0] + 1);
X    vop = &m->m_table[oprow * entries];
X    vbase = &m->m_table[baserow * entries];
X    while (entries-- > 0) {
X        mulvalue(vbase, mulval, &tmp1);
X        subvalue(vop, &tmp1, &tmp2);
X        freevalue(&tmp1);
X        freevalue(vop);
X        *vop = tmp2;
X        vop++;
X        vbase++;
X    }
X}
X
X
X/*
X * Local utility routine to multiply a row by a specified number.
X * No checks are made to verify the legality of the arguments.
X */
Xstatic void
Xmatmulrow(m, row, mulval)
X    MATRIX *m;
X    long row;
X    VALUE *mulval;
X{
X    register VALUE *val;
X    register long rows;
X    VALUE tmp;
X
X    rows = (m->m_max[0] - m->m_min[0] + 1);
X    val = &m->m_table[row * rows];
X    while (rows-- > 0) {
X        mulvalue(val, mulval, &tmp);
X        freevalue(val);
X        *val = tmp;
X        val++;
X    }
X}
X
X
X/*
X * Make a full copy of a matrix.
X */
XMATRIX *
Xmatcopy(m)
X    MATRIX *m;
X{
X    MATRIX *res;
X    register VALUE *v1, *v2;
X    register long i;
X
X    res = matalloc(m->m_size);
X    *res = *m;
X    v1 = m->m_table;
X    v2 = res->m_table;
X    i = m->m_size;
X    while (i-- > 0) {
X        if (v1->v_type == V_NUM) {
X            v2->v_type = V_NUM;
X            v2->v_num = qlink(v1->v_num);
X        } else
X            copyvalue(v1, v2);
X        v1++;
X        v2++;
X    }
X    return res;
X}
X
X
X/*
X * Allocate a matrix with the specified number of elements.
X */
XMATRIX *
Xmatalloc(size)
X    long size;
X{
X    MATRIX *m;
X
X    m = (MATRIX *) malloc(matsize(size));
X    if (m == NULL)
X        math_error("Cannot get memory to allocate matrix of size %d", size);
X    m->m_size = size;
X    return m;
X}
X
X
X/*
X * Free a matrix, along with all of its element values.
X */
Xvoid
Xmatfree(m)
X    MATRIX *m;
X{
X    register VALUE *vp;
X    register long i;
X
X    vp = m->m_table;
X    i = m->m_size;
X    while (i-- > 0) {
X        if (vp->v_type == V_NUM) {
X            vp->v_type = V_NULL;
X            qfree(vp->v_num);
X        } else
X            freevalue(vp);
X        vp++;
X    }
X    free(m);
X}
X
X
X/*
X * Test whether a matrix has any nonzero values.
X * Returns TRUE if so.
X */
XBOOL
Xmattest(m)
X    MATRIX *m;
X{
X    register VALUE *vp;
X    register long i;
X
X    vp = m->m_table;
X    i = m->m_size;
X    while (i-- > 0) {
X        if ((vp->v_type != V_NUM) || (!qiszero(vp->v_num)))
X            return TRUE;
X        vp++;
X    }
X    return FALSE;
X}
X
X
X/*
X * Test whether or not two matrices are equal.
X * Equality is determined by the shape and values of the matrices,
X * but not by their index bounds.  Returns TRUE if they differ.
X */
XBOOL
Xmatcmp(m1, m2)
X    MATRIX *m1, *m2;
X{
X    VALUE *v1, *v2;
X    long i;
X
X    if (m1 == m2)
X        return FALSE;
X    if ((m1->m_dim != m2->m_dim) || (m1->m_size != m2->m_size))
X        return TRUE;
X    for (i = 0; i < m1->m_dim; i++) {
X        if ((m1->m_max[i] - m1->m_min[i]) != (m2->m_max[i] - m2->m_min[i]))
X        return TRUE;
X    }
X    v1 = m1->m_table;
X    v2 = m2->m_table;
X    i = m1->m_size;
X    while (i-- > 0) {
X        if (comparevalue(v1++, v2++))
X            return TRUE;
X    }
X    return FALSE;
X}
X
X
X#if 0
X/*
X * Test whether or not a matrix is the identity matrix.
X * Returns TRUE if so.
X */
XBOOL
Xmatisident(m)
X    MATRIX *m;
X{
X    register VALUE *val;    /* current value */
X    long row, col;        /* row and column numbers */
X
X    if ((m->m_dim != 2) ||
X        ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])))
X            return FALSE;
X    val = m->m_table;
X    for (row = 0; row < m->m_size; row++) {
X        for (col = 0; col < m->m_size; col++) {
X            if (val->v_type != V_NUM)
X                return FALSE;
X            if (row == col) {
X                if (!qisone(val->v_num))
X                    return FALSE;
X            } else {
X                if (!qiszero(val->v_num))
X                    return FALSE;
X            }
X            val++;
X        }
X    }
X    return TRUE;
X}
X#endif
X
X
X/*
X * Return a trivial hash value for a matrix.
X */
XHASH
Xmathash(m)
X    MATRIX *m;
X{
X    HASH hash;
X    long fullsize;
X    long skip;
X    int i;
X    VALUE *vp;
X
X    hash = m->m_dim * 500009;
X    fullsize = 1;
X    for (i = m->m_dim - 1; i >= 0; i--) {
X        hash = hash * 500029 + m->m_max[i];
X        fullsize *= (m->m_max[i] - m->m_min[i] + 1);
X    }
X    hash = hash * 500041 + fullsize;
X    vp = m->m_table;
X    for (i = 0; ((i < fullsize) && (i < 16)); i++)
X        hash = hash * 500057 + hashvalue(vp++);
X    i = 16;
X    vp = &m->m_table[16];
X    skip = (fullsize / 11) + 1;
X    while (i < fullsize) {
X        hash = hash * 500069 + hashvalue(vp);
X        i += skip;
X        vp += skip;
X    }
X    return hash;
X}
X
X
X/*
X * Print a matrix and possibly few of its elements.
X * The argument supplied specifies how many elements to allow printing.
X * If any elements are printed, they are printed in short form.
X */
Xvoid
Xmatprint(m, max_print)
X    MATRIX *m;
X    long max_print;
X{
X    VALUE *vp;
X    long fullsize, count, index, num;
X    int dim, i;
X    char *msg;
X    long sizes[MAXDIM];
X
X    dim = m->m_dim;
X    fullsize = 1;
X    for (i = dim - 1; i >= 0; i--) {
X        sizes[i] = fullsize;
X        fullsize *= (m->m_max[i] - m->m_min[i] + 1);
X    }
X    msg = ((max_print > 0) ? "\nmat [" : "mat [");
X    for (i = 0; i < dim; i++) {
X        if (m->m_min[i])
X            math_fmt("%s%ld:%ld", msg, m->m_min[i], m->m_max[i]);
X        else
X            math_fmt("%s%ld", msg, m->m_max[i] + 1);
X        msg = ",";
X    }
X    if (max_print > fullsize)
X        max_print = fullsize;
X    vp = m->m_table;
X    count = 0;
X    for (index = 0; index < fullsize; index++) {
X        if ((vp->v_type != V_NUM) || !qiszero(vp->v_num))
X            count++;
X        vp++;
X    }
X    math_fmt("] (%ld element%s, %ld nonzero)",
X        fullsize, (fullsize == 1) ? "" : "s", count);
X    if (max_print <= 0)
X        return;
X
X    /*
X     * Now print the first few elements of the matrix in short
X     * and unambigous format.
X     */
X    math_str(":\n");
X    vp = m->m_table;
X    for (index = 0; index < max_print; index++) {
X        msg = "  [";
X        num = index;
X        for (i = 0; i < dim; i++) {
X            math_fmt("%s%ld", msg, m->m_min[i] + (num / sizes[i]));
X            num %= sizes[i];
X            msg = ",";
X        }
X        math_str("] = ");
X        printvalue(vp++, PRINT_SHORT | PRINT_UNAMBIG);
X        math_str("\n");
X    }
X    if (max_print < fullsize)
X        math_str("  ...\n");
X}
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/matfunc.c || echo "restore of calc2.9.0/matfunc.c fails"
set `wc -c calc2.9.0/matfunc.c`;Sum=$1
if test "$Sum" != "29631"
then echo original size 29631, current size $Sum;fi
echo "x - extracting calc2.9.0/obj.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/obj.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 * "Object" handling primatives.
X * This simply means that user-specified routines are called to perform
X * the indicated operations.
X */
X
X#include "calc.h"
X#include "opcodes.h"
X#include "func.h"
X#include "symbol.h"
X#include "string.h"
X
X
X/*
X * Types of values returned by calling object routines.
X */
X#define A_VALUE    0    /* returns arbitrary value */
X#define A_INT    1    /* returns integer value */
X#define A_UNDEF    2    /* returns no value */
X
X/*
X * Error handling actions for when the function is undefined.
X */
X#define E_NONE    0    /* no special action */
X#define E_PRINT    1    /* print element */
X#define E_CMP    2    /* compare two values */
X#define E_TEST    3    /* test value for nonzero */
X#define E_POW    4    /* call generic power routine */
X#define E_ONE    5    /* return number 1 */
X#define E_INC    6    /* increment by one */
X#define E_DEC    7    /* decrement by one */
X#define E_SQUARE 8    /* square value */
X
X
Xstatic struct objectinfo {
X    short args;    /* number of arguments */
X    short retval;    /* type of return value */
X    short error;    /* special action on errors */
X    char *name;    /* name of function to call */
X    char *comment;    /* useful comment if any */
X} objectinfo[] = {
X    1, A_UNDEF, E_PRINT, "print",    "print value, default prints elements",
X    1, A_VALUE, E_ONE,   "one",    "multiplicative identity, default is 1",
X    1, A_INT,   E_TEST,  "test",    "logical test (false,true => 0,1), default tests elements",
X    2, A_VALUE, E_NONE,  "add",    NULL,
X    2, A_VALUE, E_NONE,  "sub",    NULL,
X    1, A_VALUE, E_NONE,  "neg",    "negative",
X    2, A_VALUE, E_NONE,  "mul",    NULL,
X    2, A_VALUE, E_NONE,  "div",    "non-integral division",
X    1, A_VALUE, E_NONE,  "inv",    "multiplicative inverse",
X    2, A_VALUE, E_NONE,  "abs",    "absolute value within given error",
X    1, A_VALUE, E_NONE,  "norm",    "square of absolute value",
X    1, A_VALUE, E_NONE,  "conj",    "conjugate",
X    2, A_VALUE, E_POW,   "pow",    "integer power, default does multiply, square, inverse",
X    1, A_INT,   E_NONE,  "sgn",    "sign of value (-1, 0, 1)",
X    2, A_INT,   E_CMP,   "cmp",    "equality (equal,nonequal => 0,1), default tests elements",
X    2, A_INT,   E_NONE,  "rel",    "inequality (less,equal,greater => -1,0,1)",
X    2, A_VALUE, E_NONE,  "quo",    "integer quotient",
X    2, A_VALUE, E_NONE,  "mod",    "remainder of division",
X    1, A_VALUE, E_NONE,  "int",    "integer part",
X    1, A_VALUE, E_NONE,  "frac",    "fractional part",
X    1, A_VALUE, E_INC,   "inc",    "increment, default adds 1",
X    1, A_VALUE, E_DEC,   "dec",    "decrement, default subtracts 1",
X    1, A_VALUE, E_SQUARE,"square",    "default multiplies by itself",
X    2, A_VALUE, E_NONE,  "scale",    "multiply by power of 2",
X    2, A_VALUE, E_NONE,  "shift",    "shift left by n bits (right if negative)",
X    2, A_VALUE, E_NONE,  "round",    "round to given number of decimal places",
X    2, A_VALUE, E_NONE,  "bround",    "round to given number of binary places",
X    3, A_VALUE, E_NONE,  "root",    "root of value within given error",
X    2, A_VALUE, E_NONE,  "sqrt",    "square root within given error",
X    0, 0, 0, NULL
X};
X
X
Xstatic STRINGHEAD objectnames;    /* names of objects */
Xstatic STRINGHEAD elements;    /* element names for parts of objects */
Xstatic OBJECTACTIONS *objects[MAXOBJECTS]; /* table of actions for objects */
X
X
X/*
X * Free list of usual small objects.
X */
Xstatic FREELIST    freelist = {
X    sizeof(OBJECT),        /* size of typical objects */
X    100            /* number of free objects to keep */
X};
X
X
Xstatic VALUE objpowi MATH_PROTO((VALUE *vp, NUMBER *q));
Xstatic BOOL objtest MATH_PROTO((OBJECT *op));
Xstatic BOOL objcmp MATH_PROTO((OBJECT *op1, OBJECT *op2));
Xstatic void objprint MATH_PROTO((OBJECT *op));
X
X
X/*
X * Show all the routine names available for objects.
X */
Xvoid
Xshowobjfuncs()
X{
X    register struct objectinfo *oip;
X
X    printf("\nThe following object routines are definable.\n");
X    printf("Note: xx represents the actual object type name.\n\n");
X    printf("Name    Args    Comments\n");
X    for (oip = objectinfo; oip->name; oip++) {
X        printf("xx_%-8s %d    %s\n", oip->name, oip->args,
X            oip->comment ? oip->comment : "");
X    }
X    printf("\n");
X}
X
X
X/*
X * Call the appropriate user-defined routine to handle an object action.
X * Returns the value that the routine returned.
X */
XVALUE
Xobjcall(action, v1, v2, v3)
X    VALUE *v1, *v2, *v3;
X{
X    FUNC *fp;        /* function to call */
X    OBJECTACTIONS *oap;    /* object to call for */
X    struct objectinfo *oip;    /* information about action */
X    long index;        /* index of function (negative if undefined) */
X    VALUE val;        /* return value */
X    VALUE tmp;        /* temp value */
X    char name[SYMBOLSIZE+1];    /* full name of user routine to call */
X
X    if ((unsigned)action > OBJ_MAXFUNC)
X        math_error("Illegal action for object call");
X    oip = &objectinfo[action];
X    if (v1->v_type == V_OBJ)
X        oap = v1->v_obj->o_actions;
X    else if (v2->v_type == V_OBJ)
X        oap = v2->v_obj->o_actions;
X    else
X        math_error("Object routine called with non-object");
X    index = oap->actions[action];
X    if (index == 0) {
X        strcpy(name, oap->name);
X        strcat(name, "_");
X        strcat(name, oip->name);
X        index = adduserfunc(name);
X        oap->actions[action] = index;
X    }
X    fp = NULL;
X    if (index > 0)
X        fp = findfunc(index);
X    if (fp == NULL) {
X        switch (oip->error) {
X            case E_PRINT:
X                objprint(v1->v_obj);
X                val.v_type = V_NULL;
X                break;
X            case E_CMP:
X                val.v_type = V_INT;
X                if (v1->v_type != v2->v_type) {
X                    val.v_int = 1;
X                    return val;
X                }
X                val.v_int = objcmp(v1->v_obj, v2->v_obj);
X                break;
X            case E_TEST:
X                val.v_type = V_INT;
X                val.v_int = objtest(v1->v_obj);
X                break;
X            case E_POW:
X                if (v2->v_type != V_NUM)
X                    math_error("Non-real power");
X                val = objpowi(v1, v2->v_num);
X                break;
X            case E_ONE:
X                val.v_type = V_NUM;
X                val.v_num = qlink(&_qone_);
X                break;
X            case E_INC:
X                tmp.v_type = V_NUM;
X                tmp.v_num = &_qone_;
X                val = objcall(OBJ_ADD, v1, &tmp, NULL_VALUE);
X                break;
X            case E_DEC:
X                tmp.v_type = V_NUM;
X                tmp.v_num = &_qone_;
X                val = objcall(OBJ_SUB, v1, &tmp, NULL_VALUE);
X                break;
X            case E_SQUARE:
X                val = objcall(OBJ_MUL, v1, v1, NULL_VALUE);
X                break;
X            default:
X                math_error("Function \"%s\" is undefined", namefunc(index));
X        }
X        return val;
X    }
X    switch (oip->args) {
X        case 0:
X            break;
X        case 1:
X            ++stack;
X            stack->v_addr = v1;
X            stack->v_type = V_ADDR;
X            break;
X        case 2:
X            ++stack;
X            stack->v_addr = v1;
X            stack->v_type = V_ADDR;
X            ++stack;
X            stack->v_addr = v2;
X            stack->v_type = V_ADDR;
X            break;
X        case 3:
X            ++stack;
X            stack->v_addr = v1;
X            stack->v_type = V_ADDR;
X            ++stack;
X            stack->v_addr = v2;
X            stack->v_type = V_ADDR;
X            ++stack;
X            stack->v_addr = v3;
X            stack->v_type = V_ADDR;
X            break;
X        default:
X            math_error("Bad number of args to calculate");
X    }
X    calculate(fp, oip->args);
X    switch (oip->retval) {
X        case A_VALUE:
X            return *stack--;
X        case A_UNDEF:
X            freevalue(stack--);
X            val.v_type = V_NULL;
X            break;
X        case A_INT:
X            if ((stack->v_type != V_NUM) || qisfrac(stack->v_num))
X                math_error("Integer return value required");
X            index = qtoi(stack->v_num);
X            qfree(stack->v_num);
X            stack--;
X            val.v_type = V_INT;
X            val.v_int = index;
X            break;
X        default:
X            math_error("Bad object return");
X    }
X    return val;
X}
X
X
X/*
X * Routine called to clear the cache of known undefined functions for
X * the objects.  This changes negative indices back into positive ones
X * so that they will all be checked for existence again.
X */
Xvoid
Xobjuncache()
X{
X    register int *ip;
X    int i, j;
X
X    i = objectnames.h_count;
X    while (--i >= 0) {
X        ip = objects[i]->actions;
X        for (j = OBJ_MAXFUNC; j-- >= 0; ip++)
X            if (*ip < 0)
X                *ip = -*ip;
X    }
X}
X
X
X/*
X * Print the elements of an object in short and unambiguous format.
X * This is the default routine if the user's is not defined.
X */
Xstatic void
Xobjprint(op)
X    OBJECT *op;        /* object being printed */
X{
X    int count;        /* number of elements */
X    int i;            /* index */
X
X    count = op->o_actions->count;
X    math_fmt("obj %s {", op->o_actions->name);
X    for (i = 0; i < count; i++) {
X        if (i)
X            math_str(", ");
X        printvalue(&op->o_table[i], PRINT_SHORT | PRINT_UNAMBIG);
X    }
X    math_chr('}');
X}
X
X
X/*
X * Test an object for being "nonzero".
X * This is the default routine if the user's is not defined.
X * Returns TRUE if any of the elements are "nonzero".
X */
Xstatic BOOL
Xobjtest(op)
X    OBJECT *op;
X{
X    int i;            /* loop counter */
X
X    i = op->o_actions->count;
X    while (--i >= 0) {
X        if (testvalue(&op->o_table[i]))
X            return TRUE;
X    }
X    return FALSE;
X}
X
X
X/*
X * Compare two objects for equality, returning TRUE if they differ.
X * This is the default routine if the user's is not defined.
X * For equality, all elements must be equal.
X */
Xstatic BOOL
Xobjcmp(op1, op2)
X    OBJECT *op1, *op2;
X{
X    int i;            /* loop counter */
X
X    if (op1->o_actions != op2->o_actions)
X        return TRUE;
X    i = op1->o_actions->count;
X    while (--i >= 0) {
X        if (comparevalue(&op1->o_table[i], &op2->o_table[i]))
X            return TRUE;
X    }
X    return FALSE;
X}
X
X
X/*
X * Raise an object to an integral power.
X * This is the default routine if the user's is not defined.
X * Negative powers mean the positive power of the inverse.
X * Zero means the multiplicative identity.
X */
Xstatic VALUE
Xobjpowi(vp, q)
X    VALUE *vp;        /* value to be powered */
X    NUMBER *q;        /* power to raise number to */
X{
X    VALUE res, tmp;
X    long power;        /* power to raise to */
X    unsigned long bit;    /* current bit value */
X
X    if (qisfrac(q))
X        math_error("Raising object to non-integral power");
X    if (zisbig(q->num))
X        math_error("Raising object to very large power");
X    power = (zistiny(q->num) ? z1tol(q->num) : z2tol(q->num));
X    if (qisneg(q))
X        power = -power;
X    /*
X     * Handle some low powers specially
X     */
X    if ((power <= 2) && (power >= -2)) {
X        switch ((int) power) {
X            case 0:
X                return objcall(OBJ_ONE, vp, NULL_VALUE, NULL_VALUE);
X            case 1:
X                res.v_obj = objcopy(vp->v_obj);
X                res.v_type = V_OBJ;
X                return res;
X            case -1:
X                return objcall(OBJ_INV, vp, NULL_VALUE, NULL_VALUE);
X            case 2:
X                return objcall(OBJ_SQUARE, vp, NULL_VALUE, NULL_VALUE);
X        }
X    }
X    if (power < 0)
X        power = -power;
X    /*
X     * Compute the power by squaring and multiplying.
X     * This uses the left to right method of power raising.
X     */
X    bit = TOPFULL;
X    while ((bit & power) == 0)
X        bit >>= 1L;
X    bit >>= 1L;
X    res = objcall(OBJ_SQUARE, vp, NULL_VALUE, NULL_VALUE);
X    if (bit & power) {
X        tmp = objcall(OBJ_MUL, &res, vp, NULL_VALUE);
X        objfree(res.v_obj);
X        res = tmp;
X    }
X    bit >>= 1L;
X    while (bit) {
X        tmp = objcall(OBJ_SQUARE, &res, NULL_VALUE, NULL_VALUE);
X        objfree(res.v_obj);
X        res = tmp;
X        if (bit & power) {
X            tmp = objcall(OBJ_MUL, &res, vp, NULL_VALUE);
X            objfree(res.v_obj);
X            res = tmp;
X        }
X        bit >>= 1L;
X    }
X    if (qisneg(q)) {
X        tmp = objcall(OBJ_INV, &res, NULL_VALUE, NULL_VALUE);
X        objfree(res.v_obj);
X        return tmp;
X    }
X    return res;
X}
X
X
X/*
X * Define a (possibly) new class of objects.
X * The list of indexes for the element names is also specified here,
X * and the number of elements defined for the object.
X */
Xvoid
Xdefineobject(name, indices, count)
X    char *name;        /* name of object type */
X    int indices[];        /* table of indices for elements */
X{
X    OBJECTACTIONS *oap;    /* object definition structure */
X    STRINGHEAD *hp;
X    int index;
X
X    hp = &objectnames;
X    if (hp->h_list == NULL)
X        initstr(hp);
X    index = findstr(hp, name);
X    if (index >= 0) {
X        /*
X         * Object is already defined.  Give an error unless this
X         * new definition is exactly the same as the old one.
X         */
X        oap = objects[index];
X        if (oap->count == count) {
X            for (index = 0; ; index++) {
X                if (index >= count)
X                    return;
X                if (oap->elements[index] != indices[index])
X                    break;
X            }
X        }
X        math_error("Object type \"%s\" is already defined", name);
X    }
X
X    if (hp->h_count >= MAXOBJECTS)
X        math_error("Too many object types in use");
X    oap = (OBJECTACTIONS *) malloc(objectactionsize(count));
X    if (oap)
X        name = addstr(hp, name);
X    if ((oap == NULL) || (name == NULL))
X        math_error("Cannot allocate object type");
X    oap->name = name;
X    oap->count = count;
X    for (index = OBJ_MAXFUNC; index >= 0; index--)
X        oap->actions[index] = 0;
X    for (index = 0; index < count; index++)
X        oap->elements[index] = indices[index];
X    index = findstr(hp, name);
X    objects[index] = oap;
X    return;
X}
X
X
X/*
X * Check an object name to see if it is currently defined.
X * If so, the index for the object type is returned.
X * If the object name is currently unknown, then -1 is returned.
X */
Xint
Xcheckobject(name)
X    char *name;
X{
X    STRINGHEAD *hp;
X
X    hp = &objectnames;
X    if (hp->h_list == NULL)
X        return -1;
X    return findstr(hp, name);
X}
X
X
X/*
X * Define a (possibly) new element name for an object.
X * Returns an index which identifies the element name.
X */
Xint
Xaddelement(name)
X    char *name;
X{
X    STRINGHEAD *hp;
X    int index;
X
X    hp = &elements;
X    if (hp->h_list == NULL)
X        initstr(hp);
X    index = findstr(hp, name);
X    if (index >= 0)
X        return index;
X    if (addstr(hp, name) == NULL)
X        math_error("Cannot allocate element name");
X    return findstr(hp, name);
X}
X
X
X/*
X * Return the index which identifies an element name.
X * Returns minus one if the element name is unknown.
X */
Xint
Xfindelement(name)
X    char *name;        /* element name */
X{
X    if (elements.h_list == NULL)
X        return -1;
X    return findstr(&elements, name);
X}
X
X
X/*
X * Return the value table offset to be used for an object element name.
SHAR_EOF
echo "End of part 6"
echo "File calc2.9.0/obj.c is continued in part 7"
echo "7" > s2_seq_.tmp
exit 0