Newsgroups: comp.sources.unix
From: dbell@canb.auug.org.au (David I. Bell)
Subject: v27i133: calc-2.9.0 - arbitrary precision C-like programmable calculator, Part06/19
References: <1.755316719.21314@gw.home.vix.com>
Sender: unix-sources-moderator@gw.home.vix.com
Approved: vixie@gw.home.vix.com

Submitted-By: dbell@canb.auug.org.au (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