1 /* file: matrix.c */ /* AK 091086 */
2 #include "def.h"
3 #include "macro.h"
4
5 #ifdef MATRIXTRUE
6 static struct matrix * callocmatrix();
7 static INT scan_matrix_co();
8
transform_matrix(a,f,b)9 static INT transform_matrix(a,f,b) OP a,b; INT (*f)();
10 {
11 INT e = 0L;
12 INT i,j,erg = OK;
13 CTO(MATRIX,"transform_matrix(1)",a);
14
15 if (a==b)
16 {
17 OP c = callocobject();
18 *c = *a;
19 C_O_K(b,EMPTY);
20 e = 1L;
21 a = c;
22 }
23 m_ilih_m(S_M_LI(a),S_M_HI(a),b);
24 for (i=0L;i<S_M_HI(b);i++)
25 for (j=0L;j<S_M_LI(b);j++)
26 erg += (*f)(S_M_IJ(a,i,j),S_M_IJ(b,i,j));
27 if (e==1L)
28 erg += freeall(a);
29 ENDR("internal function:transform_matrix");
30 }
31
cast_apply_matrix(a)32 INT cast_apply_matrix(a) OP a;
33 /* AK 270295 */
34 {
35 INT i,ml,erg = OK,j;
36 OP b;
37 EOP("cast_apply_matrix(1)",a);
38
39 if (S_O_K(a) == MATRIX)
40 goto camb;
41 else if (S_O_K(a) == KOSTKA)
42 goto camb;
43 else if (S_O_K(a) == VECTOR)
44 {
45 ml = 0;
46 for (i=(INT)0;i<S_V_LI(a);i++)
47 {
48 if (not VECTORP(S_V_I(a,i)))
49 goto came;
50 if (S_V_LI(S_V_I(a,i)) > ml)
51 ml = S_V_LI(S_V_I(a,i));
52 }
53
54 /* now vector of vector == we can cast */
55 b = callocobject();
56 *b = *a;
57 C_O_K(a,EMPTY);
58 erg += m_ilih_m(ml,S_V_LI(b),a);
59 for (i=0;i<S_M_HI(a);i++)
60 for (j=0;j<S_V_LI(S_V_I(b,i));j++)
61 {
62 * S_M_IJ(a,i,j) = * S_V_I(S_V_I(b,i),j);
63 C_O_K(S_V_I(S_V_I(b,i),j),EMPTY);
64 }
65 erg += freeall(b);
66 goto camb;
67 }
68 came:
69 printobjectkind(a);
70 erg += error("cast_apply_matrix: transfer not possible");
71 camb:
72 ENDR("cast_apply_matrix");
73 }
74
mem_size_matrix(a)75 INT mem_size_matrix(a) OP a;
76 /* AK 150295 */
77 {
78 INT erg = 0,i,j;
79 OP z;
80 if (a == NULL) return 0;
81
82 if (not MATRIXP(a)) WTO("mem_size_matrix",a);
83 erg += sizeof(struct object);
84 erg += sizeof(struct matrix);
85 erg += mem_size(S_M_H(a));
86 erg += mem_size(S_M_L(a));
87 for (i=S_M_HI(a)*S_M_LI(a)-1,z=S_M_S(a);i>=0;i--,z++)
88 erg += mem_size(z);
89 return erg;
90 }
91
92
93
mod_matrix(a,b,c)94 INT mod_matrix(a,b,c) OP a,b,c;
95 /* AK 300793 */
96 {
97 INT erg = OK;
98 INT i,j;
99 CTO(MATRIX,"mod_matrix(1)",a);
100 CTO(INTEGER,"mod_matrix(2)",b);
101 CTO(EMPTY,"mod_matrix(3)",c);
102
103 erg += m_ilih_m(S_M_LI(a),S_M_HI(a),c);
104 for (i=0L;i<S_M_HI(a);i++)
105 for (j=0L;j<S_M_LI(a);j++)
106 erg += mod(S_M_IJ(a,i,j),b,S_M_IJ(c,i,j));
107 ENDR("mod_matrix");
108 }
109
matrixp(a)110 INT matrixp(a) OP a;
111 /* AK 280193 */
112 {
113 if (S_O_K(a) == MATRIX) return TRUE;
114 if (S_O_K(a) == KRANZTYPUS) return TRUE;
115 if (S_O_K(a) == INTEGERMATRIX) return TRUE;
116 return FALSE;
117 }
118
119
120
addinvers_matrix(a,b)121 INT addinvers_matrix(a,b) OP a,b;
122 /* AK 240293 */
123 /* AK 060704 V3.0 */
124 {
125 return transform_matrix(a,addinvers,b);
126 }
127
absolute_matrix(a,b)128 INT absolute_matrix(a,b) OP a,b;
129 /* AK 240293 */
130 /* AK 060704 V3.0 */
131 {
132 return transform_matrix(a,absolute,b);
133 }
134
sum_matrix(a,b)135 INT sum_matrix(a,b) OP a,b;
136 /* AK 060704 V3.0 */
137 {
138 INT erg = OK;
139 CTTO(INTEGERMATRIX,MATRIX,"sum_matrix(1)",a);
140 {
141 INT i;
142 OP z = S_M_S(a);
143 CLEVER_COPY(S_M_IJ(a,S_M_HI(a)-1,S_M_LI(a)-1),b);
144 for (i=S_M_HI(a)*S_M_LI(a);i>1;i--,z++)
145 ADD_APPLY(z, b);
146 }
147 ENDR("sum_matrix");
148 }
149
150
nullp_integermatrix(a)151 INT nullp_integermatrix(a) OP a;
152 /* AK 150802 */
153 /* AK 060704 V3.0 */
154 {
155 INT erg = OK;
156 CTO(INTEGERMATRIX,"nullp_integermatrix(1)",a);
157 {
158 INT i,j;
159 for (i=0L;i<S_M_HI(a);i++)
160 for (j=0L;j<S_M_LI(a);j++)
161 if (not NULLP_INTEGER(S_M_IJ(a,i,j))) return FALSE;
162 return TRUE;
163 }
164 ENDR("nullp_integermatrix");
165 }
166
167
nullp_matrix(a)168 INT nullp_matrix(a) OP a;
169 /* AK 110992 */
170 {
171 INT i,j;
172 for (i=0L;i<S_M_HI(a);i++)
173 for (j=0L;j<S_M_LI(a);j++)
174 if (not NULLP(S_M_IJ(a,i,j))) return FALSE;
175 return TRUE;
176 }
177
178
179
180
einsp_matrix(a)181 INT einsp_matrix(a) OP a;
182 /* AK 110992 */
183 {
184 INT i,j;
185 INT erg = OK;
186 CTTO(INTEGERMATRIX,MATRIX,"einsp_matrix(1)",a);
187 if (S_M_HI(a) != S_M_LI(a)) return FALSE;
188 for (i=0L;i<S_M_HI(a);i++)
189 for (j=0L;j<S_M_HI(a);j++)
190 if (i==j)
191 {
192 if (not EINSP(S_M_IJ(a,i,j))) return FALSE;
193 }
194 else {
195 if (not NULLP(S_M_IJ(a,i,j))) return FALSE;
196 }
197 return TRUE;
198 ENDR("einsp_matrix");
199 }
200
201
delete_row_matrix(a,index,b)202 INT delete_row_matrix(a,index,b) INT index; OP a,b;
203 /* AK 270789 */ /* AK 111289 V1.1 */ /* AK 110691 V1.2 */
204 /* AK 070891 V1.3 */ /* AK 201204 V3.0 */
205 /* AK 021106 V3.1 */
206 {
207 INT i,j;
208 INT erg = OK;
209 SYMCHECK(index >= S_M_HI(a), "delete_row_matrix: index too big");
210 SYMCHECK(index < 0,"delete_row_matrix: index < 0");
211
212
213 if (a==b) { /* AK 201204 */
214 OP z,w;
215 for (z=S_M_IJ(a,index,0),i=0;i<S_M_LI(a);i++,z++)
216 FREESELF(z);
217 z = S_M_IJ(a,index,0);
218 w = S_M_IJ(a,index+1,0);
219 for (i=index+1;i<S_M_HI(a);i++)
220 for (j=0;j<S_M_LI(a);j++,z++,w++) SWAP(z,w);
221 // now there is a row with empty objects at the end
222 C_M_S(a, (OP) SYM_realloc((char*)S_M_S(a),
223 (int) (S_M_HI(a)-1)*S_M_LI(a)*sizeof(struct object)
224 )
225 );
226 C_M_HASH(a,-1);
227 M_I_I(S_M_HI(a)-1,S_M_H(a));
228 goto endr_ende;
229 }
230
231 erg += m_ilih_m(S_M_LI(a),S_M_HI(a)-1L,b);
232 C_O_K(b,S_O_K(a));
233 for (i=0;i<index;i++)
234 for (j=0;j<S_M_LI(a);j++)
235 COPY(S_M_IJ(a,i,j),S_M_IJ(b,i,j));
236 for (i=index+1;i<S_M_HI(a);i++)
237 for (j=0;j<S_M_LI(a);j++)
238 COPY(S_M_IJ(a,i,j),S_M_IJ(b,i-1L,j));
239 ENDR("delete_row_matrix");
240 }
241
242
243
append_column_matrix(a,b,c)244 INT append_column_matrix(a,b,c) OP a,b,c;
245 /* AK 240206 V3.0 */
246 /* appends the vector object b as a last column to the matrix a */
247 {
248 INT erg = OK;
249 CTTO(MATRIX,INTEGERMATRIX,"append_column_matrix(1)",a);
250 CTTO(VECTOR,INTEGERVECTOR,"append_column_matrix(2)",b);
251 SYMCHECK(S_M_HI(a) != S_V_LI(b),"append_column_matrix: vector of wrong length");
252 {
253 INT i,j;
254 if (a==c) {
255 OP d = CALLOCOBJECT();
256 *d = *a;
257 C_O_K(c,EMPTY);
258 erg += append_column_matrix(d,b,c);
259 FREEALL(d);
260 goto endr_ende;
261 }
262 erg += m_ilih_m(S_M_LI(a)+1L,S_M_HI(a),c);
263 C_O_K(c,S_O_K(a));
264 for (i=0;i<S_M_HI(a);i++)
265 for (j=0;j<S_M_LI(a);j++) COPY(S_M_IJ(a,i,j),S_M_IJ(c,i,j));
266 for (i=0;i<S_M_HI(a);i++) COPY(S_V_I(b,i),S_M_IJ(c,i,S_M_LI(c)-1));
267 }
268 ENDR("append_column_matrix");
269 }
270
delete_column_matrix(a,index,b)271 INT delete_column_matrix(a,index,b) INT index; OP a,b;
272 /* deletes the column with index i
273 the result b is of the same type like a */
274 /* AK 270789 */ /* AK 310790 V1.1 */ /* AK 110691 V1.2 */
275 /* AK 070891 V1.3 */
276 /* AK 240206 V3.0 */
277 {
278 INT erg = OK;
279 CTTO(MATRIX,INTEGERMATRIX,"delete_column_matrix(1)",a);
280 SYMCHECK(index >= S_M_LI(a),"delete_column_matrix: index too big");
281 SYMCHECK(index <0,"delete_column_matrix: index < 0");
282 {
283 INT i,j;
284
285
286 if (a==b) {
287 OP c = CALLOCOBJECT();
288 *c = *b;
289 C_O_K(b,EMPTY);
290 erg += delete_column_matrix(c,index,b);
291 FREEALL(c);
292 goto endr_ende;
293 }
294
295 erg += m_ilih_m(S_M_LI(a)-1L,S_M_HI(a),b);
296 C_O_K(b,S_O_K(a));
297 for (j=0;j<index;j++)
298 for (i=0;i<S_M_HI(a);i++)
299 COPY(S_M_IJ(a,i,j),S_M_IJ(b,i,j));
300 for (j=index+1;j<S_M_LI(a);j++)
301 for (i=0;i<S_M_HI(a);i++)
302 COPY(S_M_IJ(a,i,j),S_M_IJ(b,i,j-1));
303 }
304
305 ENDR("delete_column_matrix");
306 }
307
308
309
det_matrix(a,b)310 INT det_matrix(a,b) OP a,b;
311 /* AK 151289 V1.1 */ /* AK 110691 V1.2 */
312 /* hier wird zum konkreten algorithmus geschaltet */
313 /* AK 210891 V1.3 */
314 /* AK 141098 V2.0 */
315 {
316 if (not quadraticp(a))
317 {
318 error("det_matrix: not quadratic matrix");
319 return(ERROR);
320 }
321 return det_mat_tri(a,b);
322 }
323
324
325
det_mat_tri(a,res)326 INT det_mat_tri(a,res) OP a,res;
327 /* AK Fri Jan 20 12:29:58 MEZ 1989
328 algorithmus zur berechnung der determinante einer matrix a
329 mittels triangulation
330 ergebnis in res
331 algorithmus 41 in CACM */
332 /* verbessert 1963 */
333 /* immer noch fehler bei zeilenvertauschung */
334 /* AK 310790 V1.1 */ /* AK 110691 V1.2 */
335 /* AK 210891 V1.3 */
336 {
337 INT r,i,j,y,count,sign=1L,n;
338 INT erg = OK;
339 OP b, factor,temp;
340
341 CTO(MATRIX,"det_mat_tri(1)",a);
342
343 n = S_M_LI(a);
344 b = callocobject();
345
346 erg += m_i_i(1L,res);
347 factor = callocobject();
348 temp = callocobject();
349
350 erg += copy(a,b);
351
352 for (r=1;r < n;r++)
353 {
354 count = r-1;
355 if (not nullp(S_M_IJ(b,r-1L,r-1))) goto det_mat_tri_resume;
356 det_mat_tri_zerocheck:
357 if (count < (n-1)) count++;
358 else goto det_mat_tri_zero;
359
360 if (not nullp(S_M_IJ(b,count,r-1))) {
361 for (y=r; y<=n; y++)
362 {
363 erg += swap(S_M_IJ(b,count,y-1),S_M_IJ(b,r-1L,y-1));
364 }
365 sign = -sign; goto det_mat_tri_resume;
366 }
367
368
369 goto det_mat_tri_zerocheck;
370 det_mat_tri_zero:
371 erg += m_i_i((INT)0,res);
372 goto det_mat_tri_return;
373 det_mat_tri_resume:
374 for (i=r+1; i<=n; i++)
375 {
376
377 erg += div(S_M_IJ(b,i-1L,r-1),S_M_IJ(b,r-1L,r-1),factor);
378 for (j=r+1;j<=n;j++)
379 {
380 erg += mult(factor,S_M_IJ(b,r-1L,j-1),temp);
381 erg += addinvers_apply(temp);
382 erg += add_apply(temp,S_M_IJ(b,i-1L,j-1L));
383 }
384 }
385 }
386
387 for (i=1;i<=n;i++)
388 erg += mult_apply(S_M_IJ(b,i-1L,i-1),res);
389
390
391 if (sign == -1L)
392 erg += addinvers_apply(res);
393
394 det_mat_tri_return:
395 FREEALL3(temp,factor,b);
396 ENDR("det_mat_tri");
397 }
398
399
400
m_ilih_nm(l,h,m)401 INT m_ilih_nm(l,h,m) OP m; INT l,h;
402 /* mit 0 vorbesetzen */
403 /* make_intlength_intheight_null_matrix */
404 {
405 INT i,erg = OK;
406 OP z;
407 SYMCHECK(l < 0,"m_ilih_nm:l<0");
408 SYMCHECK(h < 0,"m_ilih_nm:h<0");
409
410 erg += m_ilih_m(l,h,m);
411 for (z=S_M_S(m), i=S_M_HI(m) * S_M_LI(m); i>0L; i--,z++)
412 M_I_I(0L,z);
413 ENDR("m_ilih_nm");
414 }
415
m_lh_nm(l,h,m)416 INT m_lh_nm(l,h,m) OP l,h,m;
417 /* AK 110691 V1.2 */ /* mit 0 vorbesetzen */
418 /* make_length_height_null_matrix */
419 /* AK 210891 V1.3 */
420 {
421 INT i,erg = OK;
422 OP z;
423 CTO(INTEGER,"m_lh_nm(1)",l);
424 CTO(INTEGER,"m_lh_nm(2)",h);
425 SYMCHECK(S_I_I(l) < 0,"m_lh_nm:l<0");
426 SYMCHECK(S_I_I(h) < 0,"m_lh_nm:h<0");
427
428 erg += m_lh_m(l,h,m);
429 for (z=S_M_S(m), i=S_M_HI(m) * S_M_LI(m); i>0L; i--,z++)
430 M_I_I(0L,z);
431 ENDR("m_lh_nm");
432 }
433
434
b_lh_nm(l,h,m)435 INT b_lh_nm(l,h,m) OP l,h,m;
436 /* AK 110691 V1.2 */
437 /* mit 0 vorbesetzen */
438 /* build_length_height_null_matrix */
439 /* AK 210891 V1.3 */
440 {
441 INT i,erg = OK;
442 OP z;
443 CTO(INTEGER,"b_lh_nm",l);
444 CTO(INTEGER,"b_lh_nm",h);
445 erg += b_lh_m(l,h,m);
446 for (z=S_M_S(m), i=S_M_HI(m) * S_M_LI(m); i>0L; i--,z++)
447 M_I_I(0L,z);
448 ENDR("b_lh_nm");
449 }
450
451
b_lh_m(l,h,m)452 INT b_lh_m(l,h,m) OP l,h,m;
453 /* build_length_height_matrix */
454 /* height und length werden nicht kopiert */
455 /* AK 250590 V1.1 */ /* AK 110691 V1.2 */
456 /* AK 210891 V1.3 */
457 {
458 OP s;
459 INT i;
460 INT erg = OK;
461 CTO(INTEGER,"b_lh_m",l);
462 CTO(INTEGER,"b_lh_m",h);
463 i = S_I_I(l)*S_I_I(h);
464 if (i < 0)
465 {
466 erg += error("b_lh_m:negative values for dimension of a matrix");
467 }
468 else if (i==0)
469 {
470 erg += b_lhs_m(l,h,NULL,m);
471 }
472 else
473 {
474 s = (OP) SYM_malloc(S_I_I(l)*S_I_I(h)*sizeof(struct object));
475 for (i=0L;i<S_I_I(l)*S_I_I(h); i++) C_O_K(s+i,EMPTY);
476 erg += b_lhs_m(l,h,s,m);
477 }
478 ENDR("b_lh_m");
479 }
480
481
482
m_lh_m(len,height,matrix)483 INT m_lh_m(len,height,matrix) OP height, len, matrix;
484 /* make_length_height_matrix */
485 /* height und length werden kopiert */
486 /* AK 041286 */ /* AK 070789 V1.0 */ /* AK 310790 V1.1 */
487 /* AK 260291 V1.2 */ /* AK 210891 V1.3 */
488 /* AK 250398 V2.0 */ /* AK 140205 V3.0 */
489 /* parameter may be equal */
490 {
491 INT erg = OK;
492 INT li,hi;
493 CTO(INTEGER,"m_lh_m",len);
494 CTO(INTEGER,"m_lh_m",height);
495 hi = S_I_I(height);
496 li = S_I_I(len);
497 SYMCHECK(li<0,"m_lh_m:negative length");
498 SYMCHECK(hi<0,"m_lh_m:negative height");
499 SYMCHECK(hi*li<0,"m_lh_m:size too large");
500
501 erg += b_lhs_m(
502 callocobject(),
503 callocobject(),
504 (struct object *) SYM_calloc((int)(hi*li),
505 sizeof(struct object)),
506 matrix);
507 M_I_I(li,S_M_L(matrix));
508 M_I_I(hi,S_M_H(matrix));
509 ENDR("m_lh_m");
510 }
511
512
b_lhs_m(len,height,self,res)513 INT b_lhs_m(len,height,self,res) OP len, height, self, res;
514 /* AK 060789 V1.0 */ /* AK 081289 V1.1 */ /* AK 110691 V1.2 */
515 /* AK 210891 V1.3 */
516 {
517 OBJECTSELF d;
518
519 d.ob_matrix = callocmatrix();
520 b_ks_o(MATRIX, d, res);
521 C_M_L(res,len);
522 C_M_H(res,height);
523 C_M_S(res,self);
524 C_M_HASH(res,-1);
525 return(OK);
526 }
527
hash_matrix(a)528 INT hash_matrix(a) OP a;
529 /* AK 110703 */
530 {
531 INT erg = OK;
532 INT i,j;
533 SYMCHECK(not MATRIXP(a),"hash_matrix: no matrix object");
534 if (S_M_HI(a) == 0) return 4711;
535 if (S_M_LI(a) == 0) return 4711;
536 if (S_M_HASH(a) == -1)
537 {
538 INT res = 1;
539 OP z;
540 FORALL(z,a,
541 {
542 res *= 4711;
543 res += HASH(z);
544 });
545 C_M_HASH(a,res);
546 }
547 return (S_M_HASH(a));
548
549 ENDR("hash_matrix");
550 }
551
eq_matrix(a,b)552 INT eq_matrix(a,b) OP a,b;
553 /* AK 110703 */
554 {
555 INT erg = OK;
556 INT i,j; OP za,zb;
557 SYMCHECK(!MATRIXP(a),"eq_matrix(1): no matrix");
558 if (! MATRIXP(b)) return FALSE;
559 if (S_M_HI(a) != S_M_HI(b)) return FALSE;
560 if (S_M_LI(a) != S_M_LI(b)) return FALSE;
561 if (S_M_HASH(a) != -1)
562 if (S_M_HASH(b) != -1)
563 if (S_M_HASH(a) != S_M_HASH(b))
564 {
565 return FALSE;
566 }
567
568 if (S_O_K(a) == INTEGERMATRIX)
569 if (S_O_K(b) == INTEGERMATRIX)
570 return (comp_integermatrix(a,b)==0);
571
572 for (i=0,za=S_M_S(a),zb=S_M_S(b); i<S_M_HI(a); i++) /* AK 210104 */
573 for (j=0;j<S_M_LI(a);j++,za++,zb++)
574 if (not EQ(za,zb)) return FALSE;
575 return TRUE;
576 ENDR("eq_matrix");
577 }
578
579
580
581
m_ilih_m(len,height,matrix)582 INT m_ilih_m(len,height,matrix) INT height, len; OP matrix;
583 /* AK 090988 */ /* AK 070789 V1.0 */ /* AK 310790 V1.1 */
584 /* AK 210891 V1.3 */
585 {
586 INT erg=OK;
587 SYMCHECK(len < 0, "m_ilih_m:length < 0");
588 SYMCHECK(height < 0, "m_ilih_m:height < 0");
589 if (len*height == 0) /* AK 210802 */
590 {
591 erg += b_lhs_m( CALLOCOBJECT(),
592 CALLOCOBJECT(),
593 NULL,
594 matrix
595 );
596 }
597
598 else {
599 ma: erg = OK;
600 erg += b_lhs_m( CALLOCOBJECT(),
601 CALLOCOBJECT(),
602 (struct object *) SYM_calloc( (int) (height*len),
603 sizeof(struct object)
604 ),
605 matrix
606 );
607
608 if (S_M_S(matrix) == NULL)
609 {
610 INT err;
611 err = error("m_ilih_m:self == NULL ");
612 if (err==ERROR_EXPLAIN)
613 fprintf(stderr,
614 "I wanted a %" PRIINT " x %" PRIINT " matrix" ,
615 len,
616 height
617 );
618 if (err==ERROR_RETRY)
619 goto ma;
620 }
621 }
622 M_I_I(len,S_M_L(matrix));
623 M_I_I(height,S_M_H(matrix));
624 ENDR("m_ilih_m");
625 }
626
627
quadraticp(mat)628 INT quadraticp(mat) OP mat;
629 /* AK 060789 V1.0 */ /* AK 310790 V1.1 */
630 /* AK 210891 V1.3 */
631 {
632 return(S_M_LI(mat) == S_M_HI(mat));
633 }
634
635
636
637
638 #ifdef PERMTRUE
det050995(mat,perm,c)639 static INT det050995(mat,perm,c) OP mat,perm,c;
640 /* brechnet aus Matrix (a_ij) und der Permutation p_1,..,p_n
641 den wert p_1,p_2 * p_3,p_4 * .. */
642 {
643 INT i,erg = OK;
644 CTO(MATRIX,"det050995(1)",mat);
645 CTO(PERMUTATION,"det050995(2)",perm);
646
647 if (neq(S_M_L(mat),S_P_L(perm)))
648 return error("det050995:wrong lengths");
649
650 erg += copy(S_M_IJ(mat,S_P_II(perm,0L)-1,S_P_II(perm,1L)-1L),c);
651
652 for (i=2L;i<S_P_LI(perm);i+=2)
653 {
654 erg += mult_apply(S_M_IJ(mat,S_P_II(perm,i)-1L,S_P_II(perm,i+1)-1L),c);
655 }
656
657 ENDR("internal routine:det050995");
658 }
det270588(mat,perm,c)659 INT det270588(mat,perm,c) OP mat,perm,c;
660 /* AK270588 */
661 /* brechnet aus Matrix (a_ij) und der Permutation p_1,..,p_n
662 den wert a_1,p_1 * a_2,p_2 * .. a_n,p_n */
663 /* damit kann z.b. die determinante berechnet werden */
664 /* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 180691 V1.2 */
665 /* AK 210891 V1.3 */
666 {
667 INT i,erg = OK;
668 CTO(MATRIX,"det270588(1)",mat);
669 CTO(PERMUTATION,"det270588(2)",perm);
670 SYMCHECK(S_M_HI(mat) != S_M_LI(mat), "det270588:not quadratic");
671 SYMCHECK(S_M_LI(mat) != S_P_LI(perm), "det270588:wrong lengths");
672
673 FREESELF(c);
674 COPY(S_M_IJ(mat,0L,S_P_II(perm,0L)-1L),c);
675
676 for (i=1L;i<S_P_LI(perm);i++)
677 {
678 MULT_APPLY(S_M_IJ(mat,i,S_P_II(perm,i)-1L),c);
679 if (NULLP(c)) break;
680 }
681
682 CTO(ANYTYPE,"det270588(e)",c);
683 ENDR("det270588");
684 }
685
686
687
688 #ifdef CHARTRUE
det_mat_imm(mat,erg)689 INT det_mat_imm(mat,erg) OP mat,erg;
690 {
691 return det_imm_matrix(mat,erg);
692 }
693
det_imm_matrix(mat,b)694 INT det_imm_matrix(mat,b) OP mat,b;
695 /* AK 270588 neue version mit immanente */
696 /* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 180691 V1.2 */
697 /* AK 210891 V1.3 */
698 /* AK 260603 V2.0 */
699 {
700 OP part;
701 INT erg = OK; /* AK 301292 */
702 CTTO(MATRIX,INTEGERMATRIX,"det_imm_matrix(1)",mat);
703 SYMCHECK(S_M_HI(mat)!=S_M_LI(mat),"det_imm_matrix:not quadratic matrix");
704 CE2(mat,b,det_imm_matrix);
705
706
707 if (S_M_HI(mat) == 1L) /* AK 301292 */
708 {
709 COPY(S_M_IJ(mat,0L,0L),b);
710 }
711 else {
712 part = CALLOCOBJECT();
713 erg += last_partition(S_M_H(mat),part);
714 /* = 1,1,1,1,1,..,1 */
715 erg += immanente_matrix(mat,part,b);
716 /* der zugehoerige Charakter ist das signum */
717 FREEALL(part);
718 }
719 ENDR("det_imm_matrix");
720 }
721
co_050995(a)722 static INT co_050995(a) OP a;
723 {
724 INT i;
725 for (i=0; i< S_P_LI(a) ; i+= 2)
726 {
727 if (S_P_II(a,i) > S_P_II(a,i+1)) return FALSE;
728 if (i > 0)
729 if (S_P_II(a,i) < S_P_II(a,i-2)) return FALSE;
730 }
731 return TRUE;
732 }
733
pfaffian_matrix(mat,res)734 INT pfaffian_matrix(mat,res) OP mat,res;
735 /* berechnet pfaffian */
736 /* AK 050995 */
737 {
738 OP perm,zwerg,zzerg;
739 INT erg = OK;
740 CTO(MATRIX,"pfaffian_matrix(1)",mat);
741
742 SYMCHECK(not quadraticp(mat),"pfaffian:not quadratic matrix");
743 SYMCHECK(not evenp(S_M_H(mat)),"pfaffian:size of matrix not even");
744
745 perm = callocobject();
746 zwerg = callocobject();
747 zzerg = callocobject();
748 erg += first_permutation(S_M_H(mat),perm);
749
750 erg += det050995(mat,perm,res);
751 erg += signum(perm,zwerg);
752 erg += mult_apply(zwerg,res);
753
754 while (next(perm,perm))
755 {
756 if (co_050995(perm) == TRUE)
757 {
758 erg += det050995(mat,perm,zwerg);
759 erg += signum(perm,zzerg);
760 erg += mult_apply(zwerg,zzerg);
761 erg += add_apply(zzerg,res);
762 }
763 };
764 FREEALL3(perm,zwerg,zzerg);
765 ENDR("pfaffian_matrix");
766 }
767
immanente_matrix(mat,part,res)768 INT immanente_matrix(mat,part,res) OP mat,part,res;
769 /* berechnet immanente */
770 /* AK270588 */ /* AK 060789 V1.0 */ /* AK 090790 V1.1 */ /* AK 180691 V1.2 */
771 /* AK 210891 V1.3 */
772 /* AK 161098 V2.0 */
773 /* AK 130603 */
774 {
775 OP perm,nextperm,zwerg,zzerg;
776 INT i,erg = OK;
777
778 CTTO(MATRIX,INTEGERMATRIX,"immanente_matrix(1)",mat);
779 CTO(PARTITION,"immanente_matrix(2)",part);
780 SYMCHECK(S_M_HI(mat) != S_M_LI(mat),"immanente_matrix:not quadratic matrix");
781 PARTITION_WEIGHT(part,i); /* AK 130603 */
782 SYMCHECK(i != S_M_HI(mat),"immanente_matrix:wrong weight of partition");
783
784 CE3(mat,part,res,immanente_matrix);
785
786 perm = CALLOCOBJECT();
787 zwerg = CALLOCOBJECT();
788 zzerg = CALLOCOBJECT();
789 nextperm = CALLOCOBJECT();
790 erg += first_permutation(S_M_H(mat),perm);
791
792 erg += det270588(mat,perm,res);
793 erg += charvalue(part,perm,zwerg,NULL);
794 MULT_APPLY(zwerg,res);
795
796 while (next_apply(perm))
797 {
798 erg += det270588(mat,perm,zwerg);
799 erg += charvalue(part,perm,zzerg,NULL);
800 MULT_APPLY(zzerg,zwerg);
801 ADD_APPLY(zwerg,res);
802 };
803
804 FREEALL4(zzerg,zwerg,perm,nextperm);
805 ENDR("immanente_matrix");
806 }
807 #endif /* CHARTRUE */
808 #endif /* PERMTRUE */
809
810
811
inc_matrix(a)812 INT inc_matrix(a) OP a;
813 /* 250488 */ /* AK 060789 V1.0 *//* AK 130790 V1.1 */ /* AK 180691 V1.2 */
814 /* AK 210891 V1.3 */
815 /* AK 240804 V3.0 */
816 {
817 INT erg = OK;
818 CTTO(MATRIX,INTEGERMATRIX,"inc_matrix(1)",a);
819 {
820 OP l,h;
821 OP b; /* die neue matrix */
822 INT i,j;
823
824
825 CALLOCOBJECT3(l,h,b);
826
827 COPY_INTEGER(S_M_H(a),h); INC_INTEGER(h);
828 COPY_INTEGER(S_M_L(a),l); INC_INTEGER(l);
829 b_lh_m(l,h,b);C_O_K(b,S_O_K(a));
830
831 for (i=0L;i<S_M_HI(a);i++) for (j=0L;j<S_M_LI(a);j++)
832 *(S_M_IJ(b,i,j)) = *(S_M_IJ(a,i,j));
833 for (i=0L;i<S_M_HI(b);i++) C_O_K(S_M_IJ(b,i,S_M_LI(a)),EMPTY);
834 for (j=0L;j<S_M_LI(b);j++) C_O_K(S_M_IJ(b,S_M_HI(a),j),EMPTY);
835
836 SYM_free(S_M_S(a));
837 FREEALL2(S_M_H(a),S_M_L(a));
838 SYM_free(S_O_S(a).ob_matrix);
839
840 *a = *b;
841 C_O_K(b,EMPTY); FREEALL(b);
842 }
843 ENDR("inc_matrix");
844 }
845
inc_matrix_row_co(a,k)846 INT inc_matrix_row_co(a,k) OP a; INT k;
847 /* 080304 AK */ /* increase matrix by k empty row */
848 /* AK 240904 V3.0 */
849 {
850 INT erg = OK;
851 CTTO(MATRIX,INTEGERMATRIX,"inc_matrix_row_co(1)",a);
852 SYMCHECK(k<=0,"inc_matrix_row_co:parameter <0");
853 {
854 OP b; /* the new matrix */
855 OP l,h;
856 INT i,j;
857
858
859 CALLOCOBJECT3(l,h,b);
860
861 M_I_I(S_M_HI(a)+k,h);
862 M_I_I(S_M_LI(a),l);
863 b_lh_m(l,h,b);C_O_K(b,S_O_K(a));
864
865 for (i=0L;i<S_M_HI(a);i++) for (j=0L;j<S_M_LI(a);j++)
866 *(S_M_IJ(b,i,j)) = *(S_M_IJ(a,i,j));
867 for (i=0L;i<k;i++) for (j=0L;j<S_M_LI(a);j++)
868 C_O_K(S_M_IJ(b,i+S_M_HI(a),j),EMPTY);
869
870 SYM_free(S_M_S(a));
871 FREEALL2(S_M_H(a),S_M_L(a));
872 SYM_free(S_O_S(a).ob_matrix);
873
874 *a = *b;
875 C_O_K(b,EMPTY); FREEALL(b);
876 }
877 ENDR("inc_matrix_row_co");
878 }
879
inc_matrix_column_co(a,k)880 INT inc_matrix_column_co(a,k) OP a; INT k;
881 /* 080304 AK */
882 /* increase matrix by k empty columns at the end */
883 /* AK 240904 V3.0 */
884 /* AK 180607 V3.1 */
885 {
886 INT erg = OK;
887 CTTO(MATRIX,INTEGERMATRIX,"inc_matrix_column_co(1)",a);
888 SYMCHECK(k<=0,"inc_matrix_column_co:parameter <0");
889 {
890 OP l,h;
891 OP b; /* die neue matrix */
892 INT i,j;
893
894
895 CALLOCOBJECT3(l,h,b);
896
897 M_I_I(S_M_LI(a)+k,l);
898 M_I_I(S_M_HI(a),h);
899 erg += b_lh_m(l,h,b);
900 C_O_K(b,S_O_K(a));
901
902 for (i=0L;i<S_M_HI(a);i++) for (j=0L;j<S_M_LI(a);j++)
903 *(S_M_IJ(b,i,j)) = *(S_M_IJ(a,i,j));
904 for (i=0L;i<S_M_HI(a);i++) for (j=0L;j<k;j++)
905 C_O_K(S_M_IJ(b,i,j+S_M_LI(a)),EMPTY);
906
907 SYM_free(S_M_S(a));
908 FREEALL2(S_M_H(a),S_M_L(a));
909 SYM_free(S_O_S(a).ob_matrix);
910
911 *a = *b;
912 C_O_K(b,EMPTY); FREEALL(b);
913 }
914 ENDR("inc_matrix_column_co");
915 }
916
917
singularp(a)918 INT singularp(a) OP a;
919 /* AK 020304 */
920 /* return TRUE
921 if matrix singualar i.e row rank != number of rows */
922 /* AK 021106 V3.1 */
923 {
924 INT erg = OK;
925 CTTO(MATRIX,INTEGERMATRIX,"singularp(1)",a);
926 SYMCHECK(not quadraticp(a),"singularp:not quadratic");
927 {
928 OP b;INT res;
929 b=CALLOCOBJECT();
930 erg += rank(a,b);
931 if (EQ(b,S_M_H(a))) res = FALSE;
932 else res = TRUE;
933 FREEALL(b);
934 return res;
935 }
936 ENDR("singularp");
937 }
938
invers_matrix(a,b)939 INT invers_matrix(a,b) OP a,b;
940 /* AK 290388 nach stoer (dietmar) */
941 /* umgewandelt aus pascal */
942 /* AK 060789 V1.0 */ /* AK 181289 V1.1 */ /* AK 150591 V1.2 */
943 /* AK 220791 V1.3 */ /* AK 081204 V3.0 */
944 /* AK 021206 V3.1 */
945 {
946 INT erg=OK;
947 CTO(MATRIX,"invers_matrix(1)",a);
948 CTO(EMPTY,"invers_matrix(2)",b);
949 SYMCHECK(not quadraticp(a),"invers_matrix:not quadratic");
950 {
951 INT i,j,k,r;
952 /* r ist die selectierte spalte */
953 /* r = 0 ... n */
954 INT n=S_M_LI(a)-1L;
955 INT singulaer = FALSE;
956 OP p,hr,hs,hv;
957 CALLOCOBJECT4(p,hr,hs,hv);
958
959
960 erg += m_il_v(n+1L,p);
961
962 for(j=0L;j<=n;j++) M_I_I(j,S_V_I(p,j));
963
964 j= -1L;
965 CLEVER_COPY(a,b);
966 while ((j++ <n) && (! singulaer))
967 {
968 /*pivotsuche*/
969 for(r=j;r<=n;r++)
970 {
971 if (S_O_K(S_M_IJ(b,r,j))==GALOISRING) // non field
972 {
973 if (unitp_galois(S_M_IJ(b,r,j))) goto im290388;
974 }
975 else if (not NULLP(S_M_IJ(b,r,j))) goto im290388;
976 }
977
978 im290388:
979 if (r == n+1L)/* nur nullen in der spalte j */ singulaer = TRUE;
980 else {
981 /*zeilentausch*/
982 if (r>j){
983 for (k=0L;k<=n;k++)
984 SWAP(S_M_IJ(b,j,k),S_M_IJ(b,r,k));
985 SWAP(S_V_I(p,j),S_V_I(p,r));
986 };
987 /*transformation*/
988 FREESELF(hr);
989 INVERS(S_M_IJ(b,j,j),hr);
990 for (i=0L;i<=n;i++) MULT_APPLY(hr,S_M_IJ(b,i,j));
991 CLEVER_COPY(hr,S_M_IJ(b,j,j));
992 ADDINVERS_APPLY(hr);
993 for (k=0L;k<=n;k++)
994 {
995 if (k==j) k++; /* spalte j nicht anwenden */
996 if (k>n) break;
997 for (i=0L;i<=n;i++)
998 {
999 if (i==j) i++; /* auf zeile j nicht anwenden */
1000 if (i>n) break;
1001 MULT(S_M_IJ(b,i,j),S_M_IJ(b,j,k),hs);
1002 ADDINVERS_APPLY(hs);
1003 ADD_APPLY(hs,S_M_IJ(b,i,k));
1004 FREESELF(hs);
1005 };
1006 MULT_APPLY(hr,S_M_IJ(b,j,k));
1007 };
1008 }; /* end else */
1009 }; /* end while */
1010 // if (erg != OK) goto endr_ende;
1011 FREEALL2(hr,hs);
1012
1013 erg+=m_il_v(n+1L,hv);
1014 if (not singulaer)
1015 /*spaltentausch*/
1016 for (i=0L;i<=n;i++)
1017 {
1018 OP z;
1019 z = S_M_IJ(b,i,0);
1020 for (k=0;k<=n;k++,z++) CLEVER_COPY(z, S_V_I(hv,S_V_II(p,k)));
1021 z = S_M_IJ(b,i,0);
1022 for (k=0L;k<=n;k++,z++) CLEVER_COPY(S_V_I(hv,k), z);
1023 }
1024
1025 FREEALL2(p,hv);
1026 if (singulaer) {
1027 freeself(b);
1028 error("invers_matrix: singulary");
1029 return(SINGULAER);
1030 };
1031 }
1032 ENDR("invers_matrix");
1033
1034 }
1035
1036
transpose_matrix(a,b)1037 INT transpose_matrix(a,b) OP a,b;
1038 /* AK 280388 */ /* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 210891 V1.3 */
1039 /* AK 081204 V3.0 */
1040 /* AK 021106 V3.1 */
1041 {
1042 INT erg = OK;
1043 CTO(MATRIX,"transpose_matrix(1)",a);
1044 CTO(EMPTY,"transpose_matrix(2)",b);
1045 {
1046
1047 INT i,j;
1048 erg += m_ilih_m(S_M_HI(a),S_M_LI(a),b);
1049 C_O_K(b,S_O_K(a));
1050
1051 for (i=0;i<S_M_HI(b);i++)
1052 for (j=0;j<S_M_LI(b);j++)
1053 COPY(S_M_IJ(a,j,i),S_M_IJ(b,i,j));
1054 }
1055 ENDR("transpose_matrix");
1056 }
1057
transpose_second_matrix(a,b)1058 INT transpose_second_matrix(a,b) OP a,b;
1059 /* AK 280388 */ /* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 210891 V1.3 */
1060 /* not on the main diagonal */
1061 /* AK 021106 V3.1 */
1062 {
1063 INT i,j;
1064 INT erg = OK;
1065 CE2(a,b,transpose_second_matrix);
1066
1067 erg += m_ilih_m(S_M_HI(a),S_M_LI(a),b);
1068 C_O_K(b,S_O_K(a));
1069
1070 for (i=0;i<S_M_HI(b);i++)
1071 for (j=0;j<S_M_LI(b);j++)
1072 erg += copy(S_M_IJ(a,j,i),S_M_IJ(b,S_M_HI(a)-1-i,j));
1073
1074 ENDR("transpose_second_matrix");
1075 }
1076
1077
comp_kranztafel(a,b)1078 INT comp_kranztafel(a,b) OP a, b;
1079 /* AK 280390 V1.1 */ /* AK 210891 V1.3 */
1080 {
1081 INT i,j,res;
1082 OP x,y;
1083
1084 x = S_M_S(a);
1085 y = S_M_S(b);
1086 for (i=0L;i<S_M_HI(a);i++)
1087 {
1088 if (i >= S_M_HI(b)) return(1L);
1089 else {
1090 for (j=0;j<S_M_LI(a);j++)
1091 {
1092 if (j>=S_M_LI(b)) return(1L);
1093 else {
1094 res=COMP_INTEGER_INTEGER(x,y);
1095 if (res != 0) return(res);
1096 x++;y++;
1097 };
1098 }
1099 }
1100 }
1101 if ( S_M_HI(b) > S_M_HI(a) ) return -1;
1102 if ( S_M_LI(b) > S_M_LI(a) ) return -1;
1103 return 0;
1104 }
1105
comp_integermatrix(a,b)1106 INT comp_integermatrix(a,b) OP a, b;
1107 /* AK 150802 */
1108 /* in case of equal dimensions lexicographic */
1109 {
1110 INT erg = OK;
1111 CTO(INTEGERMATRIX,"comp_integermatrix(1)",a);
1112 CTO(INTEGERMATRIX,"comp_integermatrix(2)",b);
1113 {
1114 INT i,j,res;
1115 OP x,y;
1116
1117 x = S_M_S(a);
1118 y = S_M_S(b);
1119 for (i=0L;i<S_M_HI(a);i++)
1120 {
1121 if (i >= S_M_HI(b)) return(1L);
1122 else {
1123 for (j=(INT)0;j<S_M_LI(a);j++)
1124 {
1125 if (j>=S_M_LI(b)) return(1L);
1126 else {
1127 res=COMP_INTEGER_INTEGER(x,y);
1128 if (res != 0) return(res);
1129 x++;y++;
1130 };
1131 }
1132 }
1133 }
1134 if ( S_M_HI(b) > S_M_HI(a) ) return -1;
1135 if ( S_M_LI(b) > S_M_LI(a) ) return -1;
1136 return 0;
1137 }
1138 ENDR("comp_integermatrix");
1139 }
1140
1141
1142
comp_matrix(a,b)1143 INT comp_matrix(a,b) OP a, b;
1144 /* AK 060789 V1.0 */ /* AK 070290 V1.1 */ /* AK 210891 V1.3 */
1145 {
1146 INT i,j,res;
1147 OP x,y;
1148
1149 x = S_M_S(a);
1150 y = S_M_S(b);
1151 for (i=(INT)0;i<S_M_HI(a);i++)
1152 {
1153 if (i >= S_M_HI(b)) return(1L);
1154 else {
1155 for (j=(INT)0;j<S_M_LI(a);j++)
1156 {
1157 if (j>=S_M_LI(b)) return(1L);
1158 else {
1159 res=COMP(x,y);
1160 if (res != 0) return(res);
1161 x++;
1162 y++;
1163 };
1164 }
1165 }
1166 }
1167 if ( S_M_HI(b) > S_M_HI(a) ) return(-1L); /* AK 170790 */
1168 if ( S_M_LI(b) > S_M_LI(a) ) return(-1L); /* AK 241195 */
1169 return((INT)0); /* matrizen sind gleich */
1170 }
1171
1172
1173
add_apply_matrix_matrix(a,b)1174 INT add_apply_matrix_matrix(a,b) OP a,b;
1175 /* AK 220390 V1.1 */ /* AK 090891 V1.3 */
1176 /* AK 210891 V1.3 */
1177 /* b:= b += a */
1178 {
1179 OP c,d;
1180 INT erg = OK;
1181 CTTTO(MATRIX,INTEGERMATRIX,KRANZTYPUS,"add_apply_matrix_matrix(1)",a);
1182 CTTTO(MATRIX,INTEGERMATRIX,KRANZTYPUS,"add_apply_matrix_matrix(2)",b);
1183
1184 C_M_HASH(b,-1);
1185 if (
1186 (S_M_HI(a) == S_M_HI(b))&&
1187 (S_M_LI(a) == S_M_LI(b))
1188 )
1189 {
1190 INT i;
1191 i = S_M_HI(a)*S_M_LI(a);
1192 c = S_M_S(a);
1193 d=S_M_S(b);
1194 while (i-- > 0)
1195 {
1196 ADD_APPLY(c,d);
1197 c++;
1198 d++;
1199 }
1200 }
1201 else if (
1202 (S_M_HI(a) < S_M_HI(b)) &&
1203 (S_M_LI(a) < S_M_LI(b))
1204 )
1205 {
1206 INT i,j;
1207 for (i=0;i<S_M_HI(a);i++)
1208 for (j=0;j<S_M_LI(a);j++)
1209 ADD_APPLY(S_M_IJ(a,i,j),S_M_IJ(b,i,j));
1210 }
1211 else {
1212 c = CALLOCOBJECT();
1213 *c = *b;
1214 C_O_K(b,EMPTY);
1215 erg += add_matrix_matrix(a,c,b);
1216 FREEALL(c);
1217 }
1218 ENDR("add_apply_matrix_matrix");
1219 }
1220
1221
1222
add_apply_matrix(a,b)1223 INT add_apply_matrix(a,b) OP a,b;
1224 /* AK 220390 V1.1 */
1225 /* AK 210891 V1.3 */
1226 {
1227 INT erg = OK;
1228 CTTTO(MATRIX,INTEGERMATRIX,KRANZTYPUS,"add_apply_matrix(1)",a);
1229 EOP("add_apply_matrix(2)",b);
1230
1231
1232 switch (S_O_K(b)){
1233 case KRANZTYPUS:
1234 case INTEGERMATRIX:
1235 case MATRIX:
1236 erg += add_apply_matrix_matrix(a,b);
1237 break;
1238 default:
1239 WTO("add_apply_matrix(2)",b);
1240 break;
1241 }
1242 ENDR("add_apply_matrix");
1243 }
1244
1245
1246
add_matrix(a,b,c)1247 INT add_matrix(a,b,c) OP a,b,c;
1248 {
1249 INT erg=OK;
1250 if (not MATRIXP(a))
1251 { erg += WTO("add_matrix",a);goto endr_ende; }
1252 if (MATRIXP(b))
1253 erg += add_matrix_matrix(a,b,c);
1254 else
1255 erg += WTO("add_matrix",b);
1256 ENDR("add_matrix");
1257 }
1258
add_matrix_matrix(a,b,ergeb)1259 INT add_matrix_matrix(a,b,ergeb) OP a,b,ergeb;
1260 /* AK 041186 */ /* AK 171186 */ /* AK 060789 V1.0 */ /* AK 081289 V1.1 */
1261 /* AK 090891 V1.3 */
1262 /* a and b may have different sizes */
1263 {
1264 INT i,j;
1265 OP len,height;
1266 OP z;
1267 INT erg = OK;
1268
1269 CTTO(INTEGERMATRIX,MATRIX,"add_matrix_matrix(1)",a);
1270 CTTO(INTEGERMATRIX,MATRIX,"add_matrix_matrix(2)",b);
1271 CTO(EMPTY,"add_matrix_matrix(3)",ergeb);
1272
1273 len = CALLOCOBJECT(),
1274 height = CALLOCOBJECT();
1275 if (S_M_LI(a) >=S_M_LI(b))
1276 COPY_INTEGER(S_M_L(a),len);
1277 else COPY_INTEGER(S_M_L(b),len);
1278 if (S_M_HI(a) >=S_M_HI(b))
1279 COPY_INTEGER(S_M_H(a),height);
1280 else COPY_INTEGER(S_M_H(b),height);
1281
1282 erg += b_lh_m(len,height,ergeb);
1283 C_O_K(ergeb,S_O_K(a));
1284
1285 z = S_M_S(ergeb);
1286 for (i=0;i<S_M_HI(ergeb);i++)
1287 for (j=0;j<S_M_LI(ergeb);j++,z++)
1288 {
1289 if (
1290 (i<S_M_HI(a))&&(i<S_M_HI(b))&&
1291 (j<S_M_LI(a))&&(j<S_M_LI(b))
1292 )
1293
1294 {
1295 ADD(S_M_IJ(a,i,j), S_M_IJ(b,i,j), z);
1296 }
1297
1298 else if ( (i<S_M_HI(a))&&(j<S_M_LI(a)))
1299 COPY(S_M_IJ(a,i,j),z);
1300 else if ( (i<S_M_HI(b))&&(j<S_M_LI(b)))
1301 COPY(S_M_IJ(b,i,j),z);
1302 else { /* for the multiplication of matrix labelled polynomials */
1303 if (S_O_K(S_M_IJ(a,0,0)) == INTEGER)
1304 M_I_I(0,z);
1305 }
1306 }
1307 ENDR("add_matrix_matrix");
1308 }
1309
1310
1311
copy_integermatrix(a,b)1312 INT copy_integermatrix(a,b) OP a,b;
1313 {
1314 INT erg = OK;
1315 CTO(INTEGERMATRIX,"copy_integermatrix(1)",a);
1316 CTO(EMPTY,"copy_integermatrix(2)",b);
1317
1318 erg += m_ilih_m(S_M_LI(a),S_M_HI(a),b);
1319 C_O_K(b,S_O_K(a));
1320 C_M_HASH(b,S_M_HASH(a));
1321 memcpy((char *) S_M_S(b),(char *) S_M_S(a),
1322 S_M_LI(a)*S_M_HI(a)*sizeof(struct object));
1323 ENDR("copy_integermatrix");
1324 }
1325
freeself_integermatrix(a)1326 INT freeself_integermatrix(a) OP a;
1327 /* AK 281098 V2.0 */
1328 {
1329 INT erg = OK;
1330 CTO(INTEGERMATRIX,"freeself_integermatrix(1)",a);
1331 {
1332 OBJECTSELF d;
1333 d=S_O_S(a);
1334 SYM_free(S_M_S(a));
1335 FREEALL(S_M_L(a));
1336 FREEALL(S_M_H(a));
1337 SYM_free(d.ob_matrix);
1338 C_O_K(a,EMPTY);
1339 }
1340 ENDR("freeself_integermatrix");
1341 }
1342
1343
1344
1345
copy_kranztypus(a,b)1346 INT copy_kranztypus(a,b) OP a,b;
1347 /* AK 270390 V1.1 */ /* AK 210891 V1.3 */
1348 /* AK 040392 cast to char * */
1349 {
1350 m_ilih_m(S_M_LI(a),S_M_HI(a),b);
1351 C_O_K(b,S_O_K(a));
1352 memcpy((char *) S_M_S(b),(char *) S_M_S(a),
1353 S_M_LI(a)*S_M_HI(a)*sizeof(struct object));
1354 return(OK);
1355 }
1356
1357
1358
copy_matrix(von,b)1359 INT copy_matrix(von,b) OP von , b;
1360 /* AK 051186 */ /* AK 021286 */ /* AK 060789 V1.0 */ /* AK 071289 V1.1 */
1361 /* AK 210891 V1.3 */
1362 {
1363 INT k;
1364 OP z,w;
1365 INT erg = OK;
1366 CTTO(INTEGERMATRIX,MATRIX,"copy_matrix(1)",von);
1367 CTO(EMPTY,"copy_matrix(2)",b);
1368
1369 erg += m_ilih_m(S_M_LI(von),S_M_HI(von),b);
1370 C_O_K(b,S_O_K(von));
1371 C_M_HASH(b,S_M_HASH(von));
1372
1373
1374 z = S_M_IJ(von,S_M_HI(von)-1L,S_M_LI(von)-1L);
1375 w = S_M_IJ(b,S_M_HI(von)-1L,S_M_LI(von)-1L);
1376 k = S_M_HI(von) * S_M_LI(von);
1377 for (;k>0;k--,z--,w--)
1378 COPY(z,w);
1379
1380 ENDR("copy_matrix");
1381 }
1382
1383
1384
freeself_kranztypus(a)1385 INT freeself_kranztypus(a) OP a;
1386 /* AK 270390 V1.1 */ /* AK 210891 V1.3 */
1387 /* AK 281098 V2.0 */
1388 {
1389 INT erg = OK;
1390 CTO(KRANZTYPUS,"freeself_kranztypus(1)",a);
1391 {
1392 OBJECTSELF d;
1393 d=S_O_S(a);
1394 SYM_free(S_M_S(a));
1395 FREEALL(S_M_L(a));
1396 FREEALL(S_M_H(a));
1397 SYM_free(d.ob_matrix);
1398 C_O_K(a,EMPTY);
1399 }
1400 ENDR("freeself_kranztypus");
1401 }
1402
1403
1404
freeself_matrix(matrix)1405 INT freeself_matrix(matrix) OP matrix;
1406 /* AK 060789 V1.0 */ /* AK 071289 V1.1 */ /* AK 100691 V1.2 */
1407 /* AK 160891 V1.3 */
1408 {
1409 INT k;
1410 OBJECTSELF d;
1411 OP z;
1412 INT erg = OK;
1413 CTO(MATRIX,"freeself_matrix(1)",matrix);
1414
1415 d=S_O_S(matrix);
1416
1417 z = S_M_IJ(matrix,S_M_HI(matrix)-1L,S_M_LI(matrix)-1L);
1418 k = S_M_HI(matrix) * S_M_LI(matrix);
1419 for (;k>(INT)0;k--,z--)
1420 if (S_O_K(z) == INTEGER) ;
1421 else if (EMPTYP(z));
1422 else
1423 erg += freeself(z);
1424
1425 SYM_free(S_M_S(matrix));
1426 erg += freeall(S_M_L(matrix));
1427 erg += freeall(S_M_H(matrix));
1428 SYM_free(d.ob_matrix);
1429 C_O_K(matrix,EMPTY);
1430 ENDR("freeself_matrix");
1431 }
1432
1433
1434
callocmatrix()1435 static struct matrix * callocmatrix()
1436 /* AK 060789 V1.0 */ /* AK 220390 V1.1 */ /* AK 160891 V1.3 */
1437 {
1438 struct matrix *ergebnis;
1439
1440 ergebnis = (struct matrix *) SYM_calloc((int)1,sizeof(struct matrix));
1441 if (ergebnis == NULL)
1442 no_memory();
1443 return(ergebnis);
1444 }
1445
1446
1447
scan_integermatrix(ergebnis)1448 INT scan_integermatrix(ergebnis) OP ergebnis;
1449 /* AK 141293 */
1450 {
1451 return scan_matrix_co(ergebnis, INTEGER);
1452 }
1453
scan_matrix(ergebnis)1454 INT scan_matrix(ergebnis) OP ergebnis;
1455 /* AK 141293 */
1456 {
1457 return scan_matrix_co(ergebnis, EMPTY);
1458 }
1459
1460
scan_skewsymmetric_matrix(ergebnis)1461 INT scan_skewsymmetric_matrix(ergebnis) OP ergebnis;
1462 /* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 080891 V1.3 */
1463 {
1464 OP len, height;
1465 INT i,j;
1466 char a[20]; /* AK 080891 */
1467 OBJECTKIND kind;
1468
1469 len = callocobject();
1470 height = callocobject();
1471 aaa:
1472 printeingabe("height of skew symmetric matrix");
1473 scan(INTEGER,height);
1474 copy(height,len);
1475 printeingabe("enter kind of matrix elements");
1476 kind=scanobjectkind();
1477 if (S_I_I(len) <= (INT)0) /* AK 170795 */
1478 {
1479 printeingabe("you entered wrong length (<=0), do it again");
1480 goto aaa;
1481 }
1482 if (S_I_I(height) <= (INT)0) /* AK 170795 */
1483 {
1484 printeingabe("you entered wrong height (<=0), do it again");
1485 goto aaa;
1486 }
1487 b_lh_m(len,height,ergebnis);
1488 for (i=0; i<S_I_I(height); i++)
1489 m_i_i(0L,S_M_IJ(ergebnis,i,i));
1490 for (i=0; i<S_I_I(height); i++)
1491 {
1492 sprintf(a,"row nr %ld \n",(i+1L)); /* AK 080891 */
1493 printeingabe(a); /* AK 080891 */
1494 for (j=i+1;j<S_I_I(len);j++)
1495 {
1496 scan(kind,S_M_IJ(ergebnis,i,j));
1497 addinvers(S_M_IJ(ergebnis,i,j), S_M_IJ(ergebnis,j,i));
1498 }
1499 };
1500 return(OK);
1501 }
1502
scan_matrix_co(ergebnis,kind)1503 static INT scan_matrix_co(ergebnis,kind) OP ergebnis; OBJECTKIND kind;
1504 /* AK 060789 V1.0 */ /* AK 310790 V1.1 */ /* AK 080891 V1.3 */
1505 {
1506 OP len, height;
1507 INT i,j;
1508 char a[20]; /* AK 080891 */
1509
1510 len = callocobject();
1511 height = callocobject();
1512 aaa:
1513 printeingabe("height of matrix");
1514 scan(INTEGER,height);
1515 printeingabe("length of matrix");
1516 scan(INTEGER,len);
1517 if (kind == EMPTY)
1518 {
1519 printeingabe("enter kind of matrix elements");
1520 kind=scanobjectkind();
1521 }
1522 if (S_I_I(len) <= (INT)0) /* AK 170795 */
1523 {
1524 printeingabe("you entered wrong length (<=0), do it again");
1525 goto aaa;
1526 }
1527 if (S_I_I(height) <= (INT)0) /* AK 170795 */
1528 {
1529 printeingabe("you entered wrong height (<=0), do it again");
1530 goto aaa;
1531 }
1532 b_lh_m(len,height,ergebnis);
1533 for (i=0; i<S_I_I(height); i++)
1534 {
1535 sprintf(a,"row nr %ld \n",(i+1L)); /* AK 080891 */
1536 printeingabe(a); /* AK 080891 */
1537 for (j=0;j<S_I_I(len);j++)
1538 scan(kind,S_M_IJ(ergebnis,i,j));
1539 };
1540 return(OK);
1541 }
1542
1543
1544
change_column_ij(a,i,j)1545 INT change_column_ij(a,i,j) OP a; INT i,j;
1546 /* AK 301288 vertauscht spalten i und j dabei kein copy */
1547 /* AK 060789 V1.0 */ /* AK 050390 V1.1 */ /* AK 160891 V1.3 */
1548 {
1549 INT k,erg = OK;
1550 CTO(MATRIX,"change_column_ij(1)",a);
1551 SYMCHECK(i<0,"change_column_ij:i<0");
1552 SYMCHECK(j<0,"change_column_ij:j<0");
1553 SYMCHECK(i>=S_M_LI(a),"change_column_ij:i too big ");
1554 SYMCHECK(j>=S_M_LI(a),"change_column_ij:j too big ");
1555 if (i==j) goto endr_ende; /* AK 190802 */
1556
1557 for (k=0; k<S_M_HI(a); k++)
1558 SWAP(S_M_IJ(a,k,i),S_M_IJ(a,k,j));
1559 ENDR("change_column_ij");
1560 }
1561
1562
change_row_ij(a,i,j)1563 INT change_row_ij(a,i,j) OP a; INT i,j;
1564 /* AK 301288 vertauscht zeilen i und j dabei kein copy */
1565 /* AK 060789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1566 {
1567 INT k,erg = OK;
1568 CTO(MATRIX,"change_row_ij(1)",a);
1569 SYMCHECK(i<0,"change_row_ij:i<0");
1570 SYMCHECK(j<0,"change_row_ij:j<0");
1571 SYMCHECK(i>=S_M_HI(a),"change_row_ij:i too big ");
1572 SYMCHECK(j>=S_M_HI(a),"change_row_ij:j too big ");
1573 if (i==j) goto endr_ende; /* AK 190802 */
1574
1575 for (k=(INT)0; k<S_M_LI(a); k++)
1576 SWAP(S_M_IJ(a,i,k),S_M_IJ(a,j,k));
1577 ENDR("change_row_ij");
1578 }
1579
1580
s_m_s(a)1581 OP s_m_s(a) OP a;
1582 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1583 { OBJECTSELF c; c = s_o_s(a); return(c.ob_matrix->m_self); }
1584
s_m_h(a)1585 OP s_m_h(a) OP a;
1586 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1587 { OBJECTSELF c; c = s_o_s(a); return(c.ob_matrix->m_height); }
1588
s_m_l(a)1589 OP s_m_l(a) OP a;
1590 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1591 { OBJECTSELF c; c = s_o_s(a); return(c.ob_matrix->m_length); }
1592
s_m_hash(a)1593 INT s_m_hash(a) OP a;
1594 /* AK 110703 */
1595 { OBJECTSELF c; c = s_o_s(a); return(c.ob_matrix->m_hash); }
1596
s_m_hi(a)1597 INT s_m_hi(a) OP a;
1598 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1599 { return(s_i_i(s_m_h(a))); }
1600
s_m_li(a)1601 INT s_m_li(a) OP a;
1602 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1603 { return(s_i_i(s_m_l(a))); }
1604
c_m_s(a,b)1605 INT c_m_s(a,b) OP a,b;
1606 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1607 { OBJECTSELF c; c = s_o_s(a); c.ob_matrix->m_self = b; return(OK); }
1608
1609
c_m_h(a,b)1610 INT c_m_h(a,b) OP a,b;
1611 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1612 { OBJECTSELF c; c = s_o_s(a); c.ob_matrix->m_height = b; return(OK); }
1613
c_m_hash(a,b)1614 INT c_m_hash(a,b) OP a; INT b;
1615 /* AK 110703 */
1616 { OBJECTSELF c; c = s_o_s(a); c.ob_matrix->m_hash = b; return(OK); }
1617
1618
c_m_l(a,b)1619 INT c_m_l(a,b) OP a,b;
1620 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1621 { OBJECTSELF c; c = s_o_s(a); c.ob_matrix->m_length = b; return(OK); }
1622
s_m_ij(a,i,j)1623 OP s_m_ij(a,i,j) OP a; INT i,j;
1624 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1625 {
1626 if (not MATRIXP(a))
1627 {
1628 printobjectkind(a);
1629 error("s_m_ij:no matrix object");
1630 }
1631 if (i < (INT)0)
1632 {
1633 debugprint(a);
1634 fprintf(stderr, "index = %" PRIINT "\n" ,i);
1635 error("s_m_ij:row index too small");
1636 }
1637 if (i >= s_m_hi(a))
1638 {
1639 debugprint(a);
1640 fprintf(stderr, "index = %" PRIINT "\n" ,i);
1641 error("s_m_ij:row index too big");
1642 }
1643 if (j >= s_m_li(a))
1644 {
1645 debugprint(a);
1646 fprintf(stderr, "index = %" PRIINT "\n" ,j);
1647 error("s_m_ij:column index too big");
1648 }
1649 if (j < (INT)0)
1650 {
1651 debugprint(a);
1652 fprintf(stderr, "index = %" PRIINT "\n" ,j);
1653 error("s_m_ij:column index too small");
1654 }
1655 return(s_m_s(a) + (s_m_li(a)*i+j) );
1656 }
1657
s_m_iji(a,i,j)1658 INT s_m_iji(a,i,j) OP a; INT i,j;
1659 /* AK 070789 V1.0 */ /* AK 181289 V1.1 */ /* AK 160891 V1.3 */
1660 { return(s_i_i(s_m_ij(a,i,j))); }
1661
1662
1663
fprint_matrix(f,obj)1664 INT fprint_matrix(f,obj) FILE *f; OP obj;
1665 /* AK 211186 */ /* AK 070789 V1.0 */ /* AK 181289 V1.1 */
1666 /* AK 210891 V1.3 */
1667 {
1668 INT i,j;
1669 for (i=0;i<S_M_HI(obj);i++)
1670 {
1671 fprintf(f,"\n[");
1672 if (f == stdout) zeilenposition=0;
1673 for (j=0;j<S_M_LI(obj);j++)
1674 {
1675 fprint(f,S_M_IJ(obj,i,j));
1676 if (j+1 < S_M_LI(obj))
1677 {
1678 fprintf(f,":");
1679 if (f == stdout) zeilenposition++;
1680 }
1681 if ((f == stdout)&&(zeilenposition>70L))
1682 {fprintf(stdout,"\n");zeilenposition = 0L;}
1683 };
1684 fprintf(f,"]");
1685 };
1686 fprintf(f,"\n");
1687 if (f == stdout) zeilenposition=0L;
1688 return(OK);
1689 }
1690
tex_matrix(obj)1691 INT tex_matrix(obj) OP obj;
1692 {
1693 return tex_matrix_co(obj,tex);
1694 }
1695
tex_matrix_co(obj,f)1696 INT tex_matrix_co(obj,f) OP obj; INT (*f)();
1697 /* AK 150988 */ /* AK 310790 V1.1 */
1698 /* AK 070291 V1.2 texout for output */ /* AK 210891 V1.3 */
1699 {
1700 INT i,j;
1701 INT ts = texmath_yn; /* AK 190892 */
1702 INT erg = OK;
1703 CTO(MATRIX,"tex_matrix_co(1)",obj);
1704
1705 fprintf(texout,"\n");
1706 if (texmath_yn == 0L) /* d.h. not in math mode */
1707 {
1708 fprintf(texout,"$");
1709 texmath_yn=1L;
1710 }
1711 fprintf(texout,"\\matrix { \n");
1712 texposition = 0L;
1713 for (i=0;i<S_M_HI(obj);i++)
1714 {
1715 for (j=0L;j<S_M_LI(obj);j++)
1716 {
1717 (*f)(S_M_IJ(obj,i,j));
1718 fprintf(texout," & ");
1719 texposition += 3L;
1720 }
1721 fprintf(texout," \\cr\n");
1722 texposition=0L;
1723 };
1724 fprintf(texout," }");
1725 if (ts == 0L)
1726 {
1727 fprintf(texout,"$");
1728 texmath_yn=0L;
1729 }
1730 fprintf(texout," \n");
1731 texposition=0L;
1732 ENDR("tex_matrix_co");
1733 }
1734
1735
mult_scalar_matrix(scalar,mat,res)1736 INT mult_scalar_matrix(scalar,mat,res) OP scalar , mat,res;
1737 /* AK 310588 */
1738 /* AK 010889 V1.0 */ /* AK 081289 V1.1 */ /* AK 210891 V1.3 */
1739 {
1740 INT i,j,erg = OK;
1741 OP height,len;
1742 SYMCHECK( S_M_S(mat) == NULL,"mult_scalar_matrix:self == NULL");
1743
1744 height = callocobject();
1745 len = callocobject();
1746
1747 COPY_INTEGER(S_M_L(mat), len);
1748 COPY_INTEGER(S_M_H(mat), height);
1749 erg += b_lh_m(len,height,res);
1750 for (i=0L;i<S_I_I(height);i++)
1751 for(j=0L;j<S_I_I(len);j++)
1752 MULT(scalar,S_M_IJ(mat,i,j),S_M_IJ(res,i,j));
1753 ENDR("mult_scalar_matrix");
1754 }
1755
mult_apply_scalar_matrix(a,b)1756 INT mult_apply_scalar_matrix(a,b) OP a,b;
1757 /* AK 150290 V1.1 */ /* AK 210891 V1.3 */
1758 {
1759 OP z = S_M_S(b);
1760 INT grenze = S_M_LI(b)*S_M_HI(b);
1761 INT erg = OK;
1762 INT i;
1763 CTO(MATRIX,"mult_apply_scalar_matrix(2)",b);
1764 C_M_HASH(b,-1);
1765 for (i=0L;i<grenze;i++,z++)
1766 erg += mult_apply(a,z);
1767 ENDR("mult_apply_scalar_matrix");
1768 }
1769
mult_apply_matrix_matrix(a,b)1770 INT mult_apply_matrix_matrix(a,b) OP a,b;
1771 /* AK 131190 V1.1 */ /* b = a * b */
1772 /* AK 210891 V1.3 */
1773 {
1774 OP c = callocobject();
1775 INT erg = OK; /* AK 200192 */
1776 CTO(MATRIX,"mult_apply_matrix_matrix(1)",a);
1777 CTO(MATRIX,"mult_apply_matrix_matrix(2)",b);
1778 *c = *b;
1779 C_O_K(b,EMPTY);
1780 erg += mult_matrix_matrix(a,c,b);
1781 erg += freeall(c);
1782 ENDR("mult_apply_matrix_matrix");
1783 }
1784
1785
mult_matrix_matrix(a,b,c)1786 INT mult_matrix_matrix(a,b,c) OP a,b,c;
1787 /* AK 280388 */ /* AK 060789 V1.0 */ /* AK 111289 V1.1 */
1788 /* c = a * b */
1789 /* AK 210891 V1.3 */
1790 /* AK 221104 V3.0 c may be non empty */
1791 {
1792 INT erg = OK;
1793 CTO(MATRIX,"mult_matrix_matrix(1)",a);
1794 CTO(MATRIX,"mult_matrix_matrix(2)",b);
1795 {
1796 INT i,j,k,al,bl;
1797 OP z; /* zwischen ergebnis bei matrix-multiplikation */
1798 OP cp,ap,ap2,bp,bp2;
1799 al = S_M_LI(a);
1800 bl = S_M_HI(b);
1801 SYMCHECK( bl != al, "mult_matrix_matrix:different sizes");
1802 bl = S_M_LI(b);
1803
1804
1805
1806 erg += m_ilih_m(S_M_LI(b),S_M_HI(a),c);
1807 z=CALLOCOBJECT(); /* zwischensumme*/
1808 for (i=0,cp=S_M_S(c),ap=S_M_S(a);
1809 i<S_M_HI(a);i++,ap+=al) /* ueber zeilen der linken Matrix */
1810 for (j=0,bp=S_M_S(b);j<S_M_LI(b);
1811 j++,cp++,bp++) /* ueber spalten der rechten Matrix */
1812 {
1813 ap2 = ap;
1814 bp2 = bp;
1815 MULT(ap2,bp2,cp);
1816 ap2++;
1817 bp2 += bl;
1818 for (k=1;k<al;k++,ap2++,bp2+=bl)
1819 {
1820 FREESELF(z);
1821 MULT(ap2,bp2,z);
1822 ADD_APPLY(z,cp);
1823 };
1824 }
1825 FREEALL(z);
1826 }
1827 ENDR("mult_matrix_matrix");
1828 }
1829
1830
mult_matrix(a,b,d)1831 INT mult_matrix(a,b,d) OP a,b,d;
1832 /* AK 070789 V1.0 */ /* AK 310790 V1.1 */ /* AK 210891 V1.3 */
1833 {
1834 INT erg = OK;
1835 CTO(MATRIX,"mult_matrix(1)",a);
1836 EOP("mult_matrix(2)",b);
1837 CTO(EMPTY,"mult_matrix(3)",d);
1838
1839 switch(S_O_K(b))
1840 {
1841 case BRUCH:
1842 case INTEGER:
1843 case LONGINT:
1844 erg += mult_scalar_matrix(b,a,d);
1845 break;
1846
1847 case MATRIX:
1848 erg += mult_matrix_matrix(a,b,d);
1849 break;
1850 case VECTOR:
1851 erg += mult_matrix_vector(a,b,d); /* AK 200192 */
1852 break;
1853 default:
1854 WTO("mult_matrix(2)",b);
1855 break;
1856 }
1857 ENDR("mult_matrix");
1858 }
1859
1860
1861 #ifdef VECTORTRUE
mult_matrix_vector(b,a,c)1862 INT mult_matrix_vector(b,a,c) OP a, b, c;
1863 /* AK 200192 */
1864 /* AK 221104 V3.0 */
1865 {
1866 INT erg = OK;
1867 CTTO(INTEGERVECTOR,VECTOR,"mult_matrix_vector(2)",a);
1868 CTO(MATRIX,"mult_matrix_vector(1)",b);
1869 {
1870 INT i,j;
1871 OP d,h;
1872
1873 SYMCHECK (S_V_LI(a) != S_M_LI(b),"mult_matrix_vector:wrong size");
1874
1875 erg += m_il_v(S_M_HI(b),c);
1876 CALLOCOBJECT2(d,h);
1877 // erg += null(S_V_I(a,0),h);
1878 for (i=0;i<S_V_LI(c);i++)
1879 {
1880 // COPY(h,S_V_I(c,i));
1881 MULT(S_M_IJ(b,i,0),S_V_I(a,0),S_V_I(c,i));
1882 for (j=1;j<S_V_LI(a);j++)
1883 {
1884 erg += multadd_apply(S_M_IJ(b,i,j),S_V_I(a,j),S_V_I(c,i));
1885 /*
1886 FREESELF(d);
1887 MULT(S_M_IJ(b,i,j),S_V_I(a,j),d);
1888 ADD_APPLY(d,S_V_I(c,i));
1889 */
1890 }
1891 }
1892 FREEALL2(d,h);
1893 }
1894 ENDR("mult_matrix_vector");
1895 }
1896 #endif /* VECTORTRUE */
1897
1898
1899
mult_apply_matrix(a,b)1900 INT mult_apply_matrix(a,b) OP a,b;
1901 /* AK 131190 V1.1 */ /* AK 210891 V1.3 */
1902 {
1903 switch(S_O_K(b))
1904 {
1905 case MATRIX: return(mult_apply_matrix_matrix(a,b));
1906 default:
1907 {
1908 printobjectkind(b);
1909 error("mult_apply_matrix:wrong second type");
1910 return(ERROR);
1911 }
1912 }
1913 }
1914
1915
1916
1917
objectread_matrix(fp,matrix)1918 INT objectread_matrix(fp,matrix) FILE *fp; OP matrix;
1919 /* AK 300888 */ /* AK 070789 V1.0 */ /* AK 310790 V1.1 */
1920 /* AK 210891 V1.3 */
1921 {
1922 INT i,j;
1923 OP l= callocobject();
1924 OP h = callocobject();
1925 objectread(fp,h);
1926 objectread(fp,l);
1927 b_lh_m(l,h,matrix);
1928 for (i=0;i<S_M_HI(matrix); i++)
1929 for (j=0;j<S_M_LI(matrix); j++)
1930 objectread(fp,S_M_IJ(matrix,i,j));
1931 return(OK);
1932 }
1933
1934
1935
objectwrite_matrix(fp,matrix)1936 INT objectwrite_matrix(fp,matrix) FILE *fp; OP matrix;
1937 /* AK 300888 */ /* AK 070789 V1.0 */ /* AK 310790 V1.1 */
1938 /* AK 210891 V1.3 */
1939 {
1940 INT i,j;
1941
1942 fprintf(fp, " %ld ",MATRIX);
1943 objectwrite(fp,S_M_H(matrix));
1944 /* zuerst die hoehe */
1945 objectwrite(fp,S_M_L(matrix));
1946 /* dann die laenge */
1947
1948 for (i=0;i<S_M_HI(matrix); i++)
1949 for (j=0;j<S_M_LI(matrix); j++)
1950 objectwrite(fp,S_M_IJ(matrix,i,j));
1951 return(OK);
1952 }
1953
1954
1955
test_matrix()1956 INT test_matrix()
1957 /* AK 181289 V1.1 */
1958 /* AK 120891 V1.3 */
1959 {
1960 OP a = callocobject();
1961 OP b = callocobject();
1962
1963 printf("test_matrix:scan(a)");
1964 scan(MATRIX,a);
1965 println(a);
1966 printf("test_matrix:add(a,a,b)");
1967 add(a,a,b);
1968 println(b);
1969 printf("test_matrix:mult(a,b,b)");
1970 mult(a,b,b);
1971 println(b);
1972 printf("test_matrix:kronecker_product(a,b,b)");
1973 kronecker_product(a,b,b);
1974 println(b);
1975 #ifdef BRUCHTRUE
1976 printf("test_matrix:invers(b,a)");
1977 invers(b,a);
1978 println(a);
1979 #endif /* BRUCHTRUE */
1980 printf("test_matrix:delete_row_matrix(a,1L,b)");
1981 delete_row_matrix(a,1L,b);
1982 println(b);
1983 printf("test_matrix:delete_column_matrix(b,1L,b)");
1984 delete_column_matrix(b,1L,b);
1985 println(b);
1986
1987
1988 freeall(a);
1989 freeall(b);
1990 return(OK);
1991 }
1992
1993
1994
trace_matrix(a,b)1995 INT trace_matrix(a,b) OP a,b;
1996 /* AK 131289 spur einer matrix */
1997 /* AK 131289 V1.1 */ /* AK 270291 V1.2 */ /* AK 210891 V1.3 */
1998 {
1999 INT i,erg=OK;
2000 CTO(MATRIX,"trace_matrix(1)",a);
2001 CTO(EMPTY,"trace_matrix(2)",b);
2002 SYMCHECK (not quadraticp(a) ,"trace_matrix: matrix not quadratic");
2003
2004 null(S_M_S(a),b);
2005 for (i=S_M_HI(a)-1L; i>=0L; i--)
2006 ADD_APPLY(S_M_IJ(a,i,i),b);
2007 ENDR("trace_matrix");
2008 }
2009
2010
2011
symmetricp_matrix(a)2012 INT symmetricp_matrix(a) OP a;
2013 /* AK 150296 */
2014 {
2015 INT i,j;
2016 if (S_M_HI(a) != S_M_LI(a)) return FALSE;
2017 for (i=0L;i<S_M_HI(a) ; i++)
2018 for (j=0L;j<i;j++)
2019 if (neq(S_M_IJ(a,i,j),S_M_IJ(a,j,i))) return FALSE;
2020 return TRUE;
2021 }
2022
2023 #ifdef VECTORTRUE
spalten_summe(a,b)2024 INT spalten_summe(a,b) OP a,b;
2025 /* AK 281289 summe ueber die Spalten, ergebnis ist vector */
2026 /* AK 281289 V1.1 */ /* AK 240791 V1.3 */
2027 {
2028 INT i,j;
2029 INT erg = OK;
2030 CTTO(INTEGERMATRIX,MATRIX,"spalten_summe(1)",a);
2031
2032 erg += m_il_v(S_M_LI(a),b);
2033 for (j=0;j<S_M_LI(a);j++)
2034 {
2035 COPY(S_M_IJ(a,0,j),S_V_I(b,j));
2036 for (i=1;i<S_M_HI(a);i++)
2037 ADD_APPLY(S_M_IJ(a,i,j),S_V_I(b,j));
2038 }
2039 ENDR("spalten_summe");
2040 }
2041 #endif /* VECTORTRUE */
2042
2043
t_INTMATRIX_charvektor(a)2044 char * t_INTMATRIX_charvektor(a) OP a;
2045 /* AK 210891 V1.3 */
2046 {
2047 INT i,j,k=0;
2048 char *erg = (char *)
2049 SYM_malloc(S_M_HI(a)*S_M_LI(a)*sizeof(char));
2050 if(erg == NULL) error("t_INTMATRIX_charvektor:no memory");
2051 else {
2052 for (i=0;i<S_M_HI(a);i++)
2053 for (j=0;j<S_M_LI(a);j++,k++)
2054 erg[k] =(char)S_M_IJI(a,i,j);
2055 }
2056 return(erg);
2057 }
2058
2059
2060 #ifdef VECTORTRUE
m_vector_diagonalmatrix(a,b)2061 INT m_vector_diagonalmatrix(a,b) OP a,b;
2062 /* AK 171290 V1.1 */ /* AK 110691 V1.2 */ /* AK 210891 V1.3 */
2063 {
2064 INT i;
2065 m_lh_nm(S_V_L(a),S_V_L(a),b);
2066 for (i=0L; i<S_M_HI(b);i++)
2067 copy(S_V_I(a,i),S_M_IJ(b,i,i));
2068 return OK;
2069 }
2070 #endif /* VECTORTRUE */
2071
2072
max_matrix(a,b)2073 INT max_matrix(a,b) OP a,b;
2074 /* b becomes copy of the maximum entry */
2075 /* a and b may be equal */
2076 /* empty objects are not used for comparision */
2077 /* AK 110691 V1.2 */ /* AK 210891 V1.3 */
2078 /* AK 090804 V3.0 */
2079 {
2080 OP z = S_M_S(a),zb = S_M_S(a);
2081 INT i;
2082
2083 for (i=S_M_HI(a)*S_M_LI(a); i>0L; i--, z++)
2084 if (not EMPTYP(z))
2085 {
2086 if (EMPTYP(zb)) zb = z;
2087 else if (GR(z,zb)) zb = z;
2088 }
2089 return copy(zb,b);
2090 }
2091
min_matrix(a,b)2092 INT min_matrix(a,b) OP a,b;
2093 /* b becomes copy of the minimum entry */
2094 /* a and b may be equal */
2095 /* AK 140703 */
2096 /* AK 061207 */
2097 {
2098 OP z = S_M_S(a),zb = NULL;
2099 INT i;
2100
2101 for (i=S_M_HI(a)*S_M_LI(a); i>0L; i--, z++)
2102 {
2103 if (not EMPTYP(z))
2104 {
2105 if (zb==NULL) zb = z;
2106 else if (LT(z,zb)) zb = z;
2107 }
2108 }
2109 if (zb==NULL) return error("min_matrix: no entries");
2110 else return copy(zb,b);
2111 }
2112
zeilen_summe(a,b)2113 INT zeilen_summe(a,b) OP a,b;
2114 /* AK 281289 summe ueber die Zeilen, ergebnis ist vector */
2115 /* AK 281289 V1.1 */ /* AK 210891 V1.3 */
2116 {
2117 INT i,j;
2118 INT erg = OK;
2119 CTTO(INTEGERMATRIX,MATRIX,"zeilen_summe(1)",a);
2120
2121 erg += m_il_nv(S_M_HI(a),b);
2122 for (j=0;j<S_M_HI(a);j++)
2123 for (i=0;i<S_M_LI(a);i++)
2124 ADD_APPLY(S_M_IJ(a,j,i),S_V_I(b,j));
2125 ENDR("zeilen_summe");
2126 }
2127
2128
n_fold_kronecker_product(n,a,b)2129 INT n_fold_kronecker_product(n,a,b) OP n; OP a; OP b;
2130 /* RH 080891 V1.3 */
2131 {
2132 INT i;
2133
2134 if(S_I_I(n)>= 2L) kronecker_product(a,a,b);
2135 for(i=2;i<S_I_I(n);++i) kronecker_product(a,b,b);
2136 return OK;
2137 }
2138
2139
kronecker_product(a,b,c)2140 INT kronecker_product(a,b,c) OP a,b,c;
2141 /* RH 080891 V1.3 */
2142 {
2143 INT i;
2144 INT j;
2145 INT k;
2146 INT l,erg=OK;
2147
2148 OP H,dim;
2149
2150 CTTO(INTEGERMATRIX,MATRIX,"kronecker_product(1)",a);
2151 CTTO(INTEGERMATRIX,MATRIX,"kronecker_product(2)",b);
2152 CE3(a,b,c,kronecker_product);
2153
2154 H = callocobject();
2155 dim = callocobject();
2156
2157 MULT(S_M_L(a),S_M_L(b),dim);
2158
2159 erg += m_lh_m(dim,dim,c);
2160 erg += m_lh_m(S_M_L(a),S_M_L(b),H);
2161
2162 for(i=0L;i<S_M_LI(a);++i)
2163 {
2164 for(j=0L;j<S_M_HI(a);++j)
2165 {
2166 erg += mult(S_M_IJ(a,i,j),b,H);
2167 for(k=0L;k<S_M_LI(H);++k)
2168 for(l=0L;l<S_M_HI(H);++l)
2169 COPY(S_M_IJ(H,k,l), S_M_IJ(c,i*S_M_HI(b)+k,j*S_M_LI(H)+l));
2170 }
2171 }
2172
2173 FREEALL2(H,dim);
2174 ENDR("kronecker_product");
2175 }
2176
select_row(a,index,b)2177 INT select_row(a,index,b) OP a,b; INT index;
2178 /* transfer one row to vector object */
2179 /* AK 110592 */
2180 {
2181 INT erg = OK;
2182 CTTTO(INTEGERMATRIX,MATRIX,TABLEAUX,"select_row(1)",a);
2183 SYMCHECK(index < 0,"select_row: negative index");
2184 if (a == b) /* AK 210802 */
2185 {
2186 OP c;
2187 c = CALLOCOBJECT();
2188 SWAP(a,c);
2189 erg += select_row(c,index,b);
2190 FREEALL(c);
2191 }
2192 else if (TABLEAUXP(a)) /* AK 280193 */
2193 {
2194 erg += select_row_tableaux(a,index,b);
2195 }
2196 else
2197 {
2198 INT j;
2199 SYMCHECK(index >= S_M_HI(a), "select_row: index >= height");
2200 erg += m_il_v(S_M_LI(a),b);
2201 for (j=0L;j<S_M_LI(a);j++)
2202 COPY(S_M_IJ(a,index,j),S_V_I(b,j));
2203 }
2204 ENDR("select_row");
2205 }
2206
2207
select_column(a,index,b)2208 INT select_column(a,index,b) OP a,b; INT index;
2209 /* transfer one column to vector object */
2210 /* AK 110592 */
2211 /* AK 090107 V3.1 */
2212 {
2213 INT erg = OK;
2214 CTTTO(INTEGERMATRIX,MATRIX,TABLEAUX,"select_column(1)",a);
2215 SYMCHECK(index < 0,"select_column: negative index");
2216
2217 if (a == b) /* AK 210802 */
2218 {
2219 OP c;
2220 c = CALLOCOBJECT();
2221 SWAP(a,c);
2222 erg += select_column(c,index,b);
2223 FREEALL(c);
2224 }
2225 else if (TABLEAUXP(a)) /* AK 280193 */
2226 {
2227 erg += select_column_tableaux(a,index,b);
2228 }
2229 else {
2230 INT j;
2231 SYMCHECK(index >= S_M_LI(a), "select_column: index >= length");
2232 erg += m_il_v(S_M_HI(a),b);
2233 for (j=0;j<S_M_HI(a);j++)
2234 COPY(S_M_IJ(a,j,index),S_V_I(b,j));
2235 }
2236 ENDR("select_column");
2237 }
2238
2239
rank(a,b)2240 INT rank(a,b) OP a,b;
2241 /* AK 220294 */
2242 /* rank of a matrix */
2243 /* a and b may be equal */
2244 /* AK 141204 V3.0 */
2245 /* AK 2221106 V3.1 */
2246 {
2247 INT erg = OK;
2248 CTTO(MATRIX,INTEGERMATRIX,"rank(1)",a);
2249 {
2250 INT i,j,k,l,z=0;
2251 INT res = 0;
2252 OP e,c,d;
2253 CALLOCOBJECT3(c,d,e);
2254 COPY(a,e);
2255 a = e;
2256 j=0;
2257 a2:
2258 for (i=0;i<S_M_HI(a);i++)
2259 {
2260 if (not NULLP(S_M_IJ(a,i,j)))
2261 goto a1;
2262 }
2263 j++;
2264 if (j == S_M_LI(a))
2265 goto a3;
2266 goto a2;
2267 a1:
2268 /* i , j ist nicht null */
2269 if (i != z)
2270 erg += change_row_ij(a,i,z); /* AK 080394 */
2271 i = z++;
2272
2273 res ++;
2274 for (k = 0;k<S_M_HI(a);k++)
2275 if (k != i)
2276 {
2277 if (S_O_K(S_M_IJ(a,k,j))==GALOISRING) /* kein invers */
2278 {
2279 OP f1=callocobject();
2280 OP f2=callocobject();
2281 copy(S_M_IJ(a,k,j),f1); copy(S_M_IJ(a,i,j),f2);
2282 /* zeile k soll null werden in der spalte i */
2283 /* dazu f1*zeile_i - f2*zeile_k in spalte k schreiben */
2284 for (l=0;l<S_M_LI(a);l++)
2285 {
2286 CLEVER_MULT(S_M_IJ(a,i,l),f1,c);
2287 CLEVER_MULT(S_M_IJ(a,k,l),f2,d); ADDINVERS_APPLY(d);
2288 ADD(c,d,S_M_IJ(a,k,l));
2289 }
2290 FREEALL2(f1,f2);
2291 }
2292 else {
2293 erg += div(S_M_IJ(a,k,j),S_M_IJ(a,i,j),c);
2294 for (l=0;l<S_M_LI(a);l++)
2295 {
2296 CLEVER_MULT(c,S_M_IJ(a,i,l),d);
2297 ADDINVERS_APPLY(d);
2298 ADD_APPLY(d,S_M_IJ(a,k,l));
2299 }
2300 }
2301 }
2302 for (l=j+1;l<S_M_LI(a);l++)
2303 {
2304 if (S_O_K(S_M_IJ(a,i,0))==INTEGER) M_I_I(0,S_M_IJ(a,i,l));
2305 else null(S_M_IJ(a,i,0), S_M_IJ(a,i,l));
2306 }
2307 j++;
2308 if (j != S_M_LI(a))
2309 goto a2;
2310 a3:
2311 erg += m_i_i(res,b);
2312 FREEALL3(c,d,a);
2313 }
2314 ENDR("rank");
2315 }
2316
bideterminant_vector(a,b,c)2317 INT bideterminant_vector(a,b,c) OP a,b,c;
2318 /* AK 010802 */
2319 {
2320 INT i,erg=OK;
2321 CTO(VECTOR,"bideterminant_vector(1)",a);
2322 CTO(VECTOR,"bideterminant_vector(2)",b);
2323 SYMCHECK( NEQ(S_V_L(a),S_V_L(b)), "bideterminant_vector:vector of different length");
2324
2325 {
2326 OP res;
2327 INT i,j;
2328 res = CALLOCOBJECT();
2329 erg += m_ilih_m(S_V_LI(a),S_V_LI(b),res);
2330 for (i=0;i<S_M_HI(res);i++)
2331 for (j=0;j<S_M_HI(res);j++)
2332 {
2333 OP poly;
2334 poly = S_M_IJ(res,i,j);
2335 erg += b_skn_po(CALLOCOBJECT(),CALLOCOBJECT(),NULL,poly);
2336 M_I_I(1,S_PO_K(poly));
2337 erg += m_ilih_nm( S_V_II(b,j)+1, S_V_II(a,i)+1,S_PO_S(poly));
2338 C_O_K(S_PO_S(poly),INTEGERMATRIX);
2339 M_I_I(1,S_M_IJ(S_PO_S(poly),S_V_II(a,i), S_V_II(b,j)));
2340 }
2341 erg += det_mat_imm(res,c);
2342 FREEALL(res);
2343 }
2344
2345 ENDR("bideterminant_vector");
2346 }
2347
bideterminant_tableaux(a,b,c)2348 INT bideterminant_tableaux(a,b,c) OP a,b,c;
2349 /* AK 010802 */
2350 {
2351 INT i,erg=OK;
2352 CTO(TABLEAUX,"bideterminant_tableaux(1)",a);
2353 CTO(TABLEAUX,"bideterminant_tableaux(2)",b);
2354 SYMCHECK( NEQ(S_T_U(a),S_T_U(b)), "bideterminant_tableaux:tableaux of different shape");
2355
2356 erg += m_i_i(1,c);
2357 for (i=0;i<S_T_LI(a);i++)
2358 {
2359 OP row1,row2;
2360 OP res;
2361 row1 = CALLOCOBJECT();
2362 row2 = CALLOCOBJECT();
2363 erg += select_column(a,i,row1);
2364 erg += select_column(b,i,row2);
2365 res = CALLOCOBJECT();
2366 erg += bideterminant_vector(row1,row2,res);
2367 FREEALL(row1);
2368 FREEALL(row2);
2369 MULT_APPLY(res,c);
2370 FREEALL(res);
2371 }
2372 ENDR("bideterminant_tableaux");
2373 }
2374
bideterminant(a,b,c)2375 INT bideterminant(a,b,c) OP a,b,c;
2376 /* AK 010802 */
2377 {
2378 INT erg=OK;
2379 switch(S_O_K(a)){
2380 case TABLEAUX:erg += bideterminant_tableaux(a,b,c);break;
2381 case VECTOR:erg += bideterminant_vector(a,b,c);break;
2382 default:WTO("bideterminant(1)",a);break;
2383 }
2384 ENDR("bideterminant");
2385 }
2386
P_symmetrized_bideterminant(a,b,c)2387 INT P_symmetrized_bideterminant(a,b,c) OP a,b,c;
2388 /* a and b tableaux of the same shape */
2389 /* c becomes P-symmetrized bideterminant */
2390 /* AK 060802 */
2391 /* cf. Golembiowski Bayreuther Mathematische Schriften */
2392 {
2393 INT erg = OK;
2394 CTO(TABLEAUX,"P_symmetrized_bideterminant(1)",a);
2395 CTO(TABLEAUX,"P_symmetrized_bideterminant(2)",b);
2396 SYMCHECK(NEQ(S_T_U(a),S_T_U(b)),"P_symmetrized_bideterminant:different shapes of tableaux");
2397 {
2398 OP d,e;
2399 d = CALLOCOBJECT();
2400 e = CALLOCOBJECT();
2401 first_lex_tableaux(b,e);
2402 init(POLYNOM,c);
2403 do {
2404 OP f;
2405 f = CALLOCOBJECT();
2406 erg += bideterminant(a,e,f);
2407 insert(f,c,NULL,NULL);
2408 erg += copy(e,d);
2409 } while(next_lex_tableaux(d,e) == TRUE);
2410 FREEALL(d);
2411 FREEALL(e);
2412 }
2413 CTO(POLYNOM,"P_symmetrized_bideterminant(3-end)",c);
2414 ENDR("P_symmetrized_bideterminant");
2415 }
2416
operate_perm_spaltenmatrix(a,b,c)2417 INT operate_perm_spaltenmatrix(a,b,c) OP a,b,c;
2418 /* vertausch spalten gemaess der permutation */
2419 /* spalte 0 nur wenn permutation auch 0 enthaelt */
2420 /* AK 080802 */
2421 {
2422 INT erg = OK;
2423 INT i,j;
2424 CTO(PERMUTATION,"operate_perm_spaltenmatrix(1)",a);
2425 CTTO(INTEGERMATRIX,MATRIX,"operate_perm_spaltenmatrix(2)",b);
2426 CE3(a,b,c,operate_perm_spaltenmatrix);
2427 SYMCHECK(S_P_LI(a) >= S_M_LI(b),
2428 "operate_perm_spaltenmatrix: permutation degree too big");
2429 COPY(b,c);
2430 for (j=0;j<S_P_LI(a);j++)
2431 for (i=0;i<S_M_HI(b);i++)
2432 CLEVER_COPY(S_M_IJ(b,i,j+1), S_M_IJ(c,i,S_P_II(a,j)));
2433 ENDR("operate_perm_spaltenmatrix");
2434 }
2435
operate_perm_bideterminant(a,b,c)2436 INT operate_perm_bideterminant(a,b,c) OP a,b,c;
2437 /* AK 080802 */
2438 /* permutation operates on the bideterminant by exchanging columns */
2439 {
2440 INT erg = OK;
2441 OP z,m;
2442 CTO(PERMUTATION,"operate_perm_bideterminant(1)",a);
2443 CTO(POLYNOM,"operate_perm_bideterminant(2)",b);
2444 CE3(a,b,c,operate_perm_bideterminant);
2445 if (NULLP(b)) { erg += copy(b,c); goto endr_ende; }
2446 CTTO(INTEGERMATRIX,MATRIX,"operate_perm_bideterminant(2-self)",S_PO_S(b));
2447 SYMCHECK(S_P_LI(a) >= S_M_LI(S_PO_S(b)),
2448 "operate_perm_bideterminant: permutation degree too big");
2449
2450 {
2451 OP d;
2452 INT i;
2453 d = CALLOCOBJECT();
2454 t_LIST_VECTOR(b,d);
2455 for(i=0;i<S_V_LI(d);i++)
2456 erg += operate_perm_spaltenmatrix(
2457 a,S_MO_S(S_V_I(d,i)),S_MO_S(S_V_I(d,i)));
2458 qsort_vector(d);
2459 erg += t_VECTOR_POLYNOM(d,c);
2460 FREEALL(d);
2461 }
2462 ENDR("operate_perm_bideterminant");
2463 }
2464
2465
append_behind_matrix_matrix(OP a,OP b,OP c)2466 INT append_behind_matrix_matrix(OP a, OP b, OP c)
2467 /* AK 030406 V3.1 */
2468 /* a and b have same height */
2469 {
2470 INT erg = OK;
2471 CTO(MATRIX,"append_behind_matrix_matrix(1)",a);
2472 CTO(MATRIX,"append_behind_matrix_matrix(2)",b);
2473 SYMCHECK(S_M_HI(a)!=S_M_HI(b),"append_behind_matrix_matrix: different sizes");
2474 CE3(a,b,c,append_behind_matrix_matrix);
2475 {
2476 INT i,j;
2477 erg += m_ilih_m(S_M_LI(a)+S_M_LI(b),S_M_HI(a),c);
2478 for (i=0;i<S_M_HI(c);i++)
2479 for (j=0;j<S_M_LI(a);j++)
2480 COPY(S_M_IJ(a,i,j),S_M_IJ(c,i,j));
2481 for (i=0;i<S_M_HI(c);i++)
2482 for (j=0;j<S_M_LI(b);j++)
2483 COPY(S_M_IJ(b,i,j),S_M_IJ(c,i,j+S_M_LI(a)));
2484 }
2485 ENDR("append_behind_matrix_matrix");
2486 }
2487
append_below_matrix_matrix(OP a,OP b,OP c)2488 INT append_below_matrix_matrix(OP a, OP b, OP c)
2489 /* AK 120406 V3.1 */
2490 /* a and b have same length */
2491 {
2492 INT erg = OK;
2493 CTO(MATRIX,"append_below_matrix_matrix(1)",a);
2494 CTO(MATRIX,"append_below_matrix_matrix(2)",b);
2495 SYMCHECK(S_M_LI(a)!=S_M_LI(b),"append_below_matrix_matrix: different sizes");
2496 CE3(a,b,c,append_below_matrix_matrix);
2497 {
2498 INT i,j;
2499 erg += m_ilih_m(S_M_LI(a) , S_M_HI(a)+S_M_HI(b),c);
2500 for (i=0;i<S_M_HI(a);i++)
2501 for (j=0;j<S_M_LI(a);j++)
2502 COPY(S_M_IJ(a,i,j),S_M_IJ(c,i,j));
2503 for (i=0;i<S_M_HI(b);i++)
2504 for (j=0;j<S_M_LI(b);j++)
2505 COPY(S_M_IJ(b,i,j),S_M_IJ(c,i+S_M_HI(a),j));
2506 }
2507 ENDR("append_below_matrix_matrix");
2508 }
2509
m_matrix_column_vector(OP a,OP b)2510 INT m_matrix_column_vector(OP a, OP b)
2511 /* AK 231107 */
2512 /* becomes vector of columns = vector of vector */
2513 {
2514 INT erg = OK;
2515 CTO(MATRIX,"m_matrix_column_vector(1)",a)
2516 CE2(a,b,m_matrix_column_vector);
2517 {
2518 INT i;
2519 erg += m_il_v(S_M_LI(a),b);
2520 for (i=0;i<S_V_LI(b);i++) erg += select_column(a,i,S_V_I(b,i));
2521 }
2522 ENDR("m_matrix_column_vector");
2523 }
2524
2525 #endif /* MATRIXTRUE */
2526