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

Text File | 1992-05-09 | 30.3 KB | 1,368 lines

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

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

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 14 (of 21)."
# Contents:  matfunc.c
# Wrapped by dbell@elm on Tue Feb 25 15:21:12 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'matfunc.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'matfunc.c'\"
else
echo shar: Extracting \"'matfunc.c'\" \(28075 characters\)
sed "s/^X//" >'matfunc.c' <<'END_OF_FILE'
X/*
X * Copyright (c) 1992 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Extended precision rational arithmetic matrix functions.
X * Matrices can contain arbitrary types of elements.
X */
X
X#include "calc.h"
X
X
X#if 0
Xstatic char *abortmsg = "Calculation aborted";
Xstatic char *memmsg = "Not enough memory";
X#endif
X
X
Xstatic void matswaprow(), matsubrow(), matmulrow();
X#if 0
Xstatic void mataddrow();
X#endif
X
Xstatic MATRIX *matident();
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        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            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        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            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        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        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        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        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        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        error("Raising non-square matrix to a power");
X    if (qisfrac(q))
X        error("Raising matrix to non-integral power");
X    if (isbig(q->num))
X        error("Raising matrix to very large power");
X    power = (istiny(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        error("Matrix not 1d for cross product");
X    if ((m1->m_size != 3) || (m2->m_size != 3))
X        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        error("Matrix not 1d for dot product");
X    if (m1->m_size != m2->m_size)
X        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        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        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
X/*
X * Transpose the elements of a two dimensional matrix.
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 old matrix */
X    long row, col;            /* current row and column */
X    MATRIX *res;
X
X    if (m->m_dim != 2)
X        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[0] - m->m_min[0] + 1);
X    cols = (m->m_max[1] - m->m_min[1] + 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        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        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 * 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            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        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        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        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        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                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        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        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#if 0
X/*
X * Local utility routine to add one row to another one.
X * No checks are made to verify the legality of the arguments.
X */
Xstatic void
Xmataddrow(m, r1, r2)
X    MATRIX *m;
X    long r1;    /* row to be added into */
X    long r2;    /* row to add */
X{
X    register VALUE *v1, *v2;
X    register long rows;
X    VALUE tmp;
X
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        addvalue(v1, v2, &tmp);
X        freevalue(v1);
X        *v1 = tmp;
X        v1++;
X        v2++;
X    }
X}
X#endif
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        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 * 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, maxprint)
X    MATRIX *m;
X    long maxprint;
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 = ((maxprint > 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 (maxprint > fullsize)
X        maxprint = 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 (maxprint <= 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 < maxprint; 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 (maxprint < fullsize)
X        math_str("  ...\n");
X}
X
X/* END CODE */
END_OF_FILE
if test 28075 -ne `wc -c <'matfunc.c'`; then
    echo shar: \"'matfunc.c'\" unpacked with wrong size!
fi
# end of 'matfunc.c'
fi
echo shar: End of archive 14 \(of 21\).
cp /dev/null ark14isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 21 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0