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