c_cpp/etc/calc/matfunc.c

   1 /*
   2  * matfunc - extended precision rational arithmetic matrix functions
   3  *
   4  * Copyright (C) 1999-2007  David I. Bell
   5  *
   6  * Calc is open software; you can redistribute it and/or modify it under
   7  * the terms of the version 2.1 of the GNU Lesser General Public License
   8  * as published by the Free Software Foundation.
   9  *
  10  * Calc is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13  * Public License for more details.
  14  *
  15  * A copy of version 2.1 of the GNU Lesser General Public License is
  16  * distributed with calc under the filename COPYING-LGPL.  You should have
  17  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  18  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  *
  20  * @(#) $Revision: 30.2 $
  21  * @(#) $Id: matfunc.c,v 30.2 2007/07/05 13:30:38 chongo Exp $
  22  * @(#) $Source: /usr/local/src/bin/calc/RCS/matfunc.c,v $
  23  *
  24  * Under source code control:   1990/02/15 01:48:18
  25  * File existed as early as:    before 1990
  26  *
  27  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  28  */
  29
  30 /*
  31  * Extended precision rational arithmetic matrix functions.
  32  * Matrices can contain arbitrary types of elements.
  33  */
  34
  35 #include "value.h"
  36 #include "zrand.h"
  37 #include "calcerr.h"
  38
  39 #include "have_unused.h"
  40
  41 E_FUNC long irand(long s);
  42
  43 S_FUNC void matswaprow(MATRIX *m, long r1, long r2);
  44 S_FUNC void matsubrow(MATRIX *m, long oprow, long baserow, VALUE *mulval);
  45 S_FUNC void matmulrow(MATRIX *m, long row, VALUE *mulval);
  46 S_FUNC MATRIX *matident(MATRIX *m);
  47
  48
  49
  50 /*
  51  * Add two compatible matrices.
  52  */
  53 MATRIX *
  54 matadd(MATRIX *m1, MATRIX *m2)
  55 {
  56         int dim;
  57
  58         long min1, min2, max1, max2, index;
  59         VALUE *v1, *v2, *vres;
  60         MATRIX *res;
  61         MATRIX tmp;
  62
  63         if (m1->m_dim != m2->m_dim) {
  64                 math_error("Incompatible matrix dimensions for add");
  65                 /*NOTREACHED*/
  66         }
  67         tmp.m_dim = m1->m_dim;
  68         tmp.m_size = m1->m_size;
  69         for (dim = 0; dim < m1->m_dim; dim++) {
  70                 min1 = m1->m_min[dim];
  71                 max1 = m1->m_max[dim];
  72                 min2 = m2->m_min[dim];
  73                 max2 = m2->m_max[dim];
  74                 if ((min1 && min2 && (min1 != min2)) ||
  75                     ((max1-min1) != (max2-min2))) {
  76                         math_error("Incompatible matrix bounds for add");
  77                         /*NOTREACHED*/
  78                 }
  79                 tmp.m_min[dim] = (min1 ? min1 : min2);
  80                 tmp.m_max[dim] = tmp.m_min[dim] + (max1 - min1);
  81         }
  82         res = matalloc(m1->m_size);
  83         *res = tmp;
  84         v1 = m1->m_table;
  85         v2 = m2->m_table;
  86         vres = res->m_table;
  87         for (index = m1->m_size; index > 0; index--)
  88                 addvalue(v1++, v2++, vres++);
  89         return res;
  90 }
  91
  92
  93 /*
  94  * Subtract two compatible matrices.
  95  */
  96 MATRIX *
  97 matsub(MATRIX *m1, MATRIX *m2)
  98 {
  99         int dim;
 100         long min1, min2, max1, max2, index;
 101         VALUE *v1, *v2, *vres;
 102         MATRIX *res;
 103         MATRIX tmp;
 104
 105         if (m1->m_dim != m2->m_dim) {
 106                 math_error("Incompatible matrix dimensions for sub");
 107                 /*NOTREACHED*/
 108         }
 109         tmp.m_dim = m1->m_dim;
 110         tmp.m_size = m1->m_size;
 111         for (dim = 0; dim < m1->m_dim; dim++) {
 112                 min1 = m1->m_min[dim];
 113                 max1 = m1->m_max[dim];
 114                 min2 = m2->m_min[dim];
 115                 max2 = m2->m_max[dim];
 116                 if ((min1 && min2 && (min1 != min2)) ||
 117                     ((max1-min1) != (max2-min2))) {
 118                         math_error("Incompatible matrix bounds for sub");
 119                         /*NOTREACHED*/
 120                 }
 121                 tmp.m_min[dim] = (min1 ? min1 : min2);
 122                 tmp.m_max[dim] = tmp.m_min[dim] + (max1 - min1);
 123         }
 124         res = matalloc(m1->m_size);
 125         *res = tmp;
 126         v1 = m1->m_table;
 127         v2 = m2->m_table;
 128         vres = res->m_table;
 129         for (index = m1->m_size; index > 0; index--)
 130                 subvalue(v1++, v2++, vres++);
 131         return res;
 132 }
 133
 134
 135 /*
 136  * Produce the negative of a matrix.
 137  */
 138 MATRIX *
 139 matneg(MATRIX *m)
 140 {
 141         register VALUE *val, *vres;
 142         long index;
 143         MATRIX *res;
 144
 145         res = matalloc(m->m_size);
 146         *res = *m;
 147         val = m->m_table;
 148         vres = res->m_table;
 149         for (index = m->m_size; index > 0; index--)
 150                 negvalue(val++, vres++);
 151         return res;
 152 }
 153
 154
 155 /*
 156  * Multiply two compatible matrices.
 157  */
 158 MATRIX *
 159 matmul(MATRIX *m1, MATRIX *m2)
 160 {
 161         register MATRIX *res;
 162         long i1, i2, max1, max2, index, maxindex;
 163         VALUE *v1, *v2, *vres;
 164         VALUE sum, tmp1, tmp2;
 165
 166         if (m1->m_dim == 0) {
 167                 i2 = m2->m_size;
 168                 v2 = m2->m_table;
 169                 res = matalloc(i2);
 170                 *res = *m2;
 171                 vres = res->m_table;
 172                 while (i2-- > 0)
 173                         mulvalue(m1->m_table, v2++, vres++);
 174                 return res;
 175         }
 176         if (m2->m_dim == 0) {
 177                 i1 = m1->m_size;
 178                 v1 = m1->m_table;
 179                 res = matalloc(i1);
 180                 *res = *m1;
 181                 vres = res->m_table;
 182                 while (i1-- > 0)
 183                         mulvalue(v1++, m2->m_table, vres++);
 184                 return res;
 185         }
 186         if (m1->m_dim == 1 && m2->m_dim == 1) {
 187                 if (m1->m_max[0]-m1->m_min[0] != m2->m_max[0]-m2->m_min[0]) {
 188                         math_error("Incompatible bounds for 1D * 1D  matmul");
 189                         /*NOTREACHED*/
 190                 }
 191                 res = matalloc(m1->m_size);
 192                 *res = *m1;
 193                 v1 = m1->m_table;
 194                 v2 = m2->m_table;
 195                 vres = res->m_table;
 196                 for (index = m1->m_size; index > 0; index--)
 197                         mulvalue(v1++, v2++, vres++);
 198                 return res;
 199         }
 200         if (m1->m_dim == 1 && m2->m_dim == 2) {
 201                 if (m1->m_max[0]-m1->m_min[0] != m2->m_max[0]-m2->m_min[0]) {
 202                         math_error("Incompatible bounds for 1D * 2D matmul");
 203                         /*NOTREACHED*/
 204                 }
 205                 res = matalloc(m2->m_size);
 206                 *res = *m2;
 207                 i1 = m1->m_max[0] - m1->m_min[0] + 1;
 208                 max2 = m2->m_max[1] - m2->m_min[1] + 1;
 209                 v1 = m1->m_table;
 210                 v2 = m2->m_table;
 211                 vres = res->m_table;
 212                 while (i1-- > 0) {
 213                         i2 = max2;
 214                         while (i2-- > 0)
 215                                 mulvalue(v1, v2++, vres++);
 216                         v1++;
 217                 }
 218                 return res;
 219         }
 220         if (m1->m_dim == 2 && m2->m_dim == 1) {
 221                 if (m1->m_max[1]-m1->m_min[1] != m2->m_max[0]-m2->m_min[0]) {
 222                         math_error("Incompatible bounds for 2D * 1D matmul");
 223                         /*NOTREACHED*/
 224                 }
 225                 res = matalloc(m1->m_size);
 226                 *res = *m1;
 227                 i1 = m1->m_max[0] - m1->m_min[0] + 1;
 228                 max1 = m1->m_max[1] - m1->m_min[1] + 1;
 229                 v1 = m1->m_table;
 230                 vres = res->m_table;
 231                 while (i1-- > 0) {
 232                         v2 = m2->m_table;
 233                         i2 = max1;
 234                         while (i2-- > 0)
 235                                 mulvalue(v1++, v2++, vres++);
 236                 }
 237                 return res;
 238         }
 239
 240         if ((m1->m_dim != 2) || (m2->m_dim != 2)) {
 241                 math_error("Matrix dimensions not compatible for mul");
 242                 /*NOTREACHED*/
 243         }
 244         if ((m1->m_max[1]-m1->m_min[1]) != (m2->m_max[0]-m2->m_min[0])) {
 245                 math_error("Incompatible bounds for 2D * 2D matrix mul");
 246                 /*NOTREACHED*/
 247         }
 248         max1 = (m1->m_max[0] - m1->m_min[0] + 1);
 249         max2 = (m2->m_max[1] - m2->m_min[1] + 1);
 250         maxindex = (m1->m_max[1] - m1->m_min[1] + 1);
 251         res = matalloc(max1 * max2);
 252         res->m_dim = 2;
 253         res->m_min[0] = m1->m_min[0];
 254         res->m_max[0] = m1->m_max[0];
 255         res->m_min[1] = m2->m_min[1];
 256         res->m_max[1] = m2->m_max[1];
 257         for (i1 = 0; i1 < max1; i1++) {
 258                 for (i2 = 0; i2 < max2; i2++) {
 259                         sum.v_type = V_NULL;
 260                         sum.v_subtype = V_NOSUBTYPE;
 261                         v1 = &m1->m_table[i1 * maxindex];
 262                         v2 = &m2->m_table[i2];
 263                         for (index = 0; index < maxindex; index++) {
 264                                 mulvalue(v1, v2, &tmp1);
 265                                 addvalue(&sum, &tmp1, &tmp2);
 266                                 freevalue(&tmp1);
 267                                 freevalue(&sum);
 268                                 sum = tmp2;
 269                                 v1++;
 270                                 if (index+1 < maxindex)
 271                                         v2 += max2;
 272                         }
 273                         index = (i1 * max2) + i2;
 274                         res->m_table[index] = sum;
 275                 }
 276         }
 277         return res;
 278 }
 279
 280
 281 /*
 282  * Square a matrix.
 283  */
 284 MATRIX *
 285 matsquare(MATRIX *m)
 286 {
 287         register MATRIX *res;
 288         long i1, i2, max, index;
 289         VALUE *v1, *v2;
 290         VALUE sum, tmp1, tmp2;
 291
 292         if (m->m_dim < 2) {
 293                 res = matalloc(m->m_size);
 294                 *res = *m;
 295                 v1 = m->m_table;
 296                 v2 = res->m_table;
 297                 for (index = m->m_size; index > 0; index--)
 298                         squarevalue(v1++, v2++);
 299                 return res;
 300         }
 301         if (m->m_dim != 2) {
 302                 math_error("Matrix dimension exceeds two for square");
 303                 /*NOTREACHED*/
 304         }
 305         if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) {
 306                 math_error("Squaring non-square matrix");
 307                 /*NOTREACHED*/
 308         }
 309         max = (m->m_max[0] - m->m_min[0] + 1);
 310         res = matalloc(max * max);
 311         res->m_dim = 2;
 312         res->m_min[0] = m->m_min[0];
 313         res->m_max[0] = m->m_max[0];
 314         res->m_min[1] = m->m_min[1];
 315         res->m_max[1] = m->m_max[1];
 316         for (i1 = 0; i1 < max; i1++) {
 317                 for (i2 = 0; i2 < max; i2++) {
 318                         sum.v_type = V_NULL;
 319                         sum.v_subtype = V_NOSUBTYPE;
 320                         v1 = &m->m_table[i1 * max];
 321                         v2 = &m->m_table[i2];
 322                         for (index = 0; index < max; index++) {
 323                                 mulvalue(v1, v2, &tmp1);
 324                                 addvalue(&sum, &tmp1, &tmp2);
 325                                 freevalue(&tmp1);
 326                                 freevalue(&sum);
 327                                 sum = tmp2;
 328                                 v1++;
 329                                 v2 += max;
 330                         }
 331                         index = (i1 * max) + i2;
 332                         res->m_table[index] = sum;
 333                 }
 334         }
 335         return res;
 336 }
 337
 338
 339 /*
 340  * Compute the result of raising a matrix to an integer power if
 341  * dimension <= 2 and for dimension == 2, the matrix is square.
 342  * Negative powers mean the positive power of the inverse.
 343  * Note: This calculation could someday be improved for large powers
 344  * by using the characteristic polynomial of the matrix.
 345  *
 346  * given:
 347  *      m               matrix to be raised
 348  *      q               power to raise it to
 349  */
 350 MATRIX *
 351 matpowi(MATRIX *m, NUMBER *q)
 352 {
 353         MATRIX *res, *tmp;
 354         long power;             /* power to raise to */
 355         FULL bit;               /* current bit value */
 356
 357         if (m->m_dim > 2) {
 358                 math_error("Matrix dimension greater than 2 for power");
 359                 /*NOTREACHED*/
 360         }
 361         if (m->m_dim == 2 && (m->m_max[0] - m->m_min[0] !=
 362                          m->m_max[1] - m->m_min[1])) {
 363                 math_error("Raising non-square 2D matrix to a power");
 364                 /*NOTREACHED*/
 365         }
 366         if (qisfrac(q)) {
 367                 math_error("Raising matrix to non-integral power");
 368                 /*NOTREACHED*/
 369         }
 370         if (zge31b(q->num)) {
 371                 math_error("Raising matrix to very large power");
 372                 /*NOTREACHED*/
 373         }
 374         power = ztolong(q->num);
 375         if (qisneg(q))
 376                 power = -power;
 377         /*
 378          * Handle some low powers specially
 379          */
 380         if ((power <= 4) && (power >= -2)) {
 381                 switch ((int) power) {
 382                 case 0:
 383                         return matident(m);
 384                 case 1:
 385                         return matcopy(m);
 386                 case -1:
 387                         return matinv(m);
 388                 case 2:
 389                         return matsquare(m);
 390                 case -2:
 391                         tmp = matinv(m);
 392                         res = matsquare(tmp);
 393                         matfree(tmp);
 394                         return res;
 395                 case 3:
 396                         tmp = matsquare(m);
 397                         res = matmul(m, tmp);
 398                         matfree(tmp);
 399                         return res;
 400                 case 4:
 401                         tmp = matsquare(m);
 402                         res = matsquare(tmp);
 403                         matfree(tmp);
 404                         return res;
 405                 }
 406         }
 407         if (power < 0) {
 408                 m = matinv(m);
 409                 power = -power;
 410         }
 411         /*
 412          * Compute the power by squaring and multiplying.
 413          * This uses the left to right method of power raising.
 414          */
 415         bit = TOPFULL;
 416         while ((bit & power) == 0)
 417                 bit >>= 1L;
 418         bit >>= 1L;
 419         res = matsquare(m);
 420         if (bit & power) {
 421                 tmp = matmul(res, m);
 422                 matfree(res);
 423                 res = tmp;
 424         }
 425         bit >>= 1L;
 426         while (bit) {
 427                 tmp = matsquare(res);
 428                 matfree(res);
 429                 res = tmp;
 430                 if (bit & power) {
 431                         tmp = matmul(res, m);
 432                         matfree(res);
 433                         res = tmp;
 434                 }
 435                 bit >>= 1L;
 436         }
 437         if (qisneg(q))
 438                 matfree(m);
 439         return res;
 440 }
 441
 442
 443 /*
 444  * Calculate the cross product of two one dimensional matrices each
 445  * with three components.
 446  *      m3 = matcross(m1, m2);
 447  */
 448 MATRIX *
 449 matcross(MATRIX *m1, MATRIX *m2)
 450 {
 451         MATRIX *res;
 452         VALUE *v1, *v2, *vr;
 453         VALUE tmp1, tmp2;
 454
 455         res = matalloc(3L);
 456         res->m_dim = 1;
 457         res->m_min[0] = 0;
 458         res->m_max[0] = 2;
 459         v1 = m1->m_table;
 460         v2 = m2->m_table;
 461         vr = res->m_table;
 462         mulvalue(v1 + 1, v2 + 2, &tmp1);
 463         mulvalue(v1 + 2, v2 + 1, &tmp2);
 464         subvalue(&tmp1, &tmp2, vr + 0);
 465         freevalue(&tmp1);
 466         freevalue(&tmp2);
 467         mulvalue(v1 + 2, v2 + 0, &tmp1);
 468         mulvalue(v1 + 0, v2 + 2, &tmp2);
 469         subvalue(&tmp1, &tmp2, vr + 1);
 470         freevalue(&tmp1);
 471         freevalue(&tmp2);
 472         mulvalue(v1 + 0, v2 + 1, &tmp1);
 473         mulvalue(v1 + 1, v2 + 0, &tmp2);
 474         subvalue(&tmp1, &tmp2, vr + 2);
 475         freevalue(&tmp1);
 476         freevalue(&tmp2);
 477         return res;
 478 }
 479
 480
 481 /*
 482  * Return the dot product of two matrices.
 483  *      result = matdot(m1, m2);
 484  */
 485 VALUE
 486 matdot(MATRIX *m1, MATRIX *m2)
 487 {
 488         VALUE *v1, *v2;
 489         VALUE result, tmp1, tmp2;
 490         long len;
 491
 492         v1 = m1->m_table;
 493         v2 = m2->m_table;
 494         mulvalue(v1, v2, &result);
 495         len = m1->m_size;
 496         while (--len > 0) {
 497                 mulvalue(++v1, ++v2, &tmp1);
 498                 addvalue(&result, &tmp1, &tmp2);
 499                 freevalue(&tmp1);
 500                 freevalue(&result);
 501                 result = tmp2;
 502         }
 503         return result;
 504 }
 505
 506
 507 /*
 508  * Scale the elements of a matrix by a specified power of two.
 509  *
 510  * given:
 511  *      m               matrix to be scaled
 512  *      n               scale factor
 513  */
 514 MATRIX *
 515 matscale(MATRIX *m, long n)
 516 {
 517         register VALUE *val, *vres;
 518         VALUE temp;
 519         long index;
 520         MATRIX *res;            /* resulting matrix */
 521
 522         if (n == 0)
 523                 return matcopy(m);
 524         temp.v_type = V_NUM;
 525         temp.v_subtype = V_NOSUBTYPE;
 526         temp.v_num = itoq(n);
 527         res = matalloc(m->m_size);
 528         *res = *m;
 529         val = m->m_table;
 530         vres = res->m_table;
 531         for (index = m->m_size; index > 0; index--)
 532                 scalevalue(val++, &temp, vres++);
 533         qfree(temp.v_num);
 534         return res;
 535 }
 536
 537
 538 /*
 539  * Shift the elements of a matrix by the specified number of bits.
 540  * Positive shift means leftwards, negative shift rightwards.
 541  *
 542  * given:
 543  *      m               matrix to be shifted
 544  *      n               shift count
 545  */
 546 MATRIX *
 547 matshift(MATRIX *m, long n)
 548 {
 549         register VALUE *val, *vres;
 550         VALUE temp;
 551         long index;
 552         MATRIX *res;            /* resulting matrix */
 553
 554         if (n == 0)
 555                 return matcopy(m);
 556         temp.v_type = V_NUM;
 557         temp.v_subtype = V_NOSUBTYPE;
 558         temp.v_num = itoq(n);
 559         res = matalloc(m->m_size);
 560         *res = *m;
 561         val = m->m_table;
 562         vres = res->m_table;
 563         for (index = m->m_size; index > 0; index--)
 564                 shiftvalue(val++, &temp, FALSE, vres++);
 565         qfree(temp.v_num);
 566         return res;
 567 }
 568
 569
 570 /*
 571  * Multiply the elements of a matrix by a specified value.
 572  *
 573  * given:
 574  *      m               matrix to be multiplied
 575  *      vp              value to multiply by
 576  */
 577 MATRIX *
 578 matmulval(MATRIX *m, VALUE *vp)
 579 {
 580         register VALUE *val, *vres;
 581         long index;
 582         MATRIX *res;
 583
 584         res = matalloc(m->m_size);
 585         *res = *m;
 586         val = m->m_table;
 587         vres = res->m_table;
 588         for (index = m->m_size; index > 0; index--)
 589                 mulvalue(val++, vp, vres++);
 590         return res;
 591 }
 592
 593
 594 /*
 595  * Divide the elements of a matrix by a specified value, keeping
 596  * only the integer quotient.
 597  *
 598  * given:
 599  *      m               matrix to be divided
 600  *      vp              value to divide by
 601  *      v3              rounding type parameter
 602  */
 603 MATRIX *
 604 matquoval(MATRIX *m, VALUE *vp, VALUE *v3)
 605 {
 606         register VALUE *val, *vres;
 607         long index;
 608         MATRIX *res;
 609
 610         if ((vp->v_type == V_NUM) && qiszero(vp->v_num)) {
 611                 math_error("Division by zero");
 612                 /*NOTREACHED*/
 613         }
 614         res = matalloc(m->m_size);
 615         *res = *m;
 616         val = m->m_table;
 617         vres = res->m_table;
 618         for (index = m->m_size; index > 0; index--)
 619                 quovalue(val++, vp, v3, vres++);
 620         return res;
 621 }
 622
 623
 624 /*
 625  * Divide the elements of a matrix by a specified value, keeping
 626  * only the remainder of the division.
 627  *
 628  * given:
 629  *      m               matrix to be divided
 630  *      vp              value to divide by
 631  *      v3              rounding type parameter
 632  */
 633 MATRIX *
 634 matmodval(MATRIX *m, VALUE *vp, VALUE *v3)
 635 {
 636         register VALUE *val, *vres;
 637         long index;
 638         MATRIX *res;
 639
 640         if ((vp->v_type == V_NUM) && qiszero(vp->v_num)) {
 641                 math_error("Division by zero");
 642                 /*NOTREACHED*/
 643         }
 644         res = matalloc(m->m_size);
 645         *res = *m;
 646         val = m->m_table;
 647         vres = res->m_table;
 648         for (index = m->m_size; index > 0; index--)
 649                 modvalue(val++, vp, v3, vres++);
 650         return res;
 651 }
 652
 653
 654 VALUE
 655 mattrace(MATRIX *m)
 656 {
 657         VALUE *vp;
 658         VALUE sum;
 659         VALUE tmp;
 660         long i, j;
 661
 662         if (m->m_dim < 2) {
 663                 matsum(m, &sum);
 664                 return sum;
 665         }
 666         if (m->m_dim != 2)
 667                 return error_value(E_MATTRACE2);
 668         i = (m->m_max[0] - m->m_min[0] + 1);
 669         j = (m->m_max[1] - m->m_min[1] + 1);
 670         if (i != j)
 671                 return error_value(E_MATTRACE3);
 672         vp = m->m_table;
 673         copyvalue(vp, &sum);
 674         j++;
 675         while (--i > 0) {
 676                 vp += j;
 677                 addvalue(&sum, vp, &tmp);
 678                 freevalue(&sum);
 679                 sum = tmp;
 680         }
 681         return sum;
 682 }
 683
 684
 685 /*
 686  * Transpose a 2-dimensional matrix
 687  */
 688 MATRIX *
 689 mattrans(MATRIX *m)
 690 {
 691         register VALUE *v1, *v2;        /* current values */
 692         long rows, cols;                /* rows and columns in new matrix */
 693         long row, col;                  /* current row and column */
 694         MATRIX *res;
 695
 696         if (m->m_dim < 2)
 697                 return matcopy(m);
 698         res = matalloc(m->m_size);
 699         res->m_dim = 2;
 700         res->m_min[0] = m->m_min[1];
 701         res->m_max[0] = m->m_max[1];
 702         res->m_min[1] = m->m_min[0];
 703         res->m_max[1] = m->m_max[0];
 704         rows = (m->m_max[1] - m->m_min[1] + 1);
 705         cols = (m->m_max[0] - m->m_min[0] + 1);
 706         v1 = res->m_table;
 707         for (row = 0; row < rows; row++) {
 708                 v2 = &m->m_table[row];
 709                 for (col = 0; col < cols; col++) {
 710                         copyvalue(v2, v1);
 711                         v1++;
 712                         v2 += rows;
 713                 }
 714         }
 715         return res;
 716 }
 717
 718
 719 /*
 720  * Produce a matrix with values all of which are conjugated.
 721  */
 722 MATRIX *
 723 matconj(MATRIX *m)
 724 {
 725         register VALUE *val, *vres;
 726         long index;
 727         MATRIX *res;
 728
 729         res = matalloc(m->m_size);
 730         *res = *m;
 731         val = m->m_table;
 732         vres = res->m_table;
 733         for (index = m->m_size; index > 0; index--)
 734                 conjvalue(val++, vres++);
 735         return res;
 736 }
 737
 738
 739 /*
 740  * Round elements of a matrix to specified number of decimal digits
 741  */
 742 MATRIX *
 743 matround(MATRIX *m, VALUE *v2, VALUE *v3)
 744 {
 745         VALUE *p, *q;
 746         long s;
 747         MATRIX *res;
 748
 749         s = m->m_size;
 750         res = matalloc(s);
 751         *res = *m;
 752         p = m->m_table;
 753         q = res->m_table;
 754         while (s-- > 0)
 755                 roundvalue(p++, v2, v3, q++);
 756         return res;
 757 }
 758
 759
 760 /*
 761  * Round elements of a matrix to specified number of binary digits
 762  */
 763 MATRIX *
 764 matbround(MATRIX *m, VALUE *v2, VALUE *v3)
 765 {
 766         VALUE *p, *q;
 767         long s;
 768         MATRIX *res;
 769
 770         s = m->m_size;
 771         res = matalloc(s);
 772         *res = *m;
 773         p = m->m_table;
 774         q = res->m_table;
 775         while (s-- > 0)
 776                 broundvalue(p++, v2, v3, q++);
 777         return res;
 778 }
 779
 780 /*
 781  * Approximate a matrix by approximating elemenbs to be multiples of
 782  * v2, rounding type determined by v3.
 783  */
 784 MATRIX *
 785 matappr(MATRIX *m, VALUE *v2, VALUE *v3)
 786 {
 787         VALUE *p, *q;
 788         long s;
 789         MATRIX *res;
 790
 791         s = m->m_size;
 792         res = matalloc(s);
 793         *res = *m;
 794         p = m->m_table;
 795         q = res->m_table;
 796         while (s-- > 0)
 797                 apprvalue(p++, v2, v3, q++);
 798         return res;
 799 }
 800
 801
 802
 803
 804 /*
 805  * Produce a matrix with values all of which have been truncated to integers.
 806  */
 807 MATRIX *
 808 matint(MATRIX *m)
 809 {
 810         register VALUE *val, *vres;
 811         long index;
 812         MATRIX *res;
 813
 814         res = matalloc(m->m_size);
 815         *res = *m;
 816         val = m->m_table;
 817         vres = res->m_table;
 818         for (index = m->m_size; index > 0; index--)
 819                 intvalue(val++, vres++);
 820         return res;
 821 }
 822
 823
 824 /*
 825  * Produce a matrix with values all of which have only the fraction part left.
 826  */
 827 MATRIX *
 828 matfrac(MATRIX *m)
 829 {
 830         register VALUE *val, *vres;
 831         long index;
 832         MATRIX *res;
 833
 834         res = matalloc(m->m_size);
 835         *res = *m;
 836         val = m->m_table;
 837         vres = res->m_table;
 838         for (index = m->m_size; index > 0; index--)
 839                 fracvalue(val++, vres++);
 840         return res;
 841 }
 842
 843
 844 /*
 845  * Index a matrix normally by the specified set of index values.
 846  * Returns the address of the matrix element if it is valid, or generates
 847  * an error if the index values are out of range.  The create flag is TRUE
 848  * if the element is to be written, but this is ignored here.
 849  *
 850  * given:
 851  *      mp              matrix to operate on
 852  *      create          TRUE => create if element does not exist
 853  *      dim             dimension of the indexing
 854  *      indices         table of values being indexed by
 855  */
 856 /*ARGSUSED*/
 857 VALUE *
 858 matindex(MATRIX *mp, BOOL UNUSED create, long dim, VALUE *indices)
 859 {
 860         NUMBER *q;              /* index value */
 861         VALUE *vp;
 862         long index;             /* index value as an integer */
 863         long offset;            /* current offset into array */
 864         int i;                  /* loop counter */
 865
 866         if (dim < 0) {
 867                 math_error("Negative dimension %ld for matrix", dim);
 868                 /*NOTREACHED*/
 869         }
 870         for (;;) {
 871                 if (dim <  mp->m_dim) {
 872                         math_error(
 873                             "Indexing a %ldd matrix as a %ldd matrix",
 874                             mp->m_dim, dim);
 875                         /*NOTREACHED*/
 876                 }
 877                 offset = 0;
 878                 for (i = 0; i < mp->m_dim; i++) {
 879                         if (indices->v_type != V_NUM) {
 880                                 math_error("Non-numeric index for matrix");
 881                                 /*NOTREACHED*/
 882                         }
 883                         q = indices->v_num;
 884                         if (qisfrac(q)) {
 885                                 math_error("Non-integral index for matrix");
 886                                 /*NOTREACHED*/
 887                         }
 888                         index = qtoi(q);
 889                         if (zge31b(q->num) || (index < mp->m_min[i]) ||
 890                             (index > mp->m_max[i])) {
 891                                 math_error("Index out of bounds for matrix");
 892                                 /*NOTREACHED*/
 893                         }
 894                         offset *= (mp->m_max[i] - mp->m_min[i] + 1);
 895                         offset += (index - mp->m_min[i]);
 896                         indices++;
 897                 }
 898                 vp = mp->m_table + offset;
 899                 dim -= mp->m_dim;
 900                 if (dim == 0)
 901                         break;
 902                 if (vp->v_type != V_MAT) {
 903                         math_error("Non-matrix argument for matindex");
 904                         /*NOTREACHED*/
 905                 }
 906                 mp = vp->v_mat;
 907         }
 908         return vp;
 909 }
 910
 911
 912 /*
 913  * Returns the list of indices for a matrix element with specified
 914  * double-bracket index.
 915  */
 916 LIST *
 917 matindices(MATRIX *mp, long index)
 918 {
 919         LIST *lp;
 920         int j;
 921         long d;
 922         VALUE val;
 923
 924         if (index < 0 || index >= mp->m_size)
 925                 return NULL;
 926
 927         lp = listalloc();
 928         val.v_type = V_NUM;
 929         j = mp->m_dim;
 930
 931         while (--j >= 0) {
 932                 d = mp->m_max[j] - mp->m_min[j] + 1;
 933                 val.v_num = itoq(index % d + mp->m_min[j]);
 934                 insertlistfirst(lp, &val);
 935                 qfree(val.v_num);
 936                 index /= d;
 937         }
 938         return lp;
 939 }
 940
 941
 942 /*
 943  * Search a matrix for the specified value, starting with the specified index.
 944  * Returns 0 and stores index if value found; otherwise returns 1.
 945  */
 946 int
 947 matsearch(MATRIX *m, VALUE *vp, long i, long j, ZVALUE *index)
 948 {
 949         register VALUE *val;
 950
 951         val = &m->m_table[i];
 952         if (i < 0 || j > m->m_size) {
 953                 math_error("This should not happen in call to matsearch");
 954                 /*NOTREACHED*/
 955         }
 956         while (i < j) {
 957                 if (acceptvalue(val++, vp)) {
 958                         utoz(i, index);
 959                         return 0;
 960                 }
 961                 i++;
 962         }
 963         return 1;
 964 }
 965
 966
 967 /*
 968  * Search a matrix backwards for the specified value, starting with the
 969  * specified index.  Returns 0 and stores index if value found; otherwise
 970  * returns 1.
 971  */
 972 int
 973 matrsearch(MATRIX *m, VALUE *vp, long i, long j, ZVALUE *index)
 974 {
 975         register VALUE *val;
 976
 977         if (i < 0 || j > m->m_size) {
 978                 math_error("This should not happen in call to matrsearch");
 979                 /*NOTREACHED*/
 980         }
 981         val = &m->m_table[--j];
 982         while (j >= i) {
 983                 if (acceptvalue(val--, vp)) {
 984                         utoz(j, index);
 985                         return 0;
 986                 }
 987                 j--;
 988         }
 989         return 1;
 990 }
 991
 992 /*
 993  * Fill all of the elements of a matrix with one of two specified values.
 994  * All entries are filled with the first specified value, except that if
 995  * the matrix is w-dimensional and the second value pointer is non-NULL, then
 996  * all diagonal entries are filled with the second value.  This routine
 997  * affects the supplied matrix directly, and doesn't return a copy.
 998  *
 999  * given:
1000  *      m               matrix to be filled
1001  *      v1              value to fill most of matrix with
1002  *      v2              value for diagonal entries or null
1003  */
1004 void
1005 matfill(MATRIX *m, VALUE *v1, VALUE *v2)
1006 {
1007         register VALUE *val;
1008         VALUE temp1, temp2;
1009         long rows, cols;
1010         long i, j;
1011
1012         copyvalue(v1, &temp1);
1013
1014         val = m->m_table;
1015         if (m->m_dim != 2 || v2 == NULL) {
1016                 for (i = m->m_size; i > 0; i--) {
1017                         freevalue(val);
1018                         copyvalue(&temp1, val++);
1019                 }
1020                 freevalue(&temp1);
1021                 return;
1022         }
1023
1024         copyvalue(v2, &temp2);
1025         rows = m->m_max[0] - m->m_min[0] + 1;
1026         cols = m->m_max[1] - m->m_min[1] + 1;
1027
1028         for (i = 0; i < rows; i++) {
1029                 for (j = 0; j < cols; j++) {
1030                         freevalue(val);
1031                         if (i == j)
1032                                 copyvalue(&temp2, val++);
1033                         else
1034                                 copyvalue(&temp1, val++);
1035                 }
1036         }
1037         freevalue(&temp1);
1038         freevalue(&temp2);
1039 }
1040
1041
1042
1043 /*
1044  * Set a copy of a square matrix to the identity matrix.
1045  */
1046 S_FUNC MATRIX *
1047 matident(MATRIX *m)
1048 {
1049         register VALUE *val;    /* current value */
1050         long row, col;          /* current row and column */
1051         long rows;              /* number of rows */
1052         MATRIX *res;            /* resulting matrix */
1053
1054         if (m->m_dim != 2) {
1055                 math_error(
1056                     "Matrix dimension must be two for setting to identity");
1057                 /*NOTREACHED*/
1058         }
1059         if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) {
1060                 math_error("Matrix must be square for setting to identity");
1061                 /*NOTREACHED*/
1062         }
1063         res = matalloc(m->m_size);
1064         *res = *m;
1065         val = res->m_table;
1066         rows = (res->m_max[0] - res->m_min[0] + 1);
1067         for (row = 0; row < rows; row++) {
1068                 for (col = 0; col < rows; col++) {
1069                         val->v_type = V_NUM;
1070                         val->v_num = ((row == col) ? qlink(&_qone_) :
1071                                                      qlink(&_qzero_));
1072                         val++;
1073                 }
1074         }
1075         return res;
1076 }
1077
1078
1079 /*
1080  * Calculate the inverse of a matrix if it exists.
1081  * This is done by using transformations on the supplied matrix to convert
1082  * it to the identity matrix, and simultaneously applying the same set of
1083  * transformations to the identity matrix.
1084  */
1085 MATRIX *
1086 matinv(MATRIX *m)
1087 {
1088         MATRIX *res;            /* matrix to become the inverse */
1089         long rows;              /* number of rows */
1090         long cur;               /* current row being worked on */
1091         long row, col;          /* temp row and column values */
1092         VALUE *val;             /* current value in matrix*/
1093         VALUE *vres;            /* current value in result for dim < 2 */
1094         VALUE mulval;           /* value to multiply rows by */
1095         VALUE tmpval;           /* temporary value */
1096
1097         if (m->m_dim < 2) {
1098                 res = matalloc(m->m_size);
1099                 *res = *m;
1100                 val = m->m_table;
1101                 vres = res->m_table;
1102                 for (cur = m->m_size; cur > 0; cur--)
1103                         invertvalue(val++, vres++);
1104                 return res;
1105         }
1106         if (m->m_dim != 2) {
1107                 math_error("Matrix dimension exceeds two for inverse");
1108                 /*NOTREACHED*/
1109         }
1110         if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) {
1111                 math_error("Inverting non-square matrix");
1112                 /*NOTREACHED*/
1113         }
1114         /*
1115          * Begin by creating the identity matrix with the same attributes.
1116          */
1117         res = matalloc(m->m_size);
1118         *res = *m;
1119         rows = (m->m_max[0] - m->m_min[0] + 1);
1120         val = res->m_table;
1121         for (row = 0; row < rows; row++) {
1122                 for (col = 0; col < rows; col++) {
1123                         if (row == col)
1124                                 val->v_num = qlink(&_qone_);
1125                         else
1126                                 val->v_num = qlink(&_qzero_);
1127                         val->v_type = V_NUM;
1128                         val++;
1129                 }
1130         }
1131         /*
1132          * Now loop over each row, and eliminate all entries in the
1133          * corresponding column by using row operations.  Do the same
1134          * operations on the resulting matrix.  Copy the original matrix
1135          * so that we don't destroy it.
1136          */
1137         m = matcopy(m);
1138         for (cur = 0; cur < rows; cur++) {
1139                 /*
1140                  * Find the first nonzero value in the rest of the column
1141                  * downwards from [cur,cur].  If there is no such value, then
1142                  * the matrix is not invertible.  If the first nonzero entry
1143                  * is not the current row, then swap the two rows to make the
1144                  * current one nonzero.
1145                  */
1146                 row = cur;
1147                 val = &m->m_table[(row * rows) + row];
1148                 while (testvalue(val) == 0) {
1149                         if (++row >= rows) {
1150                                 matfree(m);
1151                                 matfree(res);
1152                                 math_error("Matrix is not invertible");
1153                                 /*NOTREACHED*/
1154                         }
1155                         val += rows;
1156                 }
1157                 invertvalue(val, &mulval);
1158                 if (row != cur) {
1159                         matswaprow(m, row, cur);
1160                         matswaprow(res, row, cur);
1161                 }
1162                 /*
1163                  * Now for every other nonzero entry in the current column,
1164                  * subtract the appropriate multiple of the current row to
1165                  * force that entry to become zero.
1166                  */
1167                 val = &m->m_table[cur];
1168                 for (row = 0; row < rows; row++) {
1169                         if ((row == cur) || (testvalue(val) == 0)) {
1170                                 if (row+1 < rows)
1171                                         val += rows;
1172                                 continue;
1173                         }
1174                         mulvalue(val, &mulval, &tmpval);
1175                         matsubrow(m, row, cur, &tmpval);
1176                         matsubrow(res, row, cur, &tmpval);
1177                         freevalue(&tmpval);
1178                         if (row+1 < rows)
1179                                 val += rows;
1180                 }
1181                 freevalue(&mulval);
1182         }
1183         /*
1184          * Now the original matrix has nonzero entries only on its main
1185          * diagonal.  Scale the rows of the result matrix by the inverse
1186          * of those entries.
1187          */
1188         val = m->m_table;
1189         for (row = 0; row < rows; row++) {
1190                 if ((val->v_type != V_NUM) || !qisone(val->v_num)) {
1191                         invertvalue(val, &mulval);
1192                         matmulrow(res, row, &mulval);
1193                         freevalue(&mulval);
1194                 }
1195                 if (row+1 < rows)
1196                         val += (rows + 1);
1197         }
1198         matfree(m);
1199         return res;
1200 }
1201
1202
1203 /*
1204  * Calculate the determinant of a square matrix.
1205  * This uses the fraction-free Gauss-Bareiss algorithm.
1206  */
1207 VALUE
1208 matdet(MATRIX *m)
1209 {
1210         long n;                 /* original matrix is n x n */
1211         long k;                 /* working submatrix is k x k */
1212         long i, j;
1213         VALUE *pivot, *div, *val;
1214         VALUE *vp, *vv;
1215         VALUE tmp1, tmp2, tmp3;
1216         BOOL neg;               /* whether to negate determinant */
1217
1218         if (m->m_dim < 2) {
1219                 vp = m->m_table;
1220                 i = m->m_size;
1221                 copyvalue(vp, &tmp1);
1222
1223                 while (--i > 0) {
1224                         mulvalue(&tmp1, ++vp, &tmp2);
1225                         freevalue(&tmp1);
1226                         tmp1 = tmp2;
1227                 }
1228                 return tmp1;
1229         }
1230
1231         if (m->m_dim != 2)
1232                 return error_value(E_DET2);
1233         if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
1234                 return error_value(E_DET3);
1235
1236         /*
1237          * Loop over each row, and eliminate all lower entries in the
1238          * corresponding column by using row operations.  Copy the original
1239          * matrix so that we don't destroy it.
1240          */
1241         neg = FALSE;
1242         m = matcopy(m);
1243         n = (m->m_max[0] - m->m_min[0] + 1);
1244         pivot = div = m->m_table;
1245         for (k = n; k > 0; k--) {
1246                 /*
1247                  * Find the first nonzero value in the rest of the column
1248                  * downwards from pivot.  If there is no such value, then
1249                  * the determinant is zero.  If the first nonzero entry is not
1250                  * the pivot, then swap rows in the k * k submatrix, and
1251                  * remember that the determinant changes sign.
1252                  */
1253                 val = pivot;
1254                 i = k;
1255                 while (!testvalue(val)) {
1256                         if (--i <= 0) {
1257                                 tmp1.v_type = V_NUM;
1258                                 tmp1.v_subtype = V_NOSUBTYPE;
1259                                 tmp1.v_num = qlink(&_qzero_);
1260                                 matfree(m);
1261                                 return tmp1;
1262                         }
1263                         val += n;
1264                 }
1265                 if (i < k) {
1266                         vp = pivot;
1267                         vv = val;
1268                         j = k;
1269                         while (j-- > 0) {
1270                                 tmp1 = *vp;
1271                                 *vp++ = *vv;
1272                                 *vv++ = tmp1;
1273                         }
1274                         neg = !neg;
1275                 }
1276                 /*
1277                  * Now for every val below the pivot, for each entry to
1278                  * the right of val, calculate the 2 x 2 determinant
1279                  * with corners at the pivot and the entry.  If
1280                  * k < n, divide by div (the previous pivot value).
1281                  */
1282                 val = pivot;
1283                 i = k;
1284                 while (--i > 0) {
1285                         val += n;
1286                         vp = pivot;
1287                         vv = val;
1288                         j = k;
1289                         while (--j > 0) {
1290                                 mulvalue(pivot, ++vv, &tmp1);
1291                                 mulvalue(val, ++vp, &tmp2);
1292                                 subvalue(&tmp1, &tmp2, &tmp3);
1293                                 freevalue(&tmp1);
1294                                 freevalue(&tmp2);
1295                                 freevalue(vv);
1296                                 if (k < n) {
1297                                         divvalue(&tmp3, div, vv);
1298                                         freevalue(&tmp3);
1299                                 }
1300                                 else
1301                                         *vv = tmp3;
1302                         }
1303                 }
1304                 div = pivot;
1305                 if (k > 0)
1306                         pivot += n + 1;
1307         }
1308         if (neg)
1309                 negvalue(div, &tmp1);
1310         else
1311                 copyvalue(div, &tmp1);
1312         matfree(m);
1313         return tmp1;
1314 }
1315
1316
1317 /*
1318  * Local utility routine to swap two rows of a square matrix.
1319  * No checks are made to verify the legality of the arguments.
1320  */
1321 S_FUNC void
1322 matswaprow(MATRIX *m, long r1, long r2)
1323 {
1324         register VALUE *v1, *v2;
1325         register long rows;
1326         VALUE tmp;
1327
1328         if (r1 == r2)
1329                 return;
1330         rows = (m->m_max[0] - m->m_min[0] + 1);
1331         v1 = &m->m_table[r1 * rows];
1332         v2 = &m->m_table[r2 * rows];
1333         while (rows-- > 0) {
1334                 tmp = *v1;
1335                 *v1 = *v2;
1336                 *v2 = tmp;
1337                 v1++;
1338                 v2++;
1339         }
1340 }
1341
1342
1343 /*
1344  * Local utility routine to subtract a multiple of one row to another one.
1345  * The row to be changed is oprow, the row to be subtracted is baserow.
1346  * No checks are made to verify the legality of the arguments.
1347  */
1348 S_FUNC void
1349 matsubrow(MATRIX *m, long oprow, long baserow, VALUE *mulval)
1350 {
1351         register VALUE *vop, *vbase;
1352         register long entries;
1353         VALUE tmp1, tmp2;
1354
1355         entries = (m->m_max[0] - m->m_min[0] + 1);
1356         vop = &m->m_table[oprow * entries];
1357         vbase = &m->m_table[baserow * entries];
1358         while (entries-- > 0) {
1359                 mulvalue(vbase, mulval, &tmp1);
1360                 subvalue(vop, &tmp1, &tmp2);
1361                 freevalue(&tmp1);
1362                 freevalue(vop);
1363                 *vop = tmp2;
1364                 vop++;
1365                 vbase++;
1366         }
1367 }
1368
1369
1370 /*
1371  * Local utility routine to multiply a row by a specified number.
1372  * No checks are made to verify the legality of the arguments.
1373  */
1374 S_FUNC void
1375 matmulrow(MATRIX *m, long row, VALUE *mulval)
1376 {
1377         register VALUE *val;
1378         register long rows;
1379         VALUE tmp;
1380
1381         rows = (m->m_max[0] - m->m_min[0] + 1);
1382         val = &m->m_table[row * rows];
1383         while (rows-- > 0) {
1384                 mulvalue(val, mulval, &tmp);
1385                 freevalue(val);
1386                 *val = tmp;
1387                 val++;
1388         }
1389 }
1390
1391
1392 /*
1393  * Make a full copy of a matrix.
1394  */
1395 MATRIX *
1396 matcopy(MATRIX *m)
1397 {
1398         MATRIX *res;
1399         register VALUE *v1, *v2;
1400         register long i;
1401
1402         res = matalloc(m->m_size);
1403         *res = *m;
1404         v1 = m->m_table;
1405         v2 = res->m_table;
1406         i = m->m_size;
1407         while (i-- > 0) {
1408                 copyvalue(v1, v2);
1409                 v1++;
1410                 v2++;
1411         }
1412         return res;
1413 }
1414
1415
1416 /*
1417  * Make a matrix the same size as another and filled with a fixed value.
1418  *
1419  * given:
1420  *      m               matrix to initialize
1421  *      v1              value to fill most of matrix with
1422  *      v2              value for diagonal entries (or NULL)
1423  */
1424 MATRIX *
1425 matinit(MATRIX *m, VALUE *v1, VALUE *v2)
1426 {
1427         MATRIX *res;
1428         register VALUE *v;
1429         register long i;
1430         long row;
1431         long rows;
1432
1433         /*
1434          * clone matrix size
1435          */
1436         res = matalloc(m->m_size);
1437         *res = *m;
1438
1439         /*
1440          * firewall
1441          */
1442         if (v2 && ((res->m_dim != 2) ||
1443                 ((res->m_max[0] - res->m_min[0]) !=
1444                  (res->m_max[1] - res->m_min[1])))) {
1445                         math_error("Filling diagonals of non-square matrix");
1446                         /*NOTREACHED*/
1447         }
1448
1449         /*
1450          * fill the bulk of the matrix
1451          */
1452         v = res->m_table;
1453         if (v2 == NULL) {
1454                 i = m->m_size;
1455                 while (i-- > 0) {
1456                         copyvalue(v1, v++);
1457                 }
1458                 return res;
1459         }
1460
1461         /*
1462          * fill the diaginal of a square matrix if requested
1463          */
1464         rows = res->m_max[0] - res->m_min[0] + 1;
1465         v = res->m_table;
1466         for (row = 0; row < rows; row++) {
1467                 copyvalue(v2, v+row);
1468                 v += rows;
1469         }
1470         return res;
1471 }
1472
1473
1474 /*
1475  * Allocate a matrix with the specified number of elements.
1476  */
1477 MATRIX *
1478 matalloc(long size)
1479 {
1480         MATRIX *m;
1481         long i;
1482         VALUE *vp;
1483
1484         m = (MATRIX *) malloc(matsize(size));
1485         if (m == NULL) {
1486                 math_error("Cannot get memory to allocate matrix of size %ld",
1487                            size);
1488                 /*NOTREACHED*/
1489         }
1490         m->m_size = size;
1491         for (i = size, vp = m->m_table; i > 0; i--, vp++)
1492                 vp->v_subtype = V_NOSUBTYPE;
1493         return m;
1494 }
1495
1496
1497 /*
1498  * Free a matrix, along with all of its element values.
1499  */
1500 void
1501 matfree(MATRIX *m)
1502 {
1503         register VALUE *vp;
1504         register long i;
1505
1506         vp = m->m_table;
1507         i = m->m_size;
1508         while (i-- > 0)
1509                 freevalue(vp++);
1510         free(m);
1511 }
1512
1513
1514 /*
1515  * Test whether a matrix has any "nonzero" values.
1516  * Returns TRUE if so.
1517  */
1518 BOOL
1519 mattest(MATRIX *m)
1520 {
1521         register VALUE *vp;
1522         register long i;
1523
1524         vp = m->m_table;
1525         i = m->m_size;
1526         while (i-- > 0) {
1527                 if (testvalue(vp++))
1528                         return TRUE;
1529         }
1530         return FALSE;
1531 }
1532
1533 /*
1534  * Sum the elements in a matrix.
1535  */
1536 void
1537 matsum(MATRIX *m, VALUE *vres)
1538 {
1539         VALUE *vp;
1540         VALUE tmp;              /* first sum value */
1541         VALUE sum;              /* final sum value */
1542         long i;
1543
1544         vp = m->m_table;
1545         i = m->m_size;
1546         copyvalue(vp, &sum);
1547
1548         while (--i > 0) {
1549                 addvalue(&sum, ++vp, &tmp);
1550                 freevalue(&sum);
1551                 sum = tmp;
1552         }
1553         *vres = sum;
1554 }
1555
1556
1557 /*
1558  * Test whether or not two matrices are equal.
1559  * Equality is determined by the shape and values of the matrices,
1560  * but not by their index bounds.  Returns TRUE if they differ.
1561  */
1562 BOOL
1563 matcmp(MATRIX *m1, MATRIX *m2)
1564 {
1565         VALUE *v1, *v2;
1566         long i;
1567
1568         if (m1 == m2)
1569                 return FALSE;
1570         if ((m1->m_dim != m2->m_dim) || (m1->m_size != m2->m_size))
1571                 return TRUE;
1572         for (i = 0; i < m1->m_dim; i++) {
1573                 if ((m1->m_max[i] - m1->m_min[i]) !=
1574                     (m2->m_max[i] - m2->m_min[i]))
1575                 return TRUE;
1576         }
1577         v1 = m1->m_table;
1578         v2 = m2->m_table;
1579         i = m1->m_size;
1580         while (i-- > 0) {
1581                 if (comparevalue(v1++, v2++))
1582                         return TRUE;
1583         }
1584         return FALSE;
1585 }
1586
1587
1588 void
1589 matreverse(MATRIX *m)
1590 {
1591         VALUE *p, *q;
1592         VALUE tmp;
1593
1594         p = m->m_table;
1595         q = m->m_table + m->m_size - 1;
1596         while (q > p) {
1597                 tmp = *p;
1598                 *p++ = *q;
1599                 *q-- = tmp;
1600         }
1601 }
1602
1603
1604 void
1605 matsort(MATRIX *m)
1606 {
1607         VALUE *a, *b, *next, *end;
1608         VALUE *buf, *p;
1609         VALUE *S[LONG_BITS];
1610         long len[LONG_BITS];
1611         long i, j, k;
1612
1613         buf = (VALUE *) malloc(m->m_size * sizeof(VALUE));
1614         if (buf == NULL) {
1615                 math_error("Not enough memory for matsort");
1616                 /*NOTREACHED*/
1617         }
1618         next = m->m_table;
1619         end = next + m->m_size;
1620         for (k = 0; next && k < LONG_BITS; k++) {
1621                 S[k] = next++;                  /* S[k] is start of a run */
1622                 len[k] = 1;
1623                 if (next == end)
1624                         next = NULL;
1625                 while (k > 0 && (!next || len[k] >= len[k - 1])) {/* merging */
1626                         j = len[k];
1627                         b = S[k--];
1628                         i = len[k];
1629                         a = S[k];
1630                         len[k] += j;
1631                         p = buf;
1632                         if (precvalue(b, a)) {
1633                                 do {
1634                                         *p++ = *b++;
1635                                         j--;
1636                                 } while (j > 0 && precvalue(b,a));
1637                                 if (j == 0) {
1638                                         memcpy(p, a, i * sizeof(VALUE));
1639                                         memcpy(S[k], buf,
1640                                                len[k] * sizeof(VALUE));
1641                                         continue;
1642                                 }
1643                         }
1644
1645                         do {
1646                                 do {
1647                                         *p++ = *a++;
1648                                         i--;
1649                                 } while (i > 0 && !precvalue(b,a));
1650                                 if (i == 0) {
1651                                         break;
1652                                 }
1653                                 do {
1654                                         *p++ = *b++;
1655                                         j--;
1656                                 } while (j > 0 && precvalue(b,a));
1657                         } while (j != 0);
1658
1659                         if (i == 0) {
1660                                 memcpy(S[k], buf, (p - buf) * sizeof(VALUE));
1661                         } else if (j == 0) {
1662                                 memcpy(p, a, i * sizeof(VALUE));
1663                                 memcpy(S[k], buf, len[k] * sizeof(VALUE));
1664                         }
1665                 }
1666         }
1667         free(buf);
1668         if (k >= LONG_BITS) {
1669                 /* this should never happen */
1670                 math_error("impossible k overflow in matsort!");
1671                 /*NOTREACHED*/
1672         }
1673 }
1674
1675 void
1676 matrandperm(MATRIX *m)
1677 {
1678         VALUE *vp;
1679         long s, i;
1680         VALUE val;
1681
1682         s = m->m_size;
1683         for (vp = m->m_table; s > 1; vp++, s--) {
1684                 i = irand(s);
1685                 if (i > 0) {
1686                         val = *vp;
1687                         *vp = vp[i];
1688                         vp[i] = val;
1689                 }
1690         }
1691 }
1692
1693
1694 /*
1695  * Test whether or not a matrix is the identity matrix.
1696  * Returns TRUE if so.
1697  */
1698 BOOL
1699 matisident(MATRIX *m)
1700 {
1701         register VALUE *val;    /* current value */
1702         long row, col;          /* row and column numbers */
1703
1704         val = m->m_table;
1705         if (m->m_dim == 0) {
1706                 return (val->v_type == V_NUM && qisone(val->v_num));
1707         }
1708         if (m->m_dim == 1) {
1709                 for (row = m->m_min[0]; row <= m->m_max[0]; row++, val++) {
1710                         if (val->v_type != V_NUM || !qisone(val->v_num))
1711                                 return FALSE;
1712                 }
1713                 return TRUE;
1714         }
1715         if ((m->m_dim != 2) ||
1716                 ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])))
1717                         return FALSE;
1718         for (row = m->m_min[0]; row <= m->m_max[0]; row++) {
1719                 /*
1720                  * We could use col = m->m_min[1]; col < m->m_max[1]
1721                  * but if m->m_min[0] != m->m_min[1] this won't work.
1722                  * We know that we have a square 2-dimensional matrix
1723                  * so we will pretend that m->m_min[0] == m->m_min[1].
1724                  */
1725                 for (col = m->m_min[0]; col <= m->m_max[0]; col++) {
1726                         if (val->v_type != V_NUM)
1727                                 return FALSE;
1728                         if (row == col) {
1729                                 if (!qisone(val->v_num))
1730                                         return FALSE;
1731                         } else {
1732                                 if (!qiszero(val->v_num))
1733                                         return FALSE;
1734                         }
1735                         val++;
1736                 }
1737         }
1738         return TRUE;
1739 }
1740
1741
1742 /*
1743  * Print a matrix and possibly few of its elements.
1744  * The argument supplied specifies how many elements to allow printing.
1745  * If any elements are printed, they are printed in short form.
1746  */
1747 void
1748 matprint(MATRIX *m, long max_print)
1749 {
1750         VALUE *vp;
1751         long fullsize, count, index;
1752         long dim, i, j;
1753         char *msg;
1754         long sizes[MAXDIM];
1755
1756         dim = m->m_dim;
1757         fullsize = 1;
1758         for (i = dim - 1; i >= 0; i--) {
1759                 sizes[i] = fullsize;
1760                 fullsize *= (m->m_max[i] - m->m_min[i] + 1);
1761         }
1762         msg = ((max_print > 0) ? "\nmat [" : "mat [");
1763         if (dim) {
1764                 for (i = 0; i < dim; i++) {
1765                         if (m->m_min[i]) {
1766                                 math_fmt("%s%ld:%ld", msg,
1767                                          m->m_min[i], m->m_max[i]);
1768                         } else {
1769                                 math_fmt("%s%ld", msg, m->m_max[i] + 1);
1770                         }
1771                         msg = ",";
1772                 }
1773         } else {
1774                 math_str("mat [");
1775         }
1776         if (max_print > fullsize) {
1777                 max_print = fullsize;
1778         }
1779         vp = m->m_table;
1780         count = 0;
1781         for (index = 0; index < fullsize; index++) {
1782                 if ((vp->v_type != V_NUM) || !qiszero(vp->v_num))
1783                         count++;
1784                 vp++;
1785         }
1786         math_fmt("] (%ld element%s, %ld nonzero)",
1787                 fullsize, (fullsize == 1) ? "" : "s", count);
1788         if (max_print <= 0)
1789                 return;
1790
1791         /*
1792          * Now print the first few elements of the matrix in short
1793          * and unambigous format.
1794          */
1795         math_str(":\n");
1796         vp = m->m_table;
1797         for (index = 0; index < max_print; index++) {
1798                 msg = "  [";
1799                 j = index;
1800                 if (dim) {
1801                         for (i = 0; i < dim; i++) {
1802                                 math_fmt("%s%ld", msg,
1803                                          m->m_min[i] + (j / sizes[i]));
1804                                 j %= sizes[i];
1805                                 msg = ",";
1806                         }
1807                 } else {
1808                         math_str(msg);
1809                 }
1810                 math_str("] = ");
1811                 printvalue(vp++, PRINT_SHORT | PRINT_UNAMBIG);
1812                 math_str("\n");
1813         }
1814         if (max_print < fullsize)
1815                 math_str("  ...\n");
1816 }