?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
add $(EXEEXT) to executable targets during installation for MinGW
[suif.git] / src / baseparsuif / suifmath / int_matrix.cc
blob2fd2caeac171a31eebc04b0b16fcc3b750d28953
1 /* file "int_matrix.cc" */
2
3 /* Copyright (c) 1994 Stanford University
4
5 All rights reserved.
6
7 This software is provided under the terms described in
8 the "suif_copyright.h" include file. */
9
10 #include <suif_copyright.h>
11
12 #define _MODULE_ "libsuifmath.a"
13 #pragma implementation "int_matrix.h"
14
15 #include <suif1.h>
16
17 #include "suifmath.h"
18 #include <string.h>
19
20
21 class coeff;
22
23 int compose_print_flag = FALSE;
24
25 void init_suifmath(int & /* argc */, char * /* argv */ [])
26 {
27 // k_named_sc_merge_mark_annote = lexicon->enter("named_sc_merge_mark")->sp;
28 ANNOTE(k_named_sc_merge_mark_annote, "named_sc_merge_mark", FALSE);
29
30 STRUCT_ANNOTE(code_context_annote::k_annote,
31 "code_context",
32 TRUE,
33 code_context_annote_from,
34 code_context_annote_to,
35 code_context_annote_free,
36 code_context_annote_print);
37 }
38
39 void exit_suifmath(void)
40 {
41
42 }
43
44 /********************************************************************************
45 * gcd (Gratest Common Dinominator) *
46 ********************************************************************************/
47 int gcd(int a, int b)
48 {
49 assert((a>=0)&&(b>=0));
50
51 if(a==0)
52 return b;
53 else if (b==0)
54 return a;
55 else if(b>a)
56 return gcd(b%a, a);
57 else
58 return gcd(a%b, b);
59 }
60
61 /********************************************************************************
62 * lcm (Least Common Multiple) *
63 ********************************************************************************/
64 int lcm(int a, int b)
65 {
66 assert((a>=0)&&(b>=0));
67 if(a==b) return a;
68 long gcdab = gcd(a,b);
69 long aL = a/gcdab;
70 long bL = b;
71 long out = aL*bL;
72 assert(out/bL == aL); // check for overflow
73 return out;
74 }
75
76
77
78 /* ##################################################
79 ##### integer_row #####
80 ################################################## */
81
82 /********************************************************************************
83 * initializers *
84 ********************************************************************************/
85 void integer_row::mknew(int s)
86 {
87 if(rsz < s) {
88 clear();
89 rsz = s;
90 R = new int[rsz];
91 };
92 sz = s;
93 }
94
95 integer_row::integer_row()
96 {
97 R=0;
98 sz=0;
99 rsz=0;
100 }
101
102 integer_row::integer_row(const integer_row & rw)
103 {
104 R=0;
105 sz=0;
106 rsz=0;
107 init(rw);
108 }
109
110 integer_row::integer_row(const integer_row * rw)
111 {
112 R=0;
113 sz=0;
114 rsz=0;
115 init(rw);
116 }
117
118 integer_row::integer_row(int s)
119 {
120 R=0;
121 sz=0;
122 rsz=0;
123 init(s);
124 }
125
126 void integer_row::clear()
127 {
128 if(R) {
129 delete[] R;
130 R=NULL;
131 }
132 sz=0;
133 rsz=0;
134 }
135
136 void integer_row::init(const integer_row & rw)
137 {
138 if (R) {
139 if(rw.n() <= rsz) {
140 sz = rw.n();
141 memcpy(R, rw.R, n()*sizeof(int));
142 return;
143 } else
144 clear();
145 }
146
147 rsz = sz = rw.n();
148 if(n() > 0) {
149 R = new int[n()];
150 memcpy(R, rw.R, n()*sizeof(int));
151 }
152 }
153
154
155 void integer_row::init(int s)
156 {
157 assert(s > 0);
158 if (R) {
159 if(s <= rsz) {
160 sz = s;
161 memset(R, 0, n()*sizeof(int));
162 return;
163 } else
164 clear();
165 }
166
167 rsz = sz = s;
168 R = new int[n()];
169 memset(R, 0, n()*sizeof(int));
170 }
171
172 void integer_row::init(const int * ilist, int s)
173 {
174 assert(s > 0);
175 if (R) {
176 if(s <= rsz) {
177 sz = s;
178 memcpy(R, ilist, n()*sizeof(int));
179 return;
180 } else
181 clear();
182 }
183
184 rsz = sz = s;
185 R = new int[n()];
186 memcpy(R, ilist, n()*sizeof(int));
187 }
188
189
190
191
192 /********************************************************************************
193 * row operations *
194 ********************************************************************************/
195 boolean integer_row::operator==(const integer_row & a) const
196 {
197 return ((a.n() == n()) &&
198 (memcmp(&a.R[0], &R[0], sizeof(int)*n()) == 0));
199
200
201 /* Speed hack above
202 for(int i=0; i<n(); i++)
203 if(a.R[i] != R[i]) return 0;
204 return 1;
205 */
206 }
207
208
209
210 integer_row & integer_row::operator=(const integer_row & row)
211 {
212 if(&row == this) return *this;
213
214 mknew(row.n());
215
216 memcpy(R, row.R, n()*sizeof(int));
217
218 return *this;
219 }
220
221
222 // Swap rows
223 void integer_row::operator^=(integer_row & row)
224 {
225 // assert(n() == row.n());
226 int *tmp = R;
227 R = row.R;
228 row.R = tmp;
229
230 int tmp2 = sz;
231 sz = row.n();
232 row.sz = tmp2;
233
234 tmp2 = rsz;
235 rsz = row.rsz;
236 row.rsz = tmp2;
237 }
238
239 integer_row integer_row::swap_col(int i, int j) const
240 {
241 integer_row R(this);
242 R.do_swap_col(i, j);
243 return R;
244 }
245
246
247 void integer_row::do_swap_col(int i, int j)
248 {
249 if(i == j) return;
250
251 assert((i>=0)&&(i<n()));
252 assert((j>=0)&&(j<n()));
253
254
255 int tmp = (*this)[i];
256 (*this)[i] = (*this)[j];
257 (*this)[j] = tmp;
258 }
259
260
261
262 integer_row integer_row::operator+(const integer_row & row) const
263 {
264 assert(n() == row.n());
265 integer_row result(n());
266
267 for(int i=0; i<n(); i++)
268 result[i] = R[i] + row.c(i);
269
270 return result;
271 }
272
273 integer_row integer_row::operator-(const integer_row & row) const
274 {
275 assert(n() == row.n());
276
277 integer_row result(n());
278
279 for(int i=0; i<n(); i++)
280 result[i] = R[i] - row.c(i);
281
282 return result;
283 }
284
285
286 integer_row integer_row::operator*(const integer_row & row) const
287 {
288 assert(n() == row.n());
289
290 integer_row result(n());
291
292 for(int i=0; i<n(); i++)
293 result[i] = R[i] * row.c(i);
294
295 return result;
296 }
297
298
299 integer_row integer_row::operator/(const integer_row & row) const
300 {
301 assert(n() == row.n());
302
303 integer_row result(n());
304 for(int i=0; i<n(); i++)
305 result[i] = R[i] / row.c(i);
306
307 return result;
308 }
309
310
311 // shuffle operator ret[x] = this[row[x]]
312 integer_row integer_row::operator%(const integer_row & row) const
313 {
314 assert(n() == row.n());
315
316 integer_row result(n());
317 result = 0;
318 for(int i=0; i<n(); i++) {
319 assert((row.R[i] >= 0)&&(row.R[i] < n()));
320 result[i] = R[row.R[i]];
321 }
322
323 return result;
324 }
325
326
327 integer_row & integer_row::operator+=(const integer_row & row)
328 {
329 assert(n() == row.n());
330
331 for(int i=0; i<n(); i++)
332 R[i] += row.R[i];
333
334 return *this;
335 }
336
337 integer_row & integer_row::operator-=(const integer_row & row)
338 {
339 assert(n() == row.n());
340
341 for(int i=0; i<n(); i++)
342 R[i] -= row.R[i];
343
344 return *this;
345 }
346
347
348 integer_row & integer_row::operator*=(const integer_row & row)
349 {
350 assert(n() == row.n());
351
352 for(int i=0; i<n(); i++)
353 R[i] *= row.R[i];
354
355 return *this;
356 }
357
358
359 integer_row & integer_row::operator/=(const integer_row & row)
360 {
361 assert(n() == row.n());
362
363 for(int i=0; i<n(); i++)
364 R[i] /= row.c(i);
365
366 return *this;
367 }
368
369 void integer_row::do_addmul(const integer_row & row, int mul)
370 {
371 assert(n() == row.n());
372
373 if (mul == 1)
374 *this += row;
375 else
376 for(int i=0; i<n(); i++)
377 R[i] += mul*row.c(i);
378 }
379
380 integer_row & integer_row::operator=(int val)
381 {
382 for(int i=0; i<n(); i++)
383 R[i] = val;
384 return *this;
385 }
386
387 integer_row & integer_row::operator+=(int val)
388 {
389 for(int i=0; i<n(); i++)
390 R[i] += val;
391 return *this;
392 }
393
394 integer_row & integer_row::operator-=(int val)
395 {
396 for(int i=0; i<n(); i++)
397 R[i] -= val;
398 return *this;
399 }
400
401 integer_row & integer_row::operator*=(int val)
402 {
403 for(int i=0; i<n(); i++)
404 R[i] *= val;
405 return *this;
406 }
407
408 integer_row & integer_row::operator/=(int val)
409 {
410 for(int i=0; i<n(); i++)
411 R[i] /= val;
412 return *this;
413 }
414
415
416
417 /********************************************************************************
418 * delete column *
419 ********************************************************************************/
420 integer_row integer_row::del_col(int i, int j) const
421 {
422 assert(i<= j);
423 assert((i>=0)&&(j<n()));
424 integer_row ret(n() - (j - i + 1));
425
426 int cnt = 0;
427 int a;
428 for(a = 0; a < i; a++)
429 ret[cnt++] = R[a];
430
431 for(a = j+1; a < n(); a++)
432 ret[cnt++] = R[a];
433
434 assert(cnt == ret.n());
435
436 return ret;
437 }
438
439 void integer_row::do_del_col(int i, int j)
440 {
441 assert(i<= j);
442 assert((i>=0)&&(j<n()));
443
444 int st, fn;
445 for (st = i, fn = j+1; (fn < sz); st++, fn++)
446 R[st] = R[fn];
447
448 sz = sz - (j - i + 1);
449 }
450
451 void integer_row::do_del_columns(const integer_row & mask)
452 {
453 int src, dst;
454 assert(mask.n() == n());
455
456 for (dst=0, src=0; (src < n()); src++) {
457 if (!mask.R[src]) {
458 R[dst] = R[src];
459 dst++;
460 };
461 }
462
463 sz = dst;
464 }
465
466 /********************************************************************************
467 * insert column *
468 ********************************************************************************/
469 integer_row integer_row::insert_col(int i) const
470 {
471 assert((i>=0)&&(i<=n()));
472 integer_row ret(n()+1);
473
474 int cnt = 0;
475 int a;
476 for(a = 0; a < i; a++)
477 ret[cnt++] = R[a];
478
479 ret[cnt++] = 0;
480
481 for(a = i; a < n(); a++)
482 ret[cnt++] = R[a];
483
484 assert(cnt == ret.n());
485
486 return ret;
487 }
488
489 void integer_row::do_insert_col(int i)
490 {
491 assert((i>=0)&&(i<=n()));
492 if (sz < rsz) {
493 int src, dst;
494 sz++;
495 for (src=sz-2, dst=sz-1; src >= i; --dst, --src)
496 R[dst] = R[src];
497 R[i] = 0;
498 } else {
499 assert((i>=0)&&(i<=n()));
500 integer_row ret(MAX(n()+1,n()*2));
501 ret.sz = n()+1;
502
503 int cnt = 0;
504 int a;
505 for(a = 0; a < i; a++)
506 ret[cnt++] = R[a];
507
508 ret[cnt++] = 0;
509
510 for(a = i; a < n(); a++)
511 ret[cnt++] = R[a];
512
513 assert(cnt == ret.n());
514
515 *this ^= ret;
516 }
517 return;
518 }
519
520 int integer_row::hashkey() const
521 {
522 int h = 0;
523 for(int i=0; i<n(); i++)
524 h += R[i]*i*i;
525 return h;
526 }
527
528 /* ##################################################
529 ##### integer_matrix #####
530 ################################################## */
531
532
533 integer_row *integer_matrix::alloc_rows(int num_rows)
534 {
535 realmm = num_rows;
536 return new integer_row[num_rows];
537 }
538
539
540 /********************************************************************************
541 * constructors *
542 ********************************************************************************/
543 integer_matrix::integer_matrix(int rows,int cols)
544 {
545 A=0;
546 nn = mm = realmm = 0;
547 init(rows, cols);
548 }
549
550
551 integer_matrix::integer_matrix(const integer_matrix *m)
552 {
553 A = 0;
554 nn = mm = realmm = 0;
555 init(m);
556 }
557
558
559 integer_matrix::integer_matrix(const integer_matrix & m)
560 {
561 A = 0;
562 nn = mm = realmm = 0;
563 init(m);
564 }
565
566 integer_matrix::integer_matrix(const integer_matrix *m, int rows)
567 {
568 A = 0;
569 nn = mm = realmm = 0;
570 init(m, rows);
571 }
572
573 integer_matrix::integer_matrix(const integer_matrix & m, int rows)
574 {
575 A = 0;
576 nn = mm = realmm = 0;
577 init(m, rows);
578 }
579
580 integer_matrix::integer_matrix(const coeff * c)
581 {
582 A = 0;
583 nn = mm = realmm = 0;
584 init(c);
585 }
586
587 integer_matrix::integer_matrix()
588 {
589 A = 0;
590 mm = nn = realmm = 0;
591 }
592
593 integer_matrix::~integer_matrix()
594 {
595 if(A) {
596 delete[] A;
597 A = 0;
598 mm = nn = realmm = 0;
599 }
600 }
601
602
603 void integer_matrix::clear()
604 {
605 if(A) {
606 delete[] A;
607 A = 0;
608 mm = nn = realmm = 0;
609 }
610 assert((mm == 0) && (realmm == 0));
611 }
612
613 /********************************************************************************
614 * init *
615 ********************************************************************************/
616
617 void integer_matrix::mknew(int rows, int cols)
618 {
619 assert((rows >= 0)&&(cols >= 0));
620
621 if (A) {
622 if((rows<=realmm)&&(cols<=n())) {
623 mm = rows; nn = cols;
624 for(int i=0; i<mm; i++) A[i].mknew(cols);
625 return;
626 } else
627 clear();
628 }
629
630 if(rows > 0) {
631 A = alloc_rows(rows);
632 for(int i=0; i<rows; i++) A[i].mknew(cols);
633 }
634 nn=cols;
635 mm=rows;
636 }
637
638 void integer_matrix::init(int rows,int cols)
639 {
640 assert((rows >= 0)&&(cols >= 0));
641
642 if (A) {
643 if((rows<=realmm)&&(cols<=n())) {
644 mm = rows; nn = cols;
645 for(int i=0; i<mm; i++) A[i].init(cols);
646 return;
647 } else
648 clear();
649 }
650
651 if(rows > 0) {
652 A = alloc_rows(rows);
653 for(int i=0; i<rows; i++) A[i].init(cols);
654 }
655 nn=cols;
656 mm=rows;
657 }
658
659
660 void integer_matrix::init(const integer_matrix & M)
661 {
662 assert(&M);
663
664 if (A) {
665 if((M.m()<=realmm)&&(M.n()<=n())) {
666 mm = M.m(); nn = M.n();
667 for(int i=0; i<m(); i++) A[i] = M.r(i);
668 return;
669 } else
670 clear();
671 }
672
673 if(M.m() > 0) {
674 assert(M.n()>0);
675 A = alloc_rows(M.m());
676 for(int i=0; i<M.m(); i++) A[i] = M.r(i);
677 }
678
679 nn=M.n();
680 mm=M.m();
681 }
682
683
684
685 void integer_matrix::init(const integer_matrix & M, int rows)
686 {
687 assert(&M);
688 assert(rows >= 0);
689
690 if (A) {
691 if((M.m()<=realmm)&&(M.n()<=n())) {
692 mm = M.m(); nn = M.n();
693 for(int i=0; i<m(); i++) A[i] = M.r(i);
694 return;
695 } else
696 clear();
697 }
698
699 nn=M.nn;
700 mm=rows;
701 A = alloc_rows(m());
702
703 if(m() <= M.m())
704 for(int i=0; i<m(); i++) A[i] = M.r(i);
705 else {
706 int i;
707 for(i=0; i< M.m(); i++) A[i] = M.r(i);
708 for(i=M.m(); i< m(); i++) A[i].init(n());
709 }
710 }
711
712
713 void integer_matrix::init(const coeff * c)
714 {
715 assert(c);
716 init(c->m, c->n+1);
717 for(int i=0; i<c->m; i++) {
718 (*this)[i][0] = (c->constant)[i];
719 for(int j=0; j<c->n; j++)
720 (*this)[i][j+1] = *((c->vars)+i*(c->n)+j);
721 }
722
723 }
724
725
726 void integer_matrix::init(FILE * fp)
727 {
728 int r, c;
729
730 fscanf(fp, "%d %d\n", &r, &c);
731 fprintf(stderr,"(%d x %d)\n", r, c);
732
733 init(fp, r, c);
734 }
735
736
737 void integer_matrix::init(FILE * fp, int r, int c)
738 {
739 init(r, c);
740 assert((r >= 0)&&(c >= 0));
741
742 for(int i=0; i<r; i++)
743 for(int j=0; j<c; j++)
744 fscanf(fp, "%d ", &(*this)[i][j]);
745 }
746
747 void integer_matrix::init(const int * data, int rows, int cols)
748 {
749 assert((rows >= 0)&&(cols >= 0));
750
751 for(int i=0; i<rows; i++)
752 for(int j=0; j<cols; j++)
753 A[i][j] = (*this)[i][j];
754
755 if (A) {
756 if((rows<=realmm)&&(cols<=n())) {
757 mm = rows; nn = cols;
758 for(int i=0; i<mm; i++)
759 A[i].init(data+i*cols,cols);
760 return;
761 } else
762 clear();
763 }
764
765 if(rows > 0) {
766 A = alloc_rows(rows);
767 for(int i=0; i<rows; i++)
768 A[i].init(data+i*cols,cols);
769 }
770 nn=cols;
771 mm=rows;
772 }
773
774
775
776 /********************************************************************************
777 * operators *
778 ********************************************************************************/
779 integer_matrix & integer_matrix::operator=(const integer_matrix & in)
780 {
781 if(this == &in)
782 return *this;
783 init(in);
784 return *this;
785 }
786
787
788 boolean integer_matrix::operator==(const integer_matrix & a) const
789 {
790 if((a.n() == n()) && (a.m() == m())) {
791 int noeq = 0;
792 for(int i=0; i<a.m(); i++) if(a.r(i) != r(i)) noeq = 1;
793 if(!noeq) return 1;
794 }
795 return 0;
796 }
797
798 integer_matrix integer_matrix::operator%(const integer_row & shuffle) const
799 {
800 assert(n() == shuffle.n());
801
802 integer_matrix result(m(), n());
803
804 for(int i=0; i<m(); i++) {
805 result[i] = r(i) % shuffle;
806 }
807
808 return result;
809 }
810
811
812
813 /********************************************************************************
814 * delete column *
815 ********************************************************************************/
816 integer_matrix integer_matrix::del_col(int i, int j) const
817 {
818 assert(i<= j);
819 assert((i>=0)&&(j<n()));
820 integer_matrix ret(m(), n() - (j - i + 1));
821
822 for(int a=0; a<m(); a++)
823 ret[a] = r(a).del_col(i, j);
824
825 return ret;
826 }
827
828 integer_matrix integer_matrix::insert_col(int i) const
829 {
830 assert((i>=0)&&(i<=n()));
831 integer_matrix ret(m(), n() + 1);
832
833 for(int a=0; a<m(); a++)
834 ret[a] = r(a).insert_col(i);
835
836 return ret;
837 }
838
839 integer_matrix integer_matrix::del_columns(const integer_row & mask) const
840 {
841 int cnt = 0;
842 int i, j;
843 assert(mask.n() == n());
844 for(j=0; j<mask.n(); j++)
845 if(!mask.c(j)) cnt++;
846
847 integer_matrix ret(m(), cnt);
848
849 int c = 0;
850 for(j=0; j<n(); j++)
851 if(!mask.c(j)) {
852 for(i=0; i<m(); i++)
853 ret[i][c] = this->r(i).c(j);
854 c++;
855 }
856 assert(c == cnt);
857
858 return ret;
859 }
860
861
862 integer_matrix integer_matrix::swap_col(int i, int j) const
863 {
864 assert((i>=0)&&(i<n()));
865 assert((j>=0)&&(j<n()));
866
867 integer_matrix ret(this);
868
869 for(int a=0; a<m(); a++)
870 ret[a].do_swap_col(i, j);
871
872 return ret;
873 }
874
875 /* in-place versions */
876
877 void integer_matrix::do_del_col(int i, int j)
878 {
879 assert(i<= j);
880 assert((i>=0)&&(j<n()));
881
882 for(int a=0; a<m(); a++)
883 A[a].do_del_col(i, j);
884
885 nn = nn - (j - i + 1);
886 }
887
888 void integer_matrix::do_insert_col(int i)
889 {
890 assert((i>=0)&&(i<=n()));
891
892 for(int a=0; a<m(); a++)
893 A[a].do_insert_col(i);
894 nn = nn+1;
895 }
896
897 void integer_matrix::do_del_columns(const integer_row & mask)
898 {
899 int cnt = 0;
900 int j;
901 assert(mask.n() == n());
902
903 for(j=0; j<mask.n(); j++)
904 if(!mask.c(j)) cnt++;
905
906 for(j=0; j<m(); j++)
907 A[j].do_del_columns(mask);
908 nn = cnt;
909 }
910
911 void integer_matrix::do_swap_col(int i, int j)
912 {
913 assert((i>=0)&&(i<n()));
914 assert((j>=0)&&(j<n()));
915
916 for(int a=0; a<m(); a++)
917 A[a].do_swap_col(i, j);
918 }
919
920 void integer_matrix::do_swap_row(int i, int j)
921 {
922 assert((i>=0)&&(i<m()));
923 assert((j>=0)&&(j<m()));
924
925 A[i] ^= A[j];
926 }
927
928 void integer_matrix::do_add_rows(int rws)
929 {
930 int newmm = mm+rws;
931 int i;
932 if (newmm > realmm) {
933 integer_row *A1 = alloc_rows(MAX(newmm,mm*2));
934 for (i=0; i<mm; i++) {
935 A[i] ^= A1[i];
936 }
937 delete[] A;
938 A = A1;
939 }
940 for (i=mm; i<newmm; i++) {
941 A[i].init(nn);
942 }
943 mm = newmm;
944 }
945
946 /********************************************************************************
947 * delete row *
948 ********************************************************************************/
949 integer_matrix integer_matrix::del_row(int i, int j) const
950 {
951 assert(i<= j);
952 assert((i>=0)&&(j<m()));
953 integer_matrix ret(m() - (j - i + 1), n());
954
955 int cnt = 0;
956 int a;
957 for(a = 0; a < i; a++)
958 ret[cnt++] = A[a];
959
960 for(a = j+1; a < m(); a++)
961 ret[cnt++] = A[a];
962
963 assert(cnt == ret.m());
964
965 return ret;
966 }
967
968 void integer_matrix::do_del_row(int i, int j)
969 {
970 assert(i<= j);
971 assert((i>=0)&&(j<m()));
972
973 int src, dst;
974
975 for (dst=i, src=j+1; src < m(); dst++, src++)
976 A[dst] ^= A[src];
977
978 mm = m() - (j - i + 1);
979 }
980
981 void integer_matrix::do_del_rows(const integer_row & mask)
982 {
983 int src, dst;
984 assert(mask.n() == m());
985
986 dst = 0;
987 while (dst < m() && !mask.R[dst])
988 dst++;
989 src = dst;
990 while (src < m()) {
991 if (!mask.R[src]) {
992 if (src != dst)
993 A[dst] ^= A[src];
994 dst++;
995 };
996 src++;
997 }
998
999 mm = dst;
1000 }
1001
1002 void integer_matrix::do_rowadd(int rto, int radd, int mul)
1003 {
1004 assert((radd>=0)&&(radd<m()));
1005 assert((rto>=0)&&(rto<m()));
1006 A[rto].do_addmul(A[radd],mul);
1007 }
1008
1009 void integer_matrix::do_coladd(int cto, int cadd, int mul)
1010 {
1011 assert((cto>=0)&&(cto<n()));
1012 assert((cadd>=0)&&(cadd<n()));
1013 if (mul == 1) {
1014 for (int i=0; i<m(); i++)
1015 A[i][cto] += A[i][cadd];
1016 } else {
1017 for (int i=0; i<m(); i++)
1018 A[i][cto] += A[i][cadd] * mul;
1019 }
1020 }
1021
1022
1023 /********************************************************************************
1024 * output *
1025 ********************************************************************************/
1026 void integer_matrix::print(FILE *fp, int flag) const
1027 {
1028 for (int i=0; i<mm; i++) {
1029 if(flag == 1) fprintf(fp,"%2d: ", i);
1030 for (int j=0; j<nn; j++) {
1031 fprintf(fp," %6d",r(i).c(j));
1032 }
1033 fprintf(fp,"\n");
1034 }
1035 }
1036
1037
1038
1039 void integer_matrix::check() const
1040 {
1041
1042 if((n() == 0)||(m() == 0)) return;
1043
1044 assert(A);
1045
1046 for(int i=0; i<m(); i++) {
1047 assert(A[i].n() == n());
1048 assert(A[i].R);
1049 }
1050
1051 }
1052
1053
1054 integer_matrix integer_matrix::operator+(const integer_matrix & a) const
1055 {
1056 assert(mm == a.mm);
1057 assert(nn == a.nn);
1058
1059 integer_matrix res(mm, nn);
1060
1061 for(int i=0; i<m(); i++)
1062 res.A[i] = A[i] + a.A[i];
1063
1064 return res;
1065 }
1066
1067
1068 integer_matrix integer_matrix::operator-(const integer_matrix & a) const
1069 {
1070 assert(mm == a.mm);
1071 assert(nn == a.nn);
1072
1073 integer_matrix res(mm, nn);
1074
1075 for(int i=0; i<m(); i++)
1076 res.A[i] = A[i] - a.A[i];
1077
1078 return res;
1079 }
1080
1081
1082 integer_matrix integer_matrix::operator*(const integer_matrix & a) const
1083 {
1084 assert(nn == a.mm);
1085
1086
1087 integer_matrix res(mm, a.nn);
1088
1089 for(int i=0; i<m(); i++)
1090 for(int j=0; j<a.n(); j++) {
1091 int acc = 0;
1092 for(int k=0; k<n(); k++)
1093 acc += r(i).c(k) * a.r(k).c(j);
1094 res.A[i][j] = acc;
1095 }
1096
1097 return res;
1098 }
1099
1100
1101
1102 integer_matrix integer_matrix::operator*(int val) const
1103 {
1104 integer_matrix res(this);
1105
1106 for(int i=0; i<res.m(); i++)
1107 for(int j=0; j<res.n(); j++) {
1108 res.A[i][j] *= val;
1109 }
1110
1111 return res;
1112 }
1113
1114 integer_matrix integer_matrix::operator/(int val) const
1115 {
1116 integer_matrix res(this);
1117
1118 for(int i=0; i<res.m(); i++)
1119 for(int j=0; j<res.n(); j++) {
1120 res.A[i][j] /= val;
1121 }
1122
1123 return res;
1124 }
1125
1126 integer_matrix integer_matrix::operator+(int val) const
1127 {
1128 integer_matrix res(this);
1129
1130 for(int i=0; i<res.m(); i++)
1131 for(int j=0; j<res.n(); j++) {
1132 res.A[i][j] += val;
1133 }
1134
1135 return res;
1136 }
1137
1138 integer_matrix integer_matrix::operator-(int val) const
1139 {
1140 integer_matrix res(this);
1141
1142 for(int i=0; i<res.m(); i++)
1143 for(int j=0; j<res.n(); j++) {
1144 res.A[i][j] -= val;
1145 }
1146
1147 return res;
1148 }
1149
1150 integer_matrix & integer_matrix::operator+=(const integer_matrix & in)
1151 {
1152 assert(m() == in.m());
1153 assert(n() == in.n());
1154
1155 for(int i=0; i<m(); i++)
1156 for(int j=0; j<n(); j++)
1157 (*this)[i][j] += in.r(i).c(j);
1158 return *this;
1159 }
1160
1161
1162 integer_matrix & integer_matrix::operator-=(const integer_matrix & in)
1163 {
1164 assert(m() == in.m());
1165 assert(n() == in.n());
1166
1167 for(int i=0; i<m(); i++)
1168 for(int j=0; j<n(); j++)
1169 (*this)[i][j] -= in.r(i).c(j);
1170 return *this;
1171 }
1172
1173
1174 integer_matrix & integer_matrix::operator+=(int x)
1175 {
1176 for(int i=0; i<m(); i++)
1177 for(int j=0; j<n(); j++)
1178 (*this)[i][j] += x;
1179 return *this;
1180 }
1181
1182 integer_matrix & integer_matrix::operator-=(int x)
1183 {
1184 for(int i=0; i<m(); i++)
1185 for(int j=0; j<n(); j++)
1186 (*this)[i][j] -= x;
1187 return *this;
1188 }
1189
1190 integer_matrix & integer_matrix::operator*=(int x)
1191 {
1192 for(int i=0; i<m(); i++)
1193 for(int j=0; j<n(); j++)
1194 (*this)[i][j] *= x;
1195 return *this;
1196 }
1197
1198 integer_matrix & integer_matrix::operator/=(int x)
1199 {
1200 for(int i=0; i<m(); i++)
1201 for(int j=0; j<n(); j++)
1202 (*this)[i][j] /= x;
1203 return *this;
1204 }
1205
1206 integer_matrix integer_matrix::transpose() const
1207 {
1208 integer_matrix res(nn, mm);
1209
1210 for(int i=0; i<m(); i++)
1211 for(int j=0; j<n(); j++)
1212 res.A[j][i] = A[i][j];
1213
1214 return res;
1215 }
1216
1217 int integer_matrix::determinant() const
1218 {
1219 assert(n() == m());
1220
1221 if(n() == 1) return A[0][0];
1222
1223 integer_matrix rdel;
1224 rdel = this->del_row(0);
1225
1226 int ret = 0;
1227 for(int sign=1, i=0; i<n();i++, sign *= -1)
1228 ret += sign*A[0][i]*rdel.del_col(i).determinant();
1229
1230 return ret;
1231 }
1232
1233
1234 int lcm(int, int);
1235 int gcd(int, int);
1236
1237
1238 integer_matrix integer_matrix::inverse(int * det) const
1239 {
1240 assert(nn == mm);
1241 integer_matrix curr(this);
1242
1243 integer_matrix res;
1244 res.ident(curr.m());
1245
1246 integer_matrix sw;
1247 int i;
1248 for(i=0; i<m(); i++) {
1249 if(curr[i][i] == 0) {
1250 for(int j=i+1; j<m(); j++)
1251 if(curr[j][i] != 0) {
1252 curr.do_swap_row(i,j);
1253 res.do_swap_row(i,j);
1254 break;
1255 }
1256 }
1257 assert(curr[i][i]);
1258 for(int j=0; j<m(); j++) {
1259 if((i != j)&&(curr[j][i])){
1260 int g = gcd(ABS(curr[i][i]), ABS(curr[j][i]));
1261 curr[j] *= ABS(curr[i][i])/g;
1262 res[j] *= ABS(curr[i][i])/g;
1263
1264 int mul = -curr[j][i]/curr[i][i];
1265 curr.do_rowadd(j, i, mul);
1266 res.do_rowadd(j, i, mul);
1267 }
1268 }
1269
1270 }
1271
1272 if (det) {
1273 int l = 1;
1274 for(i=0; i<m(); i++) {
1275 l = lcm(l, ABS(curr[i][i]));
1276 if(l > 10000) {
1277 fprintf(stderr, "Error lcm getting too large\n");
1278 this->print(stderr);
1279 return res;
1280 }
1281 }
1282
1283 *det = l;
1284
1285 for(i=0; i<m(); i++) {
1286 int norm = curr[i][i];
1287 curr[i] /= norm;
1288 res[i] *= l/norm;
1289 }
1290 } else {
1291 for(i=0; i<m(); i++) {
1292 int norm = curr[i][i];
1293 curr[i] /= norm;
1294 res[i] /= norm;
1295 }
1296 }
1297
1298 return res;
1299 }
1300
1301 void integer_matrix::ident(int v)
1302 {
1303 mknew(v, v);
1304
1305 for(int i=0; i<m(); i++) {
1306 A[i] = 0;
1307 A[i][i] = 1;
1308 }
1309 }
1310
1311
1312 integer_matrix integer_matrix::resize_offset(int r1, int r2, int c1, int c2,
1313 int fill) const
1314 {
1315
1316 assert(mm-r1+r2 >=0);
1317 assert(nn-c1+c2 >=0);
1318
1319 // if((mm-r1+r2 ==0)||(nn-c1+c2 ==0)) {
1320 // integer_matrix res;
1321 // return res;
1322 // }
1323
1324 integer_matrix res(mm-r1+r2, nn-c1+c2);
1325
1326 int i;
1327 for(i=0; i<res.m(); i++) res.A[i] = fill;
1328
1329 for(i=MAX(0,r1); i < mm+MIN(0,r2); i++)
1330 for(int j=MAX(0,c1); j < nn+MIN(0,c2); j++)
1331 res.A[i-r1][j-c1] = A[i][j];
1332
1333 return res;
1334 }
1335
1336
1337 integer_matrix integer_matrix::r_merge(const integer_matrix & a1,
1338 const integer_matrix & a2) const
1339 {
1340 assert(a1.n() == a2.n());
1341 integer_matrix res(a1.m()+a2.m(),0);
1342 int i,j;
1343 for (i=0; i<a1.m(); i++)
1344 res.A[i] = a1.A[i];
1345 for (j=0; j<a2.m(); j++,i++)
1346 res.A[i] = a2.A[j];
1347 res.nn=a1.n();
1348 return res;
1349 }
1350
1351 integer_matrix integer_matrix::c_merge(const integer_matrix & a1,
1352 const integer_matrix & a2) const
1353 {
1354 assert(a1.m() == a2.m());
1355 integer_matrix res(a1.m(),a1.n()+a2.n());
1356
1357 int off1 = a1.n()*sizeof(int);
1358 int off2 = a2.n()*sizeof(int);
1359 for (int i=0; i<a1.m(); i++) {
1360 memcpy(res.A[i].R, a1.A[i].R, off1);
1361 memcpy(res.A[i].R+a1.n(), a2.A[i].R, off2);
1362 }
1363
1364 return res;
1365 }
1366
1367
1368 // Permutation matricies
1369
1370 integer_matrix doswitch(int i, int x1, int x2)
1371 {
1372 integer_matrix res;
1373 res.ident(i);
1374
1375 res.A[x1][x1] = 0;
1376 res.A[x2][x2] = 0;
1377 res.A[x1][x2] = 1;
1378 res.A[x2][x1] = 1;
1379
1380 return res;
1381 }
1382
1383 integer_matrix rowswitch(const integer_matrix & A, int r1, int r2)
1384 {
1385 assert((r1>=0)&&(r1<A.m()));
1386 assert((r2>=0)&&(r2<A.m()));
1387
1388 return doswitch(A.m(), r1, r2);
1389 }
1390
1391
1392 integer_matrix colswitch(const integer_matrix & A, int c1, int c2)
1393 {
1394 assert((c1>=0)&&(c1<A.n()));
1395 assert((c2>=0)&&(c2<A.n()));
1396
1397 return doswitch(A.n(), c1, c2);
1398 }
1399
1400 integer_matrix doadd(int i, int r, int c, int val)
1401 {
1402 integer_matrix res;
1403 res.ident(i);
1404
1405 assert(r!=c);
1406
1407 res[r][c] = val;
1408
1409 return res;
1410 }
1411
1412 integer_matrix rowadd(const integer_matrix & A, int r1, int r2, int val)
1413 {
1414 assert((r1>=0)&&(r1<A.m()));
1415 assert((r2>=0)&&(r2<A.m()));
1416
1417 return doadd(A.m(), r1, r2, val);
1418 }
1419
1420 integer_matrix coladd(const integer_matrix & A, int c1, int c2, int val)
1421 {
1422 assert((c1>=0)&&(c1<A.n()));
1423 assert((c2>=0)&&(c2<A.n()));
1424
1425 return doadd(A.n(), c2, c1, val);
1426 }
1427
1428 void compose_print(int fg)
1429 {
1430 compose_print_flag = fg;
1431 }
1432
1433 integer_matrix Compose(int r_comp, int c_comp,
1434 const integer_matrix * m1, const integer_matrix * m2,
1435 const integer_matrix * m3, const integer_matrix * m4,
1436 const integer_matrix * m5, const integer_matrix * m6,
1437 const integer_matrix * m7, const integer_matrix * m8,
1438 const integer_matrix * m9, const integer_matrix * m10,
1439 const integer_matrix * m11, const integer_matrix * m12,
1440 const integer_matrix * m13, const integer_matrix * m14,
1441 const integer_matrix * m15, const integer_matrix * m16,
1442 const integer_matrix * m17, const integer_matrix * m18,
1443 const integer_matrix * m19, const integer_matrix * m20,
1444 const integer_matrix * m21, const integer_matrix * m22,
1445 const integer_matrix * m23, const integer_matrix * m24,
1446 const integer_matrix * m25, const integer_matrix * m26,
1447 const integer_matrix * m27, const integer_matrix * m28,
1448 const integer_matrix * m29, const integer_matrix * m30,
1449 const integer_matrix * m31, const integer_matrix * m32)
1450 {
1451 integer_matrix R;
1452
1453 return R.compose(r_comp, c_comp,
1454 m1, m2, m3, m4, m5, m6, m7, m8,
1455 m9, m10, m11, m12, m13, m14, m15, m16,
1456 m17, m18, m19, m20, m21, m22, m23, m24,
1457 m25, m26, m27, m28, m29, m30, m31, m32);
1458 }
1459
1460 integer_matrix integer_matrix::compose(int r_comp, int c_comp,
1461 const integer_matrix * m1,
1462 const integer_matrix * m2,
1463 const integer_matrix * m3,
1464 const integer_matrix * m4,
1465 const integer_matrix * m5,
1466 const integer_matrix * m6,
1467 const integer_matrix * m7,
1468 const integer_matrix * m8,
1469 const integer_matrix * m9,
1470 const integer_matrix * m10,
1471 const integer_matrix * m11,
1472 const integer_matrix * m12,
1473 const integer_matrix * m13,
1474 const integer_matrix * m14,
1475 const integer_matrix * m15,
1476 const integer_matrix * m16,
1477 const integer_matrix * m17,
1478 const integer_matrix * m18,
1479 const integer_matrix * m19,
1480 const integer_matrix * m20,
1481 const integer_matrix * m21,
1482 const integer_matrix * m22,
1483 const integer_matrix * m23,
1484 const integer_matrix * m24,
1485 const integer_matrix * m25,
1486 const integer_matrix * m26,
1487 const integer_matrix * m27,
1488 const integer_matrix * m28,
1489 const integer_matrix * m29,
1490 const integer_matrix * m30,
1491 const integer_matrix * m31,
1492 const integer_matrix * m32) const
1493 {
1494 assert(r_comp > 0);
1495 assert(c_comp > 0);
1496 int val = r_comp*c_comp;
1497 assert(val <= 32);
1498 const integer_matrix **list = new const integer_matrix *[val];
1499
1500 list[0] = m1;
1501 if(val>1) list[1] = m2;
1502 if(val>2) list[2] = m3;
1503 if(val>3) list[3] = m4;
1504 if(val>4) list[4] = m5;
1505 if(val>5) list[5] = m6;
1506 if(val>6) list[6] = m7;
1507 if(val>7) list[7] = m8;
1508 if(val>8) list[8] = m9;
1509 if(val>9) list[9] = m10;
1510 if(val>10) list[10] = m11;
1511 if(val>11) list[11] = m12;
1512 if(val>12) list[12] = m13;
1513 if(val>13) list[13] = m14;
1514 if(val>14) list[14] = m15;
1515 if(val>15) list[15] = m16;
1516 if(val>16) list[16] = m17;
1517 if(val>17) list[17] = m18;
1518 if(val>18) list[18] = m19;
1519 if(val>19) list[19] = m20;
1520 if(val>20) list[20] = m21;
1521 if(val>21) list[21] = m22;
1522 if(val>22) list[22] = m23;
1523 if(val>23) list[23] = m24;
1524 if(val>24) list[24] = m25;
1525 if(val>25) list[25] = m26;
1526 if(val>26) list[26] = m27;
1527 if(val>27) list[27] = m28;
1528 if(val>28) list[28] = m29;
1529 if(val>29) list[29] = m30;
1530 if(val>30) list[30] = m31;
1531 if(val>31) list[31] = m32;
1532
1533 integer_matrix ret = composeL(r_comp, c_comp, list);
1534
1535 delete[] list;
1536 return ret;
1537 }
1538
1539
1540 // DECLARE_DLIST_CLASSES(nim_list_base, nim_e,
1541 // nim_list_iter, integer_matrix *);
1542 //
1543 // class nim_list: public nim_list_base {
1544 // public:
1545 // void clear();
1546 // };
1547 //
1548 // void nim_list::clear()
1549 // {
1550 // nim_list_iter iter(this);
1551 // while(!iter.is_empty()) {
1552 // integer_matrix * m = iter.step();
1553 // delete m;
1554 // }
1555 // glist::clear();
1556 // }
1557 //
1558 // nim_list NIM_List;
1559
1560 nim_op::~nim_op()
1561 {
1562 if(m) delete m;
1563 if(next) delete next;
1564 }
1565
1566 integer_matrix * nim_op::NIM(const integer_matrix & A)
1567 {
1568 nim_op * no = new nim_op(this);
1569 no->m = new integer_matrix(A);
1570 return no->m;
1571 }
1572
1573 integer_matrix * nim_op::NIM(const integer_matrix * A)
1574 {
1575 assert(A);
1576 nim_op * no = new nim_op(this);
1577 no->m = new integer_matrix(A);
1578 return no->m;
1579 }
1580
1581 integer_matrix * nim_op::NIM(int i, int j)
1582 {
1583 nim_op * no = new nim_op(this);
1584 no->m = new integer_matrix(i, j);
1585 return no->m;
1586 }
1587
1588 integer_matrix * nim_op::NIM(int i)
1589 {
1590 nim_op * no = new nim_op(this);
1591 no->m = new integer_matrix(1, 1);
1592 (no->m)[0][0] = i;
1593 return no->m;
1594 }
1595
1596 integer_matrix * NIM(const integer_matrix & A)
1597 {
1598 integer_matrix * m = new integer_matrix(A);
1599 return m;
1600 }
1601
1602
1603 integer_matrix * NIM(const integer_matrix * A)
1604 {
1605 assert(A);
1606 integer_matrix * m = new integer_matrix(A);
1607 return m;
1608 }
1609
1610 integer_matrix * NIM(int i, int j)
1611 {
1612 integer_matrix * m = new integer_matrix(i, j);
1613 return m;
1614 }
1615
1616 integer_matrix * NIM(int i)
1617 {
1618 integer_matrix * m = new integer_matrix(1, 1);
1619 (*m)[0][0] = i;
1620 return m;
1621 }
1622
1623
1624
1625 #define ELEM(i, j) list[i*c_comp+j]
1626
1627 integer_matrix integer_matrix::composeL(int r_comp, int c_comp,
1628 const integer_matrix *const * list,
1629 int exit_on_error) const
1630 {
1631 assert(r_comp > 0);
1632 assert(c_comp > 0);
1633 // int val = r_comp*c_comp;
1634
1635 int * r_seg_sz = new int[r_comp];
1636 int * c_seg_sz = new int[c_comp];
1637 int rsz, csz;
1638 int i, j;
1639
1640 rsz = 0;
1641 for(i=0; i<r_comp; i++) {
1642 int rseg = -1;
1643 for(j=0; j<c_comp; j++)
1644 if(ELEM(i,j) != 0) {
1645 if(rseg == -1) rseg = ELEM(i,j)->m();
1646 if(rseg != ELEM(i,j)->m()) {
1647 composeL_print(TRUE, r_comp, c_comp, list);
1648 if(exit_on_error)
1649 error_line(1, NULL, "Row has unmached row count: element [%d, %d] row count %d, should be %d\n", i, j, ELEM(i,j)->m(), rseg);
1650 else return integer_matrix(0, 0);
1651 }
1652
1653
1654 }
1655 if(rseg == -1) {
1656 composeL_print(TRUE, r_comp, c_comp, list);
1657 if(exit_on_error)
1658 error_line(1, NULL, "Row %d has no row count\n", i);
1659 else return integer_matrix(0, 0);
1660 }
1661 r_seg_sz[i] = rseg;
1662 rsz += rseg;
1663 }
1664
1665 csz = 0;
1666 for(j=0; j<c_comp; j++) {
1667 int cseg = -1;
1668 for(i=0; i<r_comp; i++)
1669 if(ELEM(i,j) != 0) {
1670 if(cseg == -1) cseg = ELEM(i,j)->n();
1671 if(cseg != ELEM(i, j)->n()) {
1672 composeL_print(TRUE, r_comp, c_comp, list);
1673 if(exit_on_error)
1674 error_line(1, NULL,"Column has unmatched column count: element [%d, %d] column count %d, should be %d\n", i, j, ELEM(i,j)->n(), cseg);
1675 else return integer_matrix(0, 0);
1676 }
1677 }
1678 if(cseg == -1) {
1679 composeL_print(TRUE, r_comp, c_comp, list);
1680 if(exit_on_error)
1681 error_line(1, NULL, "Column %d has no column count \n", j);
1682 else return integer_matrix(0, 0);
1683 }
1684 c_seg_sz[j] = cseg;
1685 csz += cseg;
1686 }
1687
1688 integer_matrix ret(rsz, csz);
1689
1690 int rorg = 0;
1691 for(i=0; i<r_comp; i++) {
1692 int corg = 0;
1693 for(j=0; j<c_comp; j++) {
1694 const integer_matrix * M;
1695 if((M = ELEM(i, j)) != 0) {
1696 for(int ii = 0; ii < r_seg_sz[i]; ii++)
1697 for(int jj = 0; jj < c_seg_sz[j]; jj++)
1698 ret[ii+rorg][jj+corg] = M->r(ii).c(jj);
1699 }
1700 corg += c_seg_sz[j];
1701 }
1702 rorg += r_seg_sz[i];
1703 }
1704
1705
1706 delete r_seg_sz;
1707 delete c_seg_sz;
1708
1709 if(compose_print_flag)
1710 composeL_print(FALSE, r_comp, c_comp, list);
1711 return ret;
1712 }
1713
1714 void integer_matrix::composeL_print(int er, int r_comp, int c_comp,
1715 const integer_matrix *const * list) const
1716 {
1717 int i, j;
1718
1719 if(er) fprintf(stderr, "Composing error:\n");
1720 for(i=0; i< r_comp; i++) {
1721 for(j=0; j< c_comp; j++)
1722 if(ELEM(i,j) == 0)
1723 fprintf(stderr, " ----- ");
1724 else
1725 fprintf(stderr, "[%2dx%2d] ", ELEM(i,j)->m(), ELEM(i,j)->n());
1726 fprintf(stderr, "\n");
1727 }
1728 fprintf(stderr, "\n");
1729 }
1730
1731
1732 int integer_matrix::init(const immed_list & ilist, int c)
1733 {
1734 int xm = ilist[c++].integer();
1735 int xn = ilist[c++].integer();
1736 init(xm, xn);
1737 for(int i=0; i<m(); i++)
1738 for(int j=0; j<n(); j++)
1739 (*this)[i][j] = ilist[c++].integer();
1740
1741 return c;
1742 }
1743
1744
1745 immed_list * integer_matrix::cvt_immed_list() const
1746 {
1747 immed_list * il = new immed_list;
1748 il->append(immed(m()));
1749 il->append(immed(n()));
1750 for(int i=0; i<m(); i++)
1751 for(int j=0; j<n(); j++)
1752 il->append(immed(r(i).c(j)));
1753 return il;
1754 }
1755
1756
1757 integer_matrix Compose(int r_comp, int c_comp,
1758 int a1, int a2, int a3, int a4, int a5, int a6,
1759 int a7, int a8, int a9, int a10, int a11, int a12,
1760 int a13, int a14, int a15, int a16, int a17, int a18,
1761 int a19, int a20, int a21, int a22, int a23, int a24,
1762 int a25, int a26, int a27, int a28, int a29, int a30,
1763 int a31, int a32, int a33, int a34, int a35, int a36,
1764 int a37, int a38, int a39, int a40, int a41, int a42,
1765 int a43, int a44, int a45, int a46, int a47, int a48,
1766 int a49, int a50, int a51, int a52, int a53, int a54,
1767 int a55, int a56, int a57, int a58, int a59, int a60)
1768 {
1769 int params[60];
1770
1771 params[0] = a1;
1772 params[1] = a2;
1773 params[2] = a3;
1774 params[3] = a4;
1775 params[4] = a5;
1776 params[5] = a6;
1777 params[6] = a7;
1778 params[7] = a8;
1779 params[8] = a9;
1780 params[9] = a10;
1781 params[10] = a11;
1782 params[11] = a12;
1783 params[12] = a13;
1784 params[13] = a14;
1785 params[14] = a15;
1786 params[15] = a16;
1787 params[16] = a17;
1788 params[17] = a18;
1789 params[18] = a19;
1790 params[19] = a20;
1791 params[20] = a21;
1792 params[21] = a22;
1793 params[22] = a23;
1794 params[23] = a24;
1795 params[24] = a25;
1796 params[25] = a26;
1797 params[26] = a27;
1798 params[27] = a28;
1799 params[28] = a29;
1800 params[29] = a30;
1801 params[30] = a31;
1802 params[31] = a32;
1803 params[32] = a33;
1804 params[33] = a34;
1805 params[34] = a35;
1806 params[35] = a36;
1807 params[36] = a37;
1808 params[37] = a38;
1809 params[38] = a39;
1810 params[39] = a40;
1811 params[40] = a41;
1812 params[41] = a42;
1813 params[42] = a43;
1814 params[43] = a44;
1815 params[44] = a45;
1816 params[45] = a46;
1817 params[46] = a47;
1818 params[47] = a48;
1819 params[48] = a49;
1820 params[49] = a50;
1821 params[50] = a51;
1822 params[51] = a52;
1823 params[52] = a53;
1824 params[53] = a54;
1825 params[54] = a55;
1826 params[55] = a56;
1827 params[56] = a57;
1828 params[57] = a58;
1829 params[58] = a59;
1830 params[59] = a60;
1831
1832 assert(r_comp*c_comp <= 60);
1833 integer_matrix M(r_comp, c_comp);
1834 int * ptr = params;
1835 for(int i=0; i<r_comp; i++)
1836 for(int j=0; j<c_comp; j++)
1837 M[i][j] = *ptr++;
1838 return M;
1839 }
1840
1841
1842