My Project
Loading...
Searching...
No Matches
bigintmat.cc
Go to the documentation of this file.
1/*****************************************
2 * Computer Algebra System SINGULAR *
3 *****************************************/
4/*
5 * ABSTRACT: class bigintmat: matrices of numbers.
6 * a few functinos might be limited to bigint or euclidean rings.
7 */
8
9
10#include "misc/auxiliary.h"
11
12#include "coeffs/bigintmat.h"
13#include "misc/intvec.h"
14
15#include "coeffs/rmodulon.h"
16
17#include <cmath>
18
19#ifdef HAVE_RINGS
20///create Z/nA of type n_Zn
21static coeffs numbercoeffs(number n, coeffs c) // TODO: FIXME: replace with n_CoeffRingQuot1
22{
23 mpz_t p;
24 number2mpz(n, c, p);
25 ZnmInfo *pp = new ZnmInfo;
26 pp->base = p;
27 pp->exp = 1;
28 coeffs nc = nInitChar(n_Zn, (void*)pp);
29 mpz_clear(p);
30 delete pp;
31 return nc;
32}
33#endif
34
35//#define BIMATELEM(M,I,J) (M)[ (M).index(I,J) ]
36
38{
39 bigintmat * t = new bigintmat(col, row, basecoeffs());
40 for (int i=1; i<=row; i++)
41 {
42 for (int j=1; j<=col; j++)
43 {
44 t->set(j, i, BIMATELEM(*this,i,j));
45 }
46 }
47 return t;
48}
49
51{
52 int n = row,
53 m = col,
54 nm = n<m?n : m; // the min, describing the square part of the matrix
55 //CF: this is not optimal, but so far, it seems to work
56
57#define swap(_i, _j) \
58 int __i = (_i), __j=(_j); \
59 number c = v[__i]; \
60 v[__i] = v[__j]; \
61 v[__j] = c \
62
63 for (int i=0; i< nm; i++)
64 for (int j=i+1; j< nm; j++)
65 {
66 swap(i*m+j, j*n+i);
67 }
68 if (n<m)
69 for (int i=nm; i<m; i++)
70 for(int j=0; j<n; j++)
71 {
72 swap(j*n+i, i*m+j);
73 }
74 if (n>m)
75 for (int i=nm; i<n; i++)
76 for(int j=0; j<m; j++)
77 {
78 swap(i*m+j, j*n+i);
79 }
80#undef swap
81 row = m;
82 col = n;
83}
84
85
86// Beginnt bei [0]
87void bigintmat::set(int i, number n, const coeffs C)
88{
89 assume (C == NULL || C == basecoeffs());
90
92}
93
94// Beginnt bei [1,1]
95void bigintmat::set(int i, int j, number n, const coeffs C)
96{
97 assume (C == NULL || C == basecoeffs());
98 assume (i > 0 && j > 0);
99 assume (i <= rows() && j <= cols());
100 set(index(i, j), n, C);
101}
102
104{
105 assume (i >= 0);
106 assume (i<rows()*cols());
107
108 return n_Copy(v[i], basecoeffs());
109}
110
112{
113 assume (i >= 0);
114 assume (i<rows()*cols());
115
116 return v[i];
117}
118
119number bigintmat::get(int i, int j) const
120{
121 assume (i > 0 && j > 0);
122 assume (i <= rows() && j <= cols());
123
124 return get(index(i, j));
125}
126
127number bigintmat::view(int i, int j) const
128{
129 assume (i >= 0 && j >= 0);
130 assume (i <= rows() && j <= cols());
131
132 return view(index(i, j));
133}
134// Ueberladener *=-Operator (für int und bigint)
135// Frage hier: *= verwenden oder lieber = und * einzeln?
137{
139
141
143}
144
146{
147 assume (C == NULL || C == basecoeffs());
148
149 const int l = rows() * cols();
150
151 for (int i=0; i < l; i++)
153}
154
155// Stimmen Parameter?
156// Welche der beiden Methoden?
157// Oder lieber eine comp-Funktion?
158
159bool operator==(const bigintmat & lhr, const bigintmat & rhr)
160{
161 if (&lhr == &rhr) { return true; }
162 if (lhr.cols() != rhr.cols()) { return false; }
163 if (lhr.rows() != rhr.rows()) { return false; }
164 if (lhr.basecoeffs() != rhr.basecoeffs()) { return false; }
165
166 const int l = (lhr.rows())*(lhr.cols());
167
168 for (int i=0; i < l; i++)
169 {
170 if (!n_Equal(lhr[i], rhr[i], lhr.basecoeffs())) { return false; }
171 }
172
173 return true;
174}
175
176bool operator!=(const bigintmat & lhr, const bigintmat & rhr)
177{
178 return !(lhr==rhr);
179}
180
181// Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ?
183{
184 if (a->cols() != b->cols()) return NULL;
185 if (a->rows() != b->rows()) return NULL;
186 if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
187
188 const coeffs basecoeffs = a->basecoeffs();
189
190 int i;
191
192 bigintmat * bim = new bigintmat(a->rows(), a->cols(), basecoeffs);
193
194 for (i=a->rows()*a->cols()-1;i>=0; i--)
195 bim->rawset(i, n_Add((*a)[i], (*b)[i], basecoeffs), basecoeffs);
196
197 return bim;
198}
200{
201
202 const int mn = si_min(a->rows(),a->cols());
203
204 const coeffs basecoeffs = a->basecoeffs();
205 number bb=n_Init(b,basecoeffs);
206
207 int i;
208
209 bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
210
211 for (i=1; i<=mn; i++)
212 BIMATELEM(*bim,i,i)=n_Add(BIMATELEM(*a,i,i), bb, basecoeffs);
213
214 n_Delete(&bb,basecoeffs);
215 return bim;
216}
217
219{
220 if (a->cols() != b->cols()) return NULL;
221 if (a->rows() != b->rows()) return NULL;
222 if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
223
224 const coeffs basecoeffs = a->basecoeffs();
225
226 int i;
227
228 bigintmat * bim = new bigintmat(a->rows(), a->cols(), basecoeffs);
229
230 for (i=a->rows()*a->cols()-1;i>=0; i--)
231 bim->rawset(i, n_Sub((*a)[i], (*b)[i], basecoeffs), basecoeffs);
232
233 return bim;
234}
235
237{
238 const int mn = si_min(a->rows(),a->cols());
239
240 const coeffs basecoeffs = a->basecoeffs();
241 number bb=n_Init(b,basecoeffs);
242
243 int i;
244
245 bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
246
247 for (i=1; i<=mn; i++)
248 BIMATELEM(*bim,i,i)=n_Sub(BIMATELEM(*a,i,i), bb, basecoeffs);
249
250 n_Delete(&bb,basecoeffs);
251 return bim;
252}
253
254//TODO: make special versions for certain rings!
256{
257 const int ca = a->cols();
258 const int cb = b->cols();
259
260 const int ra = a->rows();
261 const int rb = b->rows();
262
263 if (ca != rb)
264 {
265#ifndef SING_NDEBUG
266 Werror("wrong bigintmat sizes at multiplication a * b: acols: %d != brows: %d\n", ca, rb);
267#endif
268 return NULL;
269 }
270
271 assume (ca == rb);
272
273 if (a->basecoeffs() != b->basecoeffs()) { return NULL; }
274
275 const coeffs basecoeffs = a->basecoeffs();
276
277 int i, j, k;
278
279 number sum;
280
281 bigintmat * bim = new bigintmat(ra, cb, basecoeffs);
282
283 for (i=1; i<=ra; i++)
284 for (j=1; j<=cb; j++)
285 {
286 sum = n_Init(0, basecoeffs);
287
288 for (k=1; k<=ca; k++)
289 {
290 number prod = n_Mult( BIMATELEM(*a, i, k), BIMATELEM(*b, k, j), basecoeffs);
291
292 n_InpAdd(sum, prod, basecoeffs);
293
294 n_Delete(&prod, basecoeffs);
295 }
296 bim->rawset(i, j, sum, basecoeffs);
297 }
298 return bim;
299}
300
302{
303
304 const int mn = a->rows()*a->cols();
305
306 const coeffs basecoeffs = a->basecoeffs();
307 number bb=n_Init(b,basecoeffs);
308
309 int i;
310
311 bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
312
313 for (i=0; i<mn; i++)
314 bim->rawset(i, n_Mult((*a)[i], bb, basecoeffs), basecoeffs);
315
316 n_Delete(&bb,basecoeffs);
317 return bim;
318}
319
321{
322 if (cf!=a->basecoeffs()) return NULL;
323
324 const int mn = a->rows()*a->cols();
325
326 const coeffs basecoeffs = a->basecoeffs();
327
328 int i;
329
330 bigintmat * bim = new bigintmat(a->rows(),a->cols() , basecoeffs);
331
332 for (i=0; i<mn; i++)
333 bim->rawset(i, n_Mult((*a)[i], b, basecoeffs), basecoeffs);
334
335 return bim;
336}
337
338// ----------------------------------------------------------------- //
339// Korrekt?
340
342{
343 intvec * iv = new intvec(b->rows(), b->cols(), 0);
344 for (int i=0; i<(b->rows())*(b->cols()); i++)
345 (*iv)[i] = n_Int((*b)[i], b->basecoeffs()); // Geht das so?
346 return iv;
347}
348
350{
351 const int l = (b->rows())*(b->cols());
352 bigintmat * bim = new bigintmat(b->rows(), b->cols(), C);
353
354 for (int i=0; i < l; i++)
355 bim->rawset(i, n_Init((*b)[i], C), C);
356
357 return bim;
358}
359
360// ----------------------------------------------------------------- //
361
362int bigintmat::compare(const bigintmat* op) const
363{
364 assume (basecoeffs() == op->basecoeffs() );
365
366#ifndef SING_NDEBUG
367 if (basecoeffs() != op->basecoeffs() )
368 WerrorS("wrong bigintmat comparison: different basecoeffs!\n");
369#endif
370
371 if ((col!=1) ||(op->cols()!=1))
372 {
373 if((col!=op->cols())
374 || (row!=op->rows()))
375 return -2;
376 }
377
378 int i;
379 for (i=0; i<si_min(row*col,op->rows()*op->cols()); i++)
380 {
381 if ( n_Greater(v[i], (*op)[i], basecoeffs()) )
382 return 1;
383 else if (! n_Equal(v[i], (*op)[i], basecoeffs()))
384 return -1;
385 }
386
387 for (; i<row; i++)
388 {
389 if ( n_GreaterZero(v[i], basecoeffs()) )
390 return 1;
391 else if (! n_IsZero(v[i], basecoeffs()) )
392 return -1;
393 }
394 for (; i<op->rows(); i++)
395 {
396 if ( n_GreaterZero((*op)[i], basecoeffs()) )
397 return -1;
398 else if (! n_IsZero((*op)[i], basecoeffs()) )
399 return 1;
400 }
401 return 0;
402}
403
404
406{
407 if (b == NULL)
408 return NULL;
409
410 return new bigintmat(b);
411}
412
414{
415 int n = cols(), m=rows();
416
417 for(int i=1; i<= m; i++)
418 {
419 for(int j=1; j< n; j++)
420 {
421 n_Write(v[(i-1)*n+j-1], basecoeffs());
422 StringAppendS(", ");
423 }
424 if (n) n_Write(v[i*n-1], basecoeffs());
425 if (i<m)
426 {
427 StringAppendS(", ");
428 }
429 }
430}
431
433{
434 StringSetS("");
435 Write();
436 return StringEndS();
437}
438
440{
441 char * s = String();
442 PrintS(s);
443 omFree(s);
444}
445
446
448{
449 if ((col==0) || (row==0))
450 return NULL;
451 else
452 {
453 int * colwid = getwid(80);
454 if (colwid == NULL)
455 {
456 WerrorS("not enough space to print bigintmat");
457 WerrorS("try string(...) for a unformatted output");
458 return NULL;
459 }
460 char * ps;
461 int slength = 0;
462 for (int j=0; j<col; j++)
463 slength += colwid[j]*row;
464 slength += col*row+row;
465 ps = (char*) omAlloc0(sizeof(char)*(slength));
466 int pos = 0;
467 for (int i=0; i<col*row; i++)
468 {
469 StringSetS("");
470 n_Write(v[i], basecoeffs());
471 char * ts = StringEndS();
472 const int _nl = strlen(ts);
473 int cj = i%col;
474 if (_nl > colwid[cj])
475 {
476 StringSetS("");
477 int ci = i/col;
478 StringAppend("[%d,%d]", ci+1, cj+1);
479 char * ph = StringEndS();
480 int phl = strlen(ph);
481 if (phl > colwid[cj])
482 {
483 for (int j=0; j<colwid[cj]-1; j++)
484 ps[pos+j] = ' ';
485 ps[pos+colwid[cj]-1] = '*';
486 }
487 else
488 {
489 for (int j=0; j<colwid[cj]-phl; j++)
490 ps[pos+j] = ' ';
491 for (int j=0; j<phl; j++)
492 ps[pos+colwid[cj]-phl+j] = ph[j];
493 }
494 omFree(ph);
495 }
496 else // Mit Leerzeichen auffüllen und zahl reinschreiben
497 {
498 for (int j=0; j<(colwid[cj]-_nl); j++)
499 ps[pos+j] = ' ';
500 for (int j=0; j<_nl; j++)
501 ps[pos+colwid[cj]-_nl+j] = ts[j];
502 }
503 // ", " und (evtl) "\n" einfügen
504 if ((i+1)%col == 0)
505 {
506 if (i != col*row-1)
507 {
508 ps[pos+colwid[cj]] = ',';
509 ps[pos+colwid[cj]+1] = '\n';
510 pos += colwid[cj]+2;
511 }
512 }
513 else
514 {
515 ps[pos+colwid[cj]] = ',';
516 pos += colwid[cj]+1;
517 }
518 omFree(ts); // Hier ts zerstören
519 }
520 return(ps);
521 // omFree(ps);
522}
523}
524
525static int intArrSum(int * a, int length)
526{
527 int sum = 0;
528 for (int i=0; i<length; i++)
529 sum += a[i];
530 return sum;
531}
532
533static int findLongest(int * a, int length)
534{
535 int l = 0;
536 int index;
537 for (int i=0; i<length; i++)
538 {
539 if (a[i] > l)
540 {
541 l = a[i];
542 index = i;
543 }
544 }
545 return index;
546}
547
548static int getShorter (int * a, int l, int j, int cols, int rows)
549{
550 int sndlong = 0;
551 int min;
552 for (int i=0; i<rows; i++)
553 {
554 int index = cols*i+j;
555 if ((a[index] > sndlong) && (a[index] < l))
556 {
557 min = floor(log10((double)cols))+floor(log10((double)rows))+5;
558 if ((a[index] < min) && (min < l))
559 sndlong = min;
560 else
561 sndlong = a[index];
562 }
563 }
564 if (sndlong == 0)
565 {
566 min = floor(log10((double)cols))+floor(log10((double)rows))+5;
567 if (min < l)
568 sndlong = min;
569 else
570 sndlong = 1;
571 }
572 return sndlong;
573}
574
575
577{
578 int const c = /*2**/(col-1)+1;
579 int * wv = (int*)omAlloc(sizeof(int)*col*row);
580 int * cwv = (int*)omAlloc(sizeof(int)*col);
581 for (int j=0; j<col; j++)
582 {
583 cwv[j] = 0;
584 for (int i=0; i<row; i++)
585 {
586 StringSetS("");
587 n_Write(v[col*i+j], basecoeffs());
588 char * tmp = StringEndS();
589 const int _nl = strlen(tmp);
590 wv[col*i+j] = _nl;
591 if (_nl > cwv[j]) cwv[j]=_nl;
592 omFree(tmp);
593 }
594 }
595
596 // Groesse verkleinern, bis < maxwid
597 if (intArrSum(cwv, col)+c > maxwid)
598 {
599 int j = findLongest(cwv, col);
600 cwv[j] = getShorter(wv, cwv[j], j, col, row);
601 }
602 omFree(wv);
603 return cwv;
604}
605
607{
608 if ((col==0) || (row==0))
609 PrintS("");
610 else
611 {
612 int * colwid = getwid(maxwid);
613 char * ps;
614 int slength = 0;
615 for (int j=0; j<col; j++)
616 slength += colwid[j]*row;
617 slength += col*row+row;
618 ps = (char*) omAlloc0(sizeof(char)*(slength));
619 int pos = 0;
620 for (int i=0; i<col*row; i++)
621 {
622 StringSetS("");
623 n_Write(v[i], basecoeffs());
624 char * ts = StringEndS();
625 const int _nl = strlen(ts);
626 int cj = i%col;
627 if (_nl > colwid[cj])
628 {
629 StringSetS("");
630 int ci = i/col;
631 StringAppend("[%d,%d]", ci+1, cj+1);
632 char * ph = StringEndS();
633 int phl = strlen(ph);
634 if (phl > colwid[cj])
635 {
636 for (int j=0; j<colwid[cj]-1; j++)
637 ps[pos+j] = ' ';
638 ps[pos+colwid[cj]-1] = '*';
639 }
640 else
641 {
642 for (int j=0; j<colwid[cj]-phl; j++)
643 ps[pos+j] = ' ';
644 for (int j=0; j<phl; j++)
645 ps[pos+colwid[cj]-phl+j] = ph[j];
646 }
647 omFree(ph);
648 }
649 else // Mit Leerzeichen auffüllen und zahl reinschreiben
650 {
651 for (int j=0; j<colwid[cj]-_nl; j++)
652 ps[pos+j] = ' ';
653 for (int j=0; j<_nl; j++)
654 ps[pos+colwid[cj]-_nl+j] = ts[j];
655 }
656 // ", " und (evtl) "\n" einfügen
657 if ((i+1)%col == 0)
658 {
659 if (i != col*row-1)
660 {
661 ps[pos+colwid[cj]] = ',';
662 ps[pos+colwid[cj]+1] = '\n';
663 pos += colwid[cj]+2;
664 }
665 }
666 else
667 {
668 ps[pos+colwid[cj]] = ',';
669 pos += colwid[cj]+1;
670 }
671
672 omFree(ts); // Hier ts zerstören
673 }
674 PrintS(ps);
675 omFree(ps);
676 }
677}
678
679
680//swaps columns i and j
681void bigintmat::swap(int i, int j)
682{
683 if ((i <= col) && (j <= col) && (i>0) && (j>0))
684 {
685 number tmp;
686 number t;
687 for (int k=1; k<=row; k++)
688 {
689 tmp = get(k, i);
690 t = view(k, j);
691 set(k, i, t);
692 set(k, j, tmp);
694 }
695 }
696 else
697 WerrorS("Error in swap");
698}
699
700void bigintmat::swaprow(int i, int j)
701{
702 if ((i <= row) && (j <= row) && (i>0) && (j>0))
703 {
704 number tmp;
705 number t;
706 for (int k=1; k<=col; k++)
707 {
708 tmp = get(i, k);
709 t = view(j, k);
710 set(i, k, t);
711 set(j, k, tmp);
713 }
714 }
715 else
716 WerrorS("Error in swaprow");
717}
718
720{
721 for (int j=1; j<=col; j++)
722 {
723 if (!n_IsZero(view(i,j), basecoeffs()))
724 {
725 return j;
726 }
727 }
728 return 0;
729}
730
732{
733 for (int i=row; i>=1; i--)
734 {
735 if (!n_IsZero(view(i,j), basecoeffs()))
736 {
737 return i;
738 }
739 }
740 return 0;
741}
742
744{
745 assume((j<=col) && (j>=1));
746 if (((a->rows() != row) || (a->cols() != 1)) && ((a->rows() != 1) || (a->cols() != row)))
747 {
748 assume(0);
749 WerrorS("Error in getcol. Dimensions must agree!");
750 return;
751 }
753 {
755 number t1, t2;
756 for (int i=1; i<=row;i++)
757 {
758 t1 = get(i,j);
759 t2 = f(t1, basecoeffs(), a->basecoeffs());
760 a->set(i-1,t1);
761 n_Delete(&t1, basecoeffs());
762 n_Delete(&t2, a->basecoeffs());
763 }
764 return;
765 }
766 number t1;
767 for (int i=1; i<=row;i++)
768 {
769 t1 = view(i,j);
770 a->set(i-1,t1);
771 }
772}
773
775{
776 number t1;
777 for(int ii=0; ii< no; ii++)
778 {
779 for (int i=1; i<=row;i++)
780 {
781 t1 = view(i, ii+j);
782 a->set(i, ii+1, t1);
783 }
784 }
785}
786
788{
789 if ((i>row) || (i<1))
790 {
791 WerrorS("Error in getrow: Index out of range!");
792 return;
793 }
794 if (((a->rows() != 1) || (a->cols() != col)) && ((a->rows() != col) || (a->cols() != 1)))
795 {
796 WerrorS("Error in getrow. Dimensions must agree!");
797 return;
798 }
800 {
802 number t1, t2;
803 for (int j=1; j<=col;j++)
804 {
805 t1 = get(i,j);
806 t2 = f(t1, basecoeffs(), a->basecoeffs());
807 a->set(j-1,t2);
808 n_Delete(&t1, basecoeffs());
809 n_Delete(&t2, a->basecoeffs());
810 }
811 return;
812 }
813 number t1;
814 for (int j=1; j<=col;j++)
815 {
816 t1 = get(i,j);
817 a->set(j-1,t1);
818 n_Delete(&t1, basecoeffs());
819 }
820}
821
823{
824 if ((j>col) || (j<1))
825 {
826 WerrorS("Error in setcol: Index out of range!");
827 return;
828 }
829 if (((m->rows() != row) || (m->cols() != 1)) && ((m->rows() != 1) || (m->cols() != row)))
830 {
831 WerrorS("Error in setcol. Dimensions must agree!");
832 return;
833 }
834 if (!nCoeffs_are_equal(basecoeffs(), m->basecoeffs()))
835 {
836 nMapFunc f = n_SetMap(m->basecoeffs(), basecoeffs());
837 number t1,t2;
838 for (int i=1; i<=row; i++)
839 {
840 t1 = m->get(i-1);
841 t2 = f(t1, m->basecoeffs(),basecoeffs());
842 set(i, j, t2);
843 n_Delete(&t2, basecoeffs());
844 n_Delete(&t1, m->basecoeffs());
845 }
846 return;
847 }
848 number t1;
849 for (int i=1; i<=row; i++)
850 {
851 t1 = m->view(i-1);
852 set(i, j, t1);
853 }
854}
855
857{
858 if ((j>row) || (j<1))
859 {
860 WerrorS("Error in setrow: Index out of range!");
861 return;
862 }
863 if (((m->rows() != 1) || (m->cols() != col)) && ((m->rows() != col) || (m->cols() != 1)))
864 {
865 WerrorS("Error in setrow. Dimensions must agree!");
866 return;
867 }
868 if (!nCoeffs_are_equal(basecoeffs(), m->basecoeffs()))
869 {
870 nMapFunc f = n_SetMap(m->basecoeffs(), basecoeffs());
872 for (int i=1; i<=col; i++)
873 {
874 tmp1 = m->get(i-1);
875 tmp2 = f(tmp1, m->basecoeffs(),basecoeffs());
876 set(j, i, tmp2);
878 n_Delete(&tmp1, m->basecoeffs());
879 }
880 return;
881 }
882 number tmp;
883 for (int i=1; i<=col; i++)
884 {
885 tmp = m->view(i-1);
886 set(j, i, tmp);
887 }
888}
889
891{
892 if ((b->rows() != row) || (b->cols() != col))
893 {
894 WerrorS("Error in bigintmat::add. Dimensions do not agree!");
895 return false;
896 }
897 if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
898 {
899 WerrorS("Error in bigintmat::add. coeffs do not agree!");
900 return false;
901 }
902 for (int i=1; i<=row; i++)
903 {
904 for (int j=1; j<=col; j++)
905 {
906 rawset(i, j, n_Add(b->view(i,j), view(i,j), basecoeffs()));
907 }
908 }
909 return true;
910}
911
913{
914 if ((b->rows() != row) || (b->cols() != col))
915 {
916 WerrorS("Error in bigintmat::sub. Dimensions do not agree!");
917 return false;
918 }
919 if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
920 {
921 WerrorS("Error in bigintmat::sub. coeffs do not agree!");
922 return false;
923 }
924 for (int i=1; i<=row; i++)
925 {
926 for (int j=1; j<=col; j++)
927 {
928 rawset(i, j, n_Sub(view(i,j), b->view(i,j), basecoeffs()));
929 }
930 }
931 return true;
932}
933
935{
936 if (!nCoeffs_are_equal(c, basecoeffs()))
937 {
938 WerrorS("Wrong coeffs\n");
939 return false;
940 }
941 number t1, t2;
942 if ( n_IsOne(b,c)) return true;
943 for (int i=1; i<=row; i++)
944 {
945 for (int j=1; j<=col; j++)
946 {
947 t1 = view(i, j);
948 t2 = n_Mult(t1, b, basecoeffs());
949 rawset(i, j, t2);
950 }
951 }
952 return true;
953}
954
955bool bigintmat::addcol(int i, int j, number a, coeffs c)
956{
957 if ((i>col) || (j>col) || (i<1) || (j<1))
958 {
959 WerrorS("Error in addcol: Index out of range!");
960 return false;
961 }
962 if (!nCoeffs_are_equal(c, basecoeffs()))
963 {
964 WerrorS("Error in addcol: coeffs do not agree!");
965 return false;
966 }
967 number t1, t2, t3;
968 for (int k=1; k<=row; k++)
969 {
970 t1 = view(k, j);
971 t2 = view(k, i);
972 t3 = n_Mult(t1, a, basecoeffs());
973 n_InpAdd(t3, t2, basecoeffs());
974 rawset(k, i, t3);
975 }
976 return true;
977}
978
979bool bigintmat::addrow(int i, int j, number a, coeffs c)
980{
981 if ((i>row) || (j>row) || (i<1) || (j<1))
982 {
983 WerrorS("Error in addrow: Index out of range!");
984 return false;
985 }
986 if (!nCoeffs_are_equal(c, basecoeffs()))
987 {
988 WerrorS("Error in addrow: coeffs do not agree!");
989 return false;
990 }
991 number t1, t2, t3;
992 for (int k=1; k<=col; k++)
993 {
994 t1 = view(j, k);
995 t2 = view(i, k);
996 t3 = n_Mult(t1, a, basecoeffs());
997 n_InpAdd(t3, t2, basecoeffs());
998 rawset(i, k, t3);
999 }
1000 return true;
1001}
1002
1004{
1005 if ((i>=1) && (i<=col) && (nCoeffs_are_equal(c, basecoeffs())))
1006 {
1007 number t, tmult;
1008 for (int j=1; j<=row; j++)
1009 {
1010 t = view(j, i);
1011 tmult = n_Mult(a, t, basecoeffs());
1012 rawset(j, i, tmult);
1013 }
1014 }
1015 else
1016 WerrorS("Error in colskalmult");
1017}
1018
1020{
1021 if ((i>=1) && (i<=row) && (nCoeffs_are_equal(c, basecoeffs())))
1022 {
1023 number t, tmult;
1024 for (int j=1; j<=col; j++)
1025 {
1026 t = view(i, j);
1027 tmult = n_Mult(a, t, basecoeffs());
1028 rawset(i, j, tmult);
1029 }
1030 }
1031 else
1032 WerrorS("Error in rowskalmult");
1033}
1034
1036{
1037 int ay = a->cols();
1038 int ax = a->rows();
1039 int by = b->cols();
1040 int bx = b->rows();
1041 number tmp;
1042 if (!((col == ay) && (col == by) && (ax+bx == row)))
1043 {
1044 WerrorS("Error in concatrow. Dimensions must agree!");
1045 return;
1046 }
1047 if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
1048 {
1049 WerrorS("Error in concatrow. coeffs do not agree!");
1050 return;
1051 }
1052 for (int i=1; i<=ax; i++)
1053 {
1054 for (int j=1; j<=ay; j++)
1055 {
1056 tmp = a->get(i,j);
1057 set(i, j, tmp);
1058 n_Delete(&tmp, basecoeffs());
1059 }
1060 }
1061 for (int i=1; i<=bx; i++)
1062 {
1063 for (int j=1; j<=by; j++)
1064 {
1065 tmp = b->get(i,j);
1066 set(i+ax, j, tmp);
1067 n_Delete(&tmp, basecoeffs());
1068 }
1069 }
1070}
1071
1073{
1074 bigintmat * tmp = new bigintmat(rows(), i, basecoeffs());
1075 appendCol(tmp);
1076 delete tmp;
1077}
1078
1080{
1081 coeffs R = basecoeffs();
1082 int ay = a->cols();
1083 int ax = a->rows();
1084 assume(row == ax);
1085
1087
1088 bigintmat * tmp = new bigintmat(rows(), cols() + ay, R);
1089 tmp->concatcol(this, a);
1090 this->swapMatrix(tmp);
1091 delete tmp;
1092}
1093
1095 int ay = a->cols();
1096 int ax = a->rows();
1097 int by = b->cols();
1098 int bx = b->rows();
1099 number tmp;
1100
1101 assume(row==ax && row == bx && ay+by ==col);
1102
1104
1105 for (int i=1; i<=ax; i++)
1106 {
1107 for (int j=1; j<=ay; j++)
1108 {
1109 tmp = a->view(i,j);
1110 set(i, j, tmp);
1111 }
1112 }
1113 for (int i=1; i<=bx; i++)
1114 {
1115 for (int j=1; j<=by; j++)
1116 {
1117 tmp = b->view(i,j);
1118 set(i, j+ay, tmp);
1119 }
1120 }
1121}
1122
1124{
1125 int ay = a->cols();
1126 int ax = a->rows();
1127 int by = b->cols();
1128 int bx = b->rows();
1129 number tmp;
1130 if (!(ax + bx == row))
1131 {
1132 WerrorS("Error in splitrow. Dimensions must agree!");
1133 }
1134 else if (!((col == ay) && (col == by)))
1135 {
1136 WerrorS("Error in splitrow. Dimensions must agree!");
1137 }
1138 else if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
1139 {
1140 WerrorS("Error in splitrow. coeffs do not agree!");
1141 }
1142 else
1143 {
1144 for(int i = 1; i<=ax; i++)
1145 {
1146 for(int j = 1; j<=ay;j++)
1147 {
1148 tmp = get(i,j);
1149 a->set(i,j,tmp);
1150 n_Delete(&tmp, basecoeffs());
1151 }
1152 }
1153 for (int i =1; i<=bx; i++)
1154 {
1155 for (int j=1;j<=col;j++)
1156 {
1157 tmp = get(i+ax, j);
1158 b->set(i,j,tmp);
1159 n_Delete(&tmp, basecoeffs());
1160 }
1161 }
1162 }
1163}
1164
1166{
1167 int ay = a->cols();
1168 int ax = a->rows();
1169 int by = b->cols();
1170 int bx = b->rows();
1171 number tmp;
1172 if (!((row == ax) && (row == bx)))
1173 {
1174 WerrorS("Error in splitcol. Dimensions must agree!");
1175 }
1176 else if (!(ay+by == col))
1177 {
1178 WerrorS("Error in splitcol. Dimensions must agree!");
1179 }
1180 else if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
1181 {
1182 WerrorS("Error in splitcol. coeffs do not agree!");
1183 }
1184 else
1185 {
1186 for (int i=1; i<=ax; i++)
1187 {
1188 for (int j=1; j<=ay; j++)
1189 {
1190 tmp = view(i,j);
1191 a->set(i,j,tmp);
1192 }
1193 }
1194 for (int i=1; i<=bx; i++)
1195 {
1196 for (int j=1; j<=by; j++)
1197 {
1198 tmp = view(i,j+ay);
1199 b->set(i,j,tmp);
1200 }
1201 }
1202 }
1203}
1204
1206{
1207 number tmp;
1208 if ((a->rows() != row) || (a->cols()+i-1 > col) || (i<1))
1209 {
1210 WerrorS("Error in splitcol. Dimensions must agree!");
1211 return;
1212 }
1213 if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs())))
1214 {
1215 WerrorS("Error in splitcol. coeffs do not agree!");
1216 return;
1217 }
1218 int width = a->cols();
1219 for (int j=1; j<=width; j++)
1220 {
1221 for (int k=1; k<=row; k++)
1222 {
1223 tmp = get(k, j+i-1);
1224 a->set(k, j, tmp);
1225 n_Delete(&tmp, basecoeffs());
1226 }
1227 }
1228}
1229
1231{
1232 number tmp;
1233 if ((a->cols() != col) || (a->rows()+i-1 > row) || (i<1))
1234 {
1235 WerrorS("Error in Marco-splitrow");
1236 return;
1237 }
1238
1239 if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs())))
1240 {
1241 WerrorS("Error in splitrow. coeffs do not agree!");
1242 return;
1243 }
1244 int height = a->rows();
1245 for (int j=1; j<=height; j++)
1246 {
1247 for (int k=1; k<=col; k++)
1248 {
1249 tmp = view(j+i-1, k);
1250 a->set(j, k, tmp);
1251 }
1252 }
1253}
1254
1256{
1257 if ((b->rows() != row) || (b->cols() != col))
1258 {
1259 WerrorS("Error in bigintmat::copy. Dimensions do not agree!");
1260 return false;
1261 }
1262 if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
1263 {
1264 WerrorS("Error in bigintmat::copy. coeffs do not agree!");
1265 return false;
1266 }
1267 number t1;
1268 for (int i=1; i<=row; i++)
1269 {
1270 for (int j=1; j<=col; j++)
1271 {
1272 t1 = b->view(i, j);
1273 set(i, j, t1);
1274 }
1275 }
1276 return true;
1277}
1278
1279/// copy the submatrix of b, staring at (a,b) having n rows, m cols into
1280/// the given matrix at pos. (c,d)
1281/// needs c+n, d+m <= rows, cols
1282/// a+n, b+m <= b.rows(), b.cols()
1283void bigintmat::copySubmatInto(bigintmat *B, int a, int b, int n, int m, int c, int d)
1284{
1285 number t1;
1286 for (int i=1; i<=n; i++)
1287 {
1288 for (int j=1; j<=m; j++)
1289 {
1290 t1 = B->view(a+i-1, b+j-1);
1291 set(c+i-1, d+j-1, t1);
1292 }
1293 }
1294}
1295
1297{
1298 coeffs r = basecoeffs();
1299 if (row==col)
1300 {
1301 for (int i=1; i<=row; i++)
1302 {
1303 for (int j=1; j<=col; j++)
1304 {
1305 if (i==j)
1306 {
1307 if (!n_IsOne(view(i, j), r))
1308 return 0;
1309 }
1310 else
1311 {
1312 if (!n_IsZero(view(i,j), r))
1313 return 0;
1314 }
1315 }
1316 }
1317 }
1318 return 1;
1319}
1320
1322{
1323 if (row==col)
1324 {
1325 number one = n_Init(1, basecoeffs()),
1326 zero = n_Init(0, basecoeffs());
1327 for (int i=1; i<=row; i++)
1328 {
1329 for (int j=1; j<=col; j++)
1330 {
1331 if (i==j)
1332 {
1333 set(i, j, one);
1334 }
1335 else
1336 {
1337 set(i, j, zero);
1338 }
1339 }
1340 }
1341 n_Delete(&one, basecoeffs());
1343 }
1344}
1345
1347{
1348 number tmp = n_Init(0, basecoeffs());
1349 for (int i=1; i<=row; i++)
1350 {
1351 for (int j=1; j<=col; j++)
1352 {
1353 set(i, j, tmp);
1354 }
1355 }
1357}
1358
1360{
1361 for (int i=1; i<=row; i++) {
1362 for (int j=1; j<=col; j++) {
1363 if (!n_IsZero(view(i,j), basecoeffs()))
1364 return FALSE;
1365 }
1366 }
1367 return TRUE;
1368}
1369//****************************************************************************
1370//
1371//****************************************************************************
1372
1373//used in the det function. No idea what it does.
1374//looks like it return the submatrix where the i-th row
1375//and j-th column has been removed in the LaPlace generic
1376//determinant algorithm
1378{
1379 if ((i<=0) || (i>row) || (j<=0) || (j>col))
1380 return NULL;
1381 int cx, cy;
1382 cx=1;
1383 cy=1;
1384 number t;
1385 bigintmat *b = new bigintmat(row-1, col-1, basecoeffs());
1386 for (int k=1; k<=row; k++) {
1387 if (k!=i)
1388 {
1389 cy=1;
1390 for (int l=1; l<=col; l++)
1391 {
1392 if (l!=j)
1393 {
1394 t = get(k, l);
1395 b->set(cx, cy, t);
1396 n_Delete(&t, basecoeffs());
1397 cy++;
1398 }
1399 }
1400 cx++;
1401 }
1402 }
1403 return b;
1404}
1405
1406
1407//returns d such that a/d is the inverse of the input
1408//TODO: make work for p not prime using the euc stuff.
1409//long term: rewrite for Z using p-adic lifting
1410//and Dixon. Possibly even the sparse recent Storjohann stuff
1412
1413 // Falls Matrix über reellen Zahlen nicht invertierbar, breche ab
1414 assume((a->rows() == row) && (a->rows() == a->cols()) && (row == col));
1415
1416 number det = this->det(); //computes the HNF, so should e reused.
1417 if ((n_IsZero(det, basecoeffs())))
1418 return det;
1419
1420 // Hänge Einheitsmatrix über Matrix und wendet HNF an. An Stelle der Einheitsmatrix steht im Ergebnis die Transformationsmatrix dazu
1421 a->one();
1422 bigintmat *m = new bigintmat(2*row, col, basecoeffs());
1423 m->concatrow(a,this);
1424 m->hnf();
1425 // Arbeite weiterhin mit der zusammengehängten Matrix
1426 // Laufe durch die Diagonalelemente, und multipliziere jede Spalte rechts davon damit, speichere aber den alten Eintrag der Spalte, temp, der in der Zeile des Diagonalelements liegt, zwischen. Dann addiere das -temp-Fache der Diagonalspalte zur entsprechenenden Spalte rechts davon. Dadurch entsteht überall rechts der Diagonalen eine 0
1427 number diag;
1428 number temp, ttemp;
1429 for (int i=1; i<=col; i++) {
1430 diag = m->get(row+i, i);
1431 for (int j=i+1; j<=col; j++) {
1432 temp = m->get(row+i, j);
1433 m->colskalmult(j, diag, basecoeffs());
1435 m->addcol(j, i, temp, basecoeffs());
1437 }
1439 }
1440 // Falls wir nicht modulo n arbeiten, können wir die Spalten durch den ggT teilen, um die Einträge kleiner zu bekommen
1441 // Bei Z/n sparen wir uns das, da es hier sinnlos ist
1442 number g;
1443 number gcd;
1444 for (int j=1; j<=col; j++) {
1445 g = n_Init(0, basecoeffs());
1446 for (int i=1; i<=2*row; i++) {
1447 temp = m->get(i,j);
1448 gcd = n_Gcd(g, temp, basecoeffs());
1449 n_Delete(&g, basecoeffs());
1451 g = n_Copy(gcd, basecoeffs());
1452 n_Delete(&gcd, basecoeffs());
1453 }
1454 if (!(n_IsOne(g, basecoeffs())))
1455 m->colskaldiv(j, g);
1456 n_Delete(&g, basecoeffs());
1457 }
1458
1459 // Nun müssen die Diagonalelemente durch Spaltenmultiplikation gleich gesett werden. Bei Z können wir mit dem kgV arbeiten, bei Z/n bringen wir jedes Diagonalelement auf 1 (wir arbeiten immer mit n = Primzahl. Für n != Primzahl muss noch an anderen Stellen etwas geändert werden)
1460
1461 g = n_Init(0, basecoeffs());
1462 number prod = n_Init(1, basecoeffs());
1463 for (int i=1; i<=col; i++) {
1464 gcd = n_Gcd(g, m->get(row+i, i), basecoeffs());
1465 n_Delete(&g, basecoeffs());
1466 g = n_Copy(gcd, basecoeffs());
1467 n_Delete(&gcd, basecoeffs());
1469 temp = m->get(row+i, i);
1474 }
1475 number lcm;
1476 lcm = n_Div(prod, g, basecoeffs());
1477 for (int j=1; j<=col; j++) {
1478 ttemp = m->get(row+j,j);
1480 m->colskalmult(j, temp, basecoeffs());
1483 }
1484 n_Delete(&lcm, basecoeffs());
1486
1487 number divisor = m->get(row+1, 1);
1488 m->splitrow(a, 1);
1489 delete m;
1490 n_Delete(&det, basecoeffs());
1491 return divisor;
1492}
1493
1495{
1496 assume (col == row);
1497 number t = get(1,1),
1498 h;
1499 coeffs r = basecoeffs();
1500 for(int i=2; i<= col; i++) {
1501 h = n_Add(t, view(i,i), r);
1502 n_Delete(&t, r);
1503 t = h;
1504 }
1505 return t;
1506}
1507
1509{
1510 assume (row==col);
1511
1512 if (col == 1)
1513 return get(1, 1);
1514 // should work as well in Z/pZ of type n_Zp?
1515 // relies on XExtGcd and the other euc. functinos.
1517 return hnfdet();
1518 }
1519 number sum = n_Init(0, basecoeffs());
1520 number t1, t2, t3, t4;
1521 bigintmat *b;
1522 for (int i=1; i<=row; i++) {
1523 b = elim(i, 1);
1524 t1 = get(i, 1);
1525 t2 = b->det();
1526 t3 = n_Mult(t1, t2, basecoeffs());
1527 t4 = n_Copy(sum, basecoeffs());
1528 n_Delete(&sum, basecoeffs());
1529 if ((i+1)>>1<<1==(i+1))
1530 sum = n_Add(t4, t3, basecoeffs());
1531 else
1532 sum = n_Sub(t4, t3, basecoeffs());
1533 n_Delete(&t1, basecoeffs());
1534 n_Delete(&t2, basecoeffs());
1535 n_Delete(&t3, basecoeffs());
1536 n_Delete(&t4, basecoeffs());
1537 }
1538 return sum;
1539}
1540
1542{
1543 assume (col == row);
1544
1545 if (col == 1)
1546 return get(1, 1);
1547 bigintmat *m = new bigintmat(this);
1548 m->hnf();
1549 number prod = n_Init(1, basecoeffs());
1550 number temp, temp2;
1551 for (int i=1; i<=col; i++) {
1552 temp = m->get(i, i);
1555 prod = temp2;
1557 }
1558 delete m;
1559 return prod;
1560}
1561
1563{
1564 int n = rows(), m = cols();
1565 row = a->rows();
1566 col = a->cols();
1567 number * V = v;
1568 v = a->v;
1569 a->v = V;
1570 a->row = n;
1571 a->col = m;
1572}
1574{
1575 coeffs R = basecoeffs();
1576 for(int i=1; i<=rows(); i++)
1577 if (!n_IsZero(view(i, j), R)) return FALSE;
1578 return TRUE;
1579}
1580
1582{
1583 coeffs R = basecoeffs();
1584 hnf(); // as a starting point...
1585 if (getCoeffType(R)== n_Z) return; //wrong, need to prune!
1586
1587 int n = cols(), m = rows(), i, j, k;
1588
1589 //make sure, the matrix has enough space. We need no rows+1 columns.
1590 //The resulting Howell form will be pruned to be at most square.
1591 bigintmat * t = new bigintmat(m, m+1, R);
1592 t->copySubmatInto(this, 1, n>m ? n-m+1 : 1, m, n>m ? m : n, 1, n>m ? 2 : m+2-n );
1593 swapMatrix(t);
1594 delete t;
1595 for(i=1; i<= cols(); i++) {
1596 if (!colIsZero(i)) break;
1597 }
1598 assume (i>1);
1599 if (i>cols()) {
1600 t = new bigintmat(rows(), 0, R);
1601 swapMatrix(t);
1602 delete t;
1603 return; // zero matrix found, clearly normal.
1604 }
1605
1606 int last_zero_col = i-1;
1607 for (int c = cols(); c>0; c--) {
1608 for(i=rows(); i>0; i--) {
1609 if (!n_IsZero(view(i, c), R)) break;
1610 }
1611 if (i==0) break; // matrix SHOULD be zero from here on
1612 number a = n_Ann(view(i, c), R);
1613 addcol(last_zero_col, c, a, R);
1614 n_Delete(&a, R);
1615 for(j = c-1; j>last_zero_col; j--) {
1616 for(k=rows(); k>0; k--) {
1617 if (!n_IsZero(view(k, j), R)) break;
1618 if (!n_IsZero(view(k, last_zero_col), R)) break;
1619 }
1620 if (k==0) break;
1621 if (!n_IsZero(view(k, last_zero_col), R)) {
1622 number gcd, co1, co2, co3, co4;
1623 gcd = n_XExtGcd(view(k, last_zero_col), view(k, j), &co1, &co2, &co3, &co4, R);
1624 if (n_Equal(gcd, view(k, j), R)) {
1625 number q = n_Div(view(k, last_zero_col), gcd, R);
1626 q = n_InpNeg(q, R);
1627 addcol(last_zero_col, j, q, R);
1628 n_Delete(&q, R);
1629 } else if (n_Equal(gcd, view(k, last_zero_col), R)) {
1631 number q = n_Div(view(k, last_zero_col), gcd, R);
1632 q = n_InpNeg(q, R);
1633 addcol(last_zero_col, j, q, R);
1634 n_Delete(&q, R);
1635 } else {
1637 }
1638 n_Delete(&gcd, R);
1639 n_Delete(&co1, R);
1640 n_Delete(&co2, R);
1641 n_Delete(&co3, R);
1642 n_Delete(&co4, R);
1643 }
1644 }
1645 for(k=rows(); k>0; k--) {
1646 if (!n_IsZero(view(k, last_zero_col), R)) break;
1647 }
1648 if (k) last_zero_col--;
1649 }
1650 t = new bigintmat(rows(), cols()-last_zero_col, R);
1651 t->copySubmatInto(this, 1, last_zero_col+1, rows(), cols()-last_zero_col, 1, 1);
1652 swapMatrix(t);
1653 delete t;
1654}
1655
1657{
1658 // Laufen von unten nach oben und von links nach rechts
1659 // CF: TODO: for n_Z: write a recursive version. This one will
1660 // have exponential blow-up. Look at Michianchio
1661 // Alternatively, do p-adic det and modular method
1662
1663#if 0
1664 char * s;
1665 ::PrintS("mat over Z is \n");
1666 ::Print("%s\n", s = nCoeffString(basecoeffs()));
1667 omFree(s);
1668 Print();
1669 ::Print("\n(%d x %d)\n", rows(), cols());
1670#endif
1671
1672 int i = rows();
1673 int j = cols();
1674 number q = n_Init(0, basecoeffs());
1675 number one = n_Init(1, basecoeffs());
1677 number tmp1 = n_Init(0, basecoeffs());
1678 number tmp2 = n_Init(0, basecoeffs());
1679 number co1, co2, co3, co4;
1680 number ggt = n_Init(0, basecoeffs());
1681
1682 while ((i>0) && (j>0))
1683 {
1684 // Falls erstes Nicht-Null-Element in Zeile i nicht existiert, oder hinter Spalte j vorkommt, gehe in nächste Zeile
1685 if ((findnonzero(i)==0) || (findnonzero(i)>j))
1686 {
1687 i--;
1688 }
1689 else
1690 {
1691 // Laufe von links nach rechts durch die Zeile:
1692 for (int l=1; l<=j-1; l++)
1693 {
1695 tmp1 = get(i, l);
1696 // Falls Eintrag (im folgenden x genannt) gleich 0, gehe eine Spalte weiter. Ansonsten...
1697 if (!n_IsZero(tmp1, basecoeffs()))
1698 {
1700 tmp2 = get(i, l+1);
1701 // Falls Eintrag (i.f. y g.) rechts daneben gleich 0, tausche beide Spalten, sonst...
1702 if (!n_IsZero(tmp2, basecoeffs()))
1703 {
1704 n_Delete(&ggt, basecoeffs());
1705 ggt = n_XExtGcd(tmp1, tmp2, &co1, &co2, &co3, &co4, basecoeffs());
1706 // Falls x=ggT(x, y), tausche die beiden Spalten und ziehe die (neue) rechte Spalte so häufig von der linken ab, dass an der ehemaligen Stelle von x nun eine 0 steht. Dazu:
1707 if (n_Equal(tmp1, ggt, basecoeffs()))
1708 {
1709 swap(l, l+1);
1710 n_Delete(&q, basecoeffs());
1711 q = n_Div(tmp2, ggt, basecoeffs());
1712 q = n_InpNeg(q, basecoeffs());
1713 // Dann addiere das -q-fache der (neuen) rechten Spalte zur linken dazu. Damit erhalten wir die gewünschte 0
1714
1715 addcol(l, l+1, q, basecoeffs());
1716 n_Delete(&q, basecoeffs());
1717 }
1718 else if (n_Equal(tmp1, minusone, basecoeffs()))
1719 {
1720 // Falls x=-1, so ist x=-ggt(x, y). Dann gehe wie oben vor, multipliziere aber zuerst die neue rechte Spalte (die mit x) mit -1
1721 // Die Berechnung von q (=y/ggt) entfällt, da ggt=1
1722 swap(l, l+1);
1725 addcol(l, l+1, tmp2, basecoeffs());
1726 }
1727 else
1728 {
1729 // CF: use the 2x2 matrix (co1, co2)(co3, co4) to
1730 // get the gcd in position and the 0 in the other:
1731#ifdef CF_DEB
1732 ::PrintS("applying trafo\n");
1733 StringSetS("");
1738 ::Print("%s\nfor l=%d\n", StringEndS(), l);
1739 {char * s = String();
1740 ::Print("to %s\n", s);omFree(s);};
1741#endif
1742 coltransform(l, l+1, co3, co4, co1, co2);
1743#ifdef CF_DEB
1744 {char * s = String();
1745 ::Print("gives %s\n", s);}
1746#endif
1747 }
1748 n_Delete(&co1, basecoeffs());
1749 n_Delete(&co2, basecoeffs());
1750 n_Delete(&co3, basecoeffs());
1751 n_Delete(&co4, basecoeffs());
1752 }
1753 else
1754 {
1755 swap(l, l+1);
1756 }
1757 // Dann betrachte die vormals rechte Spalte als neue linke, und die rechts daneben als neue rechte.
1758 }
1759 }
1760
1761 #ifdef HAVE_RINGS
1762 // normalize by units:
1763 if (!n_IsZero(view(i, j), basecoeffs()))
1764 {
1765 number u = n_GetUnit(view(i, j), basecoeffs());
1766 if (!n_IsOne(u, basecoeffs()))
1767 {
1768 colskaldiv(j, u);
1769 }
1770 n_Delete(&u, basecoeffs());
1771 }
1772 #endif
1773 // Zum Schluss mache alle Einträge rechts vom Diagonalelement betragsmäßig kleiner als dieses
1774 for (int l=j+1; l<=col; l++)
1775 {
1776 n_Delete(&q, basecoeffs());
1777 q = n_QuotRem(view(i, l), view(i, j), NULL, basecoeffs());
1778 q = n_InpNeg(q, basecoeffs());
1779 addcol(l, j, q, basecoeffs());
1780 }
1781 i--;
1782 j--;
1783 // Dann betrachte die Zeile darüber und gehe dort wie vorher vor
1784 }
1785 }
1786 n_Delete(&q, basecoeffs());
1789 n_Delete(&ggt, basecoeffs());
1790 n_Delete(&one, basecoeffs());
1792
1793#if 0
1794 ::PrintS("hnf over Z is \n");
1795 Print();
1796 ::Print("\n(%d x %d)\n", rows(), cols());
1797#endif
1798}
1799
1801{
1802 coeffs cold = a->basecoeffs();
1803 bigintmat *b = new bigintmat(a->rows(), a->cols(), cnew);
1804 // Erzeugt Karte von alten coeffs nach neuen
1806 number t1;
1807 number t2;
1808 // apply map to all entries.
1809 for (int i=1; i<=a->rows(); i++)
1810 {
1811 for (int j=1; j<=a->cols(); j++)
1812 {
1813 t1 = a->get(i, j);
1814 t2 = f(t1, cold, cnew);
1815 b->set(i, j, t2);
1816 n_Delete(&t1, cold);
1817 n_Delete(&t2, cnew);
1818 }
1819 }
1820 return b;
1821}
1822
1823#ifdef HAVE_RINGS
1824//OK: a HNF of (this | p*I)
1825//so the result will always have FULL rank!!!!
1826//(This is different form a lift of the HNF mod p: consider the matrix (p)
1827//to see the difference. It CAN be computed as HNF mod p^2 usually..)
1829{
1830 coeffs Rp = numbercoeffs(p, R); // R/pR
1831 bigintmat *m = bimChangeCoeff(this, Rp);
1832 m->howell();
1833 bigintmat *a = bimChangeCoeff(m, R);
1834 delete m;
1835 bigintmat *C = new bigintmat(rows(), rows(), R);
1836 int piv = rows(), i = a->cols();
1837 while (piv)
1838 {
1839 if (!i || n_IsZero(a->view(piv, i), R))
1840 {
1841 C->set(piv, piv, p, R);
1842 }
1843 else
1844 {
1845 C->copySubmatInto(a, 1, i, rows(), 1, 1, piv);
1846 i--;
1847 }
1848 piv--;
1849 }
1850 delete a;
1851 return C;
1852}
1853#endif
1854
1855
1856//exactly divide matrix by b
1858{
1859 number tmp1, tmp2;
1860 for (int i=1; i<=row; i++)
1861 {
1862 for (int j=1; j<=col; j++)
1863 {
1864 tmp1 = view(i, j);
1865 tmp2 = n_Div(tmp1, b, basecoeffs());
1866 rawset(i, j, tmp2);
1867 }
1868 }
1869}
1870
1871//exactly divide col j by b
1873{
1874 number tmp1, tmp2;
1875 for (int i=1; i<=row; i++)
1876 {
1877 tmp1 = view(i, j);
1878 tmp2 = n_Div(tmp1, b, basecoeffs());
1879 rawset(i, j, tmp2);
1880 }
1881}
1882
1883// col(j, k) <- col(j,k)*matrix((a, c)(b, d))
1884// mostly used internally in the hnf and Howell stuff
1886{
1888 for (int i=1; i<=row; i++)
1889 {
1890 tmp1 = get(i, j);
1891 tmp2 = get(i, k);
1892 tmp3 = n_Mult(tmp1, a, basecoeffs());
1893 tmp4 = n_Mult(tmp2, b, basecoeffs());
1896
1897 n_InpMult(tmp1, c, basecoeffs());
1898 n_InpMult(tmp2, d, basecoeffs());
1901
1902 set(i, j, tmp3);
1903 set(i, k, tmp1);
1906 }
1907}
1908
1909
1910
1911//reduce all entries mod p. Does NOT change the coeffs type
1913{
1914 // produce the matrix in Z/pZ
1915 number tmp1, tmp2;
1916 for (int i=1; i<=row; i++)
1917 {
1918 for (int j=1; j<=col; j++)
1919 {
1920 tmp1 = get(i, j);
1921 tmp2 = n_IntMod(tmp1, p, basecoeffs());
1923 set(i, j, tmp2);
1924 }
1925 }
1926}
1927
1929{
1930 if (!nCoeffs_are_equal(a->basecoeffs(), b->basecoeffs()))
1931 {
1932 WerrorS("Error in bimMult. Coeffs do not agree!");
1933 return;
1934 }
1935 if ((a->rows() != c->rows()) || (b->cols() != c->cols()) || (a->cols() != b->rows()))
1936 {
1937 WerrorS("Error in bimMult. Dimensions do not agree!");
1938 return;
1939 }
1940 bigintmat *tmp = bimMult(a, b);
1941 c->copy(tmp);
1942
1943 delete tmp;
1944}
1945
1947{
1948 //write b = Ax + eps where eps is "small" in the sense of bounded by the
1949 //pivot entries in H. H does not need to be Howell (or HNF) but need
1950 //to be triagonal in the same direction.
1951 //b can have multiple columns.
1952#if 0
1953 PrintS("reduce_mod_howell: A:\n");
1954 A->Print();
1955 PrintS("\nb:\n");
1956 b->Print();
1957#endif
1958
1959 coeffs R = A->basecoeffs();
1960 assume(x->basecoeffs() == R);
1961 assume(b->basecoeffs() == R);
1962 assume(eps->basecoeffs() == R);
1963 if (!A->cols())
1964 {
1965 x->zero();
1966 eps->copy(b);
1967
1968#if 0
1969 PrintS("\nx:\n");
1970 x->Print();
1971 PrintS("\neps:\n");
1972 eps->Print();
1973 PrintS("\n****************************************\n");
1974#endif
1975 return;
1976 }
1977
1978 bigintmat * B = new bigintmat(b->rows(), 1, R);
1979 for(int i=1; i<= b->cols(); i++)
1980 {
1981 int A_col = A->cols();
1982 b->getcol(i, B);
1983 for(int j = B->rows(); j>0; j--)
1984 {
1985 number Ai = A->view(A->rows() - B->rows() + j, A_col);
1986 if (n_IsZero(Ai, R) &&
1987 n_IsZero(B->view(j, 1), R))
1988 {
1989 continue; //all is fine: 0*x = 0
1990 }
1991 else if (n_IsZero(B->view(j, 1), R))
1992 {
1993 x->rawset(x->rows() - B->rows() + j, i, n_Init(0, R));
1994 A_col--;
1995 }
1996 else if (n_IsZero(Ai, R))
1997 {
1998 A_col--;
1999 }
2000 else
2001 {
2002 // "solve" ax=b, possibly enlarging d
2003 number Bj = B->view(j, 1);
2004 number q = n_Div(Bj, Ai, R);
2005 x->rawset(x->rows() - B->rows() + j, i, q);
2006 for(int k=j; k>B->rows() - A->rows(); k--)
2007 {
2008 //B[k] = B[k] - x[k]A[k][j]
2009 number s = n_Mult(q, A->view(A->rows() - B->rows() + k, A_col), R);
2010 B->rawset(k, 1, n_Sub(B->view(k, 1), s, R));
2011 n_Delete(&s, R);
2012 }
2013 A_col--;
2014 }
2015 if (!A_col)
2016 {
2017 break;
2018 }
2019 }
2020 eps->setcol(i, B);
2021 }
2022 delete B;
2023#if 0
2024 PrintS("\nx:\n");
2025 x->Print();
2026 PrintS("\neps:\n");
2027 eps->Print();
2028 PrintS("\n****************************************\n");
2029#endif
2030}
2031
2033{
2034 coeffs R = A->basecoeffs();
2035 bigintmat *m = new bigintmat(A->rows()+A->cols(), A->cols(), R);
2036 m->copySubmatInto(A, 1, 1, A->rows(), A->cols(), A->cols()+1, 1);
2037 number one = n_Init(1, R);
2038 for(int i=1; i<= A->cols(); i++)
2039 m->set(i,i,one);
2040 n_Delete(&one, R);
2041 return m;
2042}
2043
2045{
2046 coeffs Z = A->basecoeffs(),
2047 Q = nInitChar(n_Q, 0);
2048 number den = n_Init(1, Z);
2049 nMapFunc f = n_SetMap(Q, Z);
2050
2051 for(int i=1; i<= A->rows(); i++)
2052 {
2053 for(int j=1; j<= A->cols(); j++)
2054 {
2055 number ad = n_Mult(den, A->view(i, j), Z);
2056 number re = n_IntMod(ad, N, Z);
2057 n_Delete(&ad, Z);
2058 number q = n_Farey(re, N, Z);
2059 n_Delete(&re, Z);
2060 if (!q)
2061 {
2062 n_Delete(&ad, Z);
2063 n_Delete(&den, Z);
2064 return NULL;
2065 }
2066
2067 number d = n_GetDenom(q, Q),
2068 n = n_GetNumerator(q, Q);
2069
2070 n_Delete(&q, Q);
2071 n_Delete(&ad, Z);
2072 number dz = f(d, Q, Z),
2073 nz = f(n, Q, Z);
2074 n_Delete(&d, Q);
2075 n_Delete(&n, Q);
2076
2077 if (!n_IsOne(dz, Z))
2078 {
2079 L->skalmult(dz, Z);
2080 n_InpMult(den, dz, Z);
2081#if 0
2082 PrintS("den increasing to ");
2083 n_Print(den, Z);
2084 PrintLn();
2085#endif
2086 }
2087 n_Delete(&dz, Z);
2088 L->rawset(i, j, nz);
2089 }
2090 }
2091
2092 nKillChar(Q);
2093 PrintS("bimFarey worked\n");
2094#if 0
2095 L->Print();
2096 PrintS("\n * 1/");
2097 n_Print(den, Z);
2098 PrintLn();
2099#endif
2100 return den;
2101}
2102
2103#ifdef HAVE_RINGS
2105 coeffs R = A->basecoeffs();
2106
2107 assume(getCoeffType(R) == n_Z);
2108
2109 number p = n_Init(536870909, R); // PreviousPrime(2^29); not clever
2110 coeffs Rp = numbercoeffs(p, R); // R/pR
2112 *m = prependIdentity(Ap),
2113 *Tp, *Hp;
2114 delete Ap;
2115
2116 m->howell();
2117 Hp = new bigintmat(A->rows(), A->cols(), Rp);
2118 Hp->copySubmatInto(m, A->cols()+1, 1, A->rows(), A->cols(), 1, 1);
2119 Tp = new bigintmat(A->cols(), A->cols(), Rp);
2120 Tp->copySubmatInto(m, 1, 1, A->cols(), A->cols(), 1, 1);
2121
2122 int i, j;
2123
2124 for(i=1; i<= A->cols(); i++)
2125 {
2126 for(j=m->rows(); j>A->cols(); j--)
2127 {
2128 if (!n_IsZero(m->view(j, i), Rp)) break;
2129 }
2130 if (j>A->cols()) break;
2131 }
2132// Print("Found nullity (kern dim) of %d\n", i-1);
2133 bigintmat * kp = new bigintmat(A->cols(), i-1, Rp);
2134 kp->copySubmatInto(Tp, 1, 1, A->cols(), i-1, 1, 1);
2135 kp->howell();
2136
2137 delete m;
2138
2139 //Hp is the mod-p howell form
2140 //Tp the transformation, mod p
2141 //kp a basis for the kernel, in howell form, mod p
2142
2143 bigintmat * eps_p = new bigintmat(B->rows(), B->cols(), Rp),
2144 * x_p = new bigintmat(A->cols(), B->cols(), Rp),
2145 * fps_p = new bigintmat(kp->cols(), B->cols(), Rp);
2146
2147 //initial solution
2148
2149 number zero = n_Init(0, R);
2150 x->skalmult(zero, R);
2151 n_Delete(&zero, R);
2152
2153 bigintmat * b = new bigintmat(B);
2154 number pp = n_Init(1, R);
2155 i = 1;
2156 do
2157 {
2158 bigintmat * b_p = bimChangeCoeff(b, Rp), * s;
2159 bigintmat * t1, *t2;
2161 delete b_p;
2162 if (!eps_p->isZero())
2163 {
2164 PrintS("no solution, since no modular solution\n");
2165
2166 delete eps_p;
2167 delete x_p;
2168 delete Hp;
2169 delete kp;
2170 delete Tp;
2171 delete b;
2172 n_Delete(&pp, R);
2173 n_Delete(&p, R);
2174 nKillChar(Rp);
2175
2176 return NULL;
2177 }
2178 t1 = bimMult(Tp, x_p);
2179 delete x_p;
2180 x_p = t1;
2181 reduce_mod_howell(kp, x_p, x_p, fps_p); //we're not all interested in fps_p
2182 s = bimChangeCoeff(x_p, R);
2183 t1 = bimMult(A, s);
2184 t2 = bimSub(b, t1);
2185 t2->skaldiv(p);
2186 delete b;
2187 delete t1;
2188 b = t2;
2189 s->skalmult(pp, R);
2190 t1 = bimAdd(x, s);
2191 delete s;
2192 x->swapMatrix(t1);
2193 delete t1;
2194
2195 if(kern && i==1)
2196 {
2198 t1 = bimMult(A, ker);
2199 t1->skaldiv(p);
2200 t1->skalmult(n_Init(-1, R), R);
2201 b->appendCol(t1);
2202 delete t1;
2203 x->appendCol(ker);
2204 delete ker;
2205 x_p->extendCols(kp->cols());
2206 eps_p->extendCols(kp->cols());
2207 fps_p->extendCols(kp->cols());
2208 }
2209
2210 n_InpMult(pp, p, R);
2211
2212 if (b->isZero())
2213 {
2214 //exact solution found, stop
2215 delete eps_p;
2216 delete fps_p;
2217 delete x_p;
2218 delete Hp;
2219 delete kp;
2220 delete Tp;
2221 delete b;
2222 n_Delete(&pp, R);
2223 n_Delete(&p, R);
2224 nKillChar(Rp);
2225
2226 return n_Init(1, R);
2227 }
2228 else
2229 {
2230 bigintmat *y = new bigintmat(x->rows(), x->cols(), R);
2231 number d = bimFarey(x, pp, y);
2232 if (d)
2233 {
2234 bigintmat *c = bimMult(A, y);
2235 bigintmat *bd = new bigintmat(B);
2236 bd->skalmult(d, R);
2237 if (kern)
2238 {
2239 bd->extendCols(kp->cols());
2240 }
2241 if (*c == *bd)
2242 {
2243 x->swapMatrix(y);
2244 delete y;
2245 delete c;
2246 if (kern)
2247 {
2248 y = new bigintmat(x->rows(), B->cols(), R);
2249 c = new bigintmat(x->rows(), kp->cols(), R);
2250 x->splitcol(y, c);
2251 x->swapMatrix(y);
2252 delete y;
2253 kern->swapMatrix(c);
2254 delete c;
2255 }
2256
2257 delete bd;
2258
2259 delete eps_p;
2260 delete fps_p;
2261 delete x_p;
2262 delete Hp;
2263 delete kp;
2264 delete Tp;
2265 delete b;
2266 n_Delete(&pp, R);
2267 n_Delete(&p, R);
2268 nKillChar(Rp);
2269
2270 return d;
2271 }
2272 delete c;
2273 delete bd;
2274 n_Delete(&d, R);
2275 }
2276 delete y;
2277 }
2278 i++;
2279 } while (1);
2280 delete eps_p;
2281 delete fps_p;
2282 delete x_p;
2283 delete Hp;
2284 delete kp;
2285 delete Tp;
2286 n_Delete(&pp, R);
2287 n_Delete(&p, R);
2288 nKillChar(Rp);
2289 return NULL;
2290}
2291#endif
2292
2293//TODO: re-write using reduce_mod_howell
2295{
2296 // try to solve Ax=b, more precisely, find
2297 // number d
2298 // bigintmat x
2299 // sth. Ax=db
2300 // where d is small-ish (divides the determinant of A if this makes sense)
2301 // return 0 if there is no solution.
2302 //
2303 // if kern is non-NULL, return a basis for the kernel
2304
2305 //Algo: we do row-howell (triangular matrix). The idea is
2306 // Ax = b <=> AT T^-1x = b
2307 // y := T^-1 x, solve AT y = b
2308 // and return Ty.
2309 //Howell does not compute the trafo, hence we need to cheat:
2310 //B := (I_n | A^t)^t, then the top part of the Howell form of
2311 //B will give a useful trafo
2312 //Then we can find x by back-substitution and lcm/gcd to find the denominator
2313 //The defining property of Howell makes this work.
2314
2315 coeffs R = A->basecoeffs();
2317 m->howell(); // since m contains the identity, we'll have A->cols()
2318 // many cols.
2319 number den = n_Init(1, R);
2320
2321 bigintmat * B = new bigintmat(A->rows(), 1, R);
2322 for(int i=1; i<= b->cols(); i++)
2323 {
2324 int A_col = A->cols();
2325 b->getcol(i, B);
2326 B->skalmult(den, R);
2327 for(int j = B->rows(); j>0; j--)
2328 {
2329 number Ai = m->view(m->rows()-B->rows() + j, A_col);
2330 if (n_IsZero(Ai, R) &&
2331 n_IsZero(B->view(j, 1), R))
2332 {
2333 continue; //all is fine: 0*x = 0
2334 }
2335 else if (n_IsZero(B->view(j, 1), R))
2336 {
2337 x->rawset(x->rows() - B->rows() + j, i, n_Init(0, R));
2338 A_col--;
2339 }
2340 else if (n_IsZero(Ai, R))
2341 {
2342 delete m;
2343 delete B;
2344 n_Delete(&den, R);
2345 return 0;
2346 }
2347 else
2348 {
2349 // solve ax=db, possibly enlarging d
2350 // so x = db/a
2351 number Bj = B->view(j, 1);
2352 number g = n_Gcd(Bj, Ai, R);
2353 number xi;
2354 if (n_Equal(Ai, g, R))
2355 { //good: den stable!
2356 xi = n_Div(Bj, Ai, R);
2357 }
2358 else
2359 { //den <- den * (a/g), so old sol. needs to be adjusted
2360 number inc_d = n_Div(Ai, g, R);
2361 n_InpMult(den, inc_d, R);
2362 x->skalmult(inc_d, R);
2363 B->skalmult(inc_d, R);
2364 xi = n_Div(Bj, g, R);
2365 n_Delete(&inc_d, R);
2366 } //now for the back-substitution:
2367 x->rawset(x->rows() - B->rows() + j, i, xi);
2368 for(int k=j; k>0; k--)
2369 {
2370 //B[k] = B[k] - x[k]A[k][j]
2371 number s = n_Mult(xi, m->view(m->rows()-B->rows() + k, A_col), R);
2372 B->rawset(k, 1, n_Sub(B->view(k, 1), s, R));
2373 n_Delete(&s, R);
2374 }
2375 n_Delete(&g, R);
2376 A_col--;
2377 }
2378 if (!A_col)
2379 {
2380 if (B->isZero()) break;
2381 else
2382 {
2383 delete m;
2384 delete B;
2385 n_Delete(&den, R);
2386 return 0;
2387 }
2388 }
2389 }
2390 }
2391 delete B;
2392 bigintmat *T = new bigintmat(A->cols(), A->cols(), R);
2393 T->copySubmatInto(m, 1, 1, A->cols(), A->cols(), 1, 1);
2394 if (kern)
2395 {
2396 int i, j;
2397 for(i=1; i<= A->cols(); i++)
2398 {
2399 for(j=m->rows(); j>A->cols(); j--)
2400 {
2401 if (!n_IsZero(m->view(j, i), R)) break;
2402 }
2403 if (j>A->cols()) break;
2404 }
2405 Print("Found nullity (kern dim) of %d\n", i-1);
2406 bigintmat * ker = new bigintmat(A->rows(), i-1, R);
2407 ker->copySubmatInto(T, 1, 1, A->rows(), i-1, 1, 1);
2408 kern->swapMatrix(ker);
2409 delete ker;
2410 }
2411 delete m;
2412 bigintmat * y = bimMult(T, x);
2413 x->swapMatrix(y);
2414 delete y;
2415 x->simplifyContentDen(&den);
2416#if 0
2417 PrintS("sol = 1/");
2418 n_Print(den, R);
2419 PrintS(" *\n");
2420 x->Print();
2421 PrintLn();
2422#endif
2423 return den;
2424}
2425
2427{
2428#if 0
2429 PrintS("Solve Ax=b for A=\n");
2430 A->Print();
2431 PrintS("\nb = \n");
2432 b->Print();
2433 PrintS("\nx = \n");
2434 x->Print();
2435 PrintLn();
2436#endif
2437
2438 coeffs R = A->basecoeffs();
2439 assume (R == b->basecoeffs());
2440 assume (R == x->basecoeffs());
2441 assume ((x->cols() == b->cols()) && (x->rows() == A->cols()) && (A->rows() == b->rows()));
2442
2443 switch (getCoeffType(R))
2444 {
2445 #ifdef HAVE_RINGS
2446 case n_Z:
2447 return solveAx_dixon(A, b, x, NULL);
2448 case n_Zn:
2449 case n_Znm:
2450 case n_Z2m:
2451 return solveAx_howell(A, b, x, NULL);
2452 #endif
2453 case n_Zp:
2454 case n_Q:
2455 case n_GF:
2456 case n_algExt:
2457 case n_transExt:
2458 WarnS("have field, should use Gauss or better");
2459 break;
2460 default:
2461 if (R->cfXExtGcd && R->cfAnn)
2462 { //assume it's Euclidean
2463 return solveAx_howell(A, b, x, NULL);
2464 }
2465 WerrorS("have no solve algorithm");
2466 break;
2467 }
2468 return NULL;
2469}
2470
2472{
2473 bigintmat * t, *s, *a=A;
2474 coeffs R = a->basecoeffs();
2475 if (T)
2476 {
2477 *T = new bigintmat(a->cols(), a->cols(), R),
2478 (*T)->one();
2479 t = new bigintmat(*T);
2480 }
2481 else
2482 {
2483 t = *T;
2484 }
2485
2486 if (S)
2487 {
2488 *S = new bigintmat(a->rows(), a->rows(), R);
2489 (*S)->one();
2490 s = new bigintmat(*S);
2491 }
2492 else
2493 {
2494 s = *S;
2495 }
2496
2497 int flip=0;
2498 do
2499 {
2500 bigintmat * x, *X;
2501 if (flip)
2502 {
2503 x = s;
2504 X = *S;
2505 }
2506 else
2507 {
2508 x = t;
2509 X = *T;
2510 }
2511
2512 if (x)
2513 {
2514 x->one();
2515 bigintmat * r = new bigintmat(a->rows()+a->cols(), a->cols(), R);
2516 bigintmat * rw = new bigintmat(1, a->cols(), R);
2517 for(int i=0; i<a->cols(); i++)
2518 {
2519 x->getrow(i+1, rw);
2520 r->setrow(i+1, rw);
2521 }
2522 for (int i=0; i<a->rows(); i++)
2523 {
2524 a->getrow(i+1, rw);
2525 r->setrow(i+a->cols()+1, rw);
2526 }
2527 r->hnf();
2528 for(int i=0; i<a->cols(); i++)
2529 {
2530 r->getrow(i+1, rw);
2531 x->setrow(i+1, rw);
2532 }
2533 for(int i=0; i<a->rows(); i++)
2534 {
2535 r->getrow(i+a->cols()+1, rw);
2536 a->setrow(i+1, rw);
2537 }
2538 delete rw;
2539 delete r;
2540
2541#if 0
2542 Print("X: %ld\n", X);
2543 X->Print();
2544 Print("\nx: %ld\n", x);
2545 x->Print();
2546#endif
2547 bimMult(X, x, X);
2548#if 0
2549 Print("\n2:X: %ld %ld %ld\n", X, *S, *T);
2550 X->Print();
2551 Print("\n2:x: %ld\n", x);
2552 x->Print();
2553 PrintLn();
2554#endif
2555 }
2556 else
2557 {
2558 a->hnf();
2559 }
2560
2561 int diag = 1;
2562 for(int i=a->rows(); diag && i>0; i--)
2563 {
2564 for(int j=a->cols(); j>0; j--)
2565 {
2566 if ((a->rows()-i)!=(a->cols()-j) && !n_IsZero(a->view(i, j), R))
2567 {
2568 diag = 0;
2569 break;
2570 }
2571 }
2572 }
2573#if 0
2574 PrintS("Diag ? %d\n", diag);
2575 a->Print();
2576 PrintLn();
2577#endif
2578 if (diag) break;
2579
2580 a = a->transpose(); // leaks - I need to write inpTranspose
2581 flip = 1-flip;
2582 } while (1);
2583 if (flip)
2584 a = a->transpose();
2585
2586 if (S) *S = (*S)->transpose();
2587 if (s) delete s;
2588 if (t) delete t;
2589 A->copy(a);
2590}
2591
2592#ifdef HAVE_RINGS
2593//a "q-base" for the kernel of a.
2594//Should be re-done to use Howell rather than smith.
2595//can be done using solveAx now.
2597{
2598#if 0
2599 PrintS("Kernel of ");
2600 a->Print();
2601 PrintS(" modulo ");
2602 n_Print(p, q);
2603 PrintLn();
2604#endif
2605
2606 coeffs coe = numbercoeffs(p, q);
2607 bigintmat *m = bimChangeCoeff(a, coe), *U, *V;
2608 diagonalForm(m, &U, &V);
2609#if 0
2610 PrintS("\ndiag form: ");
2611 m->Print();
2612 PrintS("\nU:\n");
2613 U->Print();
2614 PrintS("\nV:\n");
2615 V->Print();
2616 PrintLn();
2617#endif
2618
2619 int rg = 0;
2620#undef MIN
2621#define MIN(a,b) (a < b ? a : b)
2622 for(rg=0; rg<MIN(m->rows(), m->cols()) && !n_IsZero(m->view(m->rows()-rg,m->cols()-rg), coe); rg++);
2623
2624 bigintmat * k = new bigintmat(m->cols(), m->rows(), coe);
2625 for(int i=0; i<rg; i++)
2626 {
2627 number A = n_Ann(m->view(m->rows()-i, m->cols()-i), coe);
2628 k->set(m->cols()-i, i+1, A);
2629 n_Delete(&A, coe);
2630 }
2631 for(int i=rg; i<m->cols(); i++)
2632 {
2633 k->set(m->cols()-i, i+1-rg, n_Init(1, coe));
2634 }
2635 bimMult(V, k, k);
2636 c->copy(bimChangeCoeff(k, q));
2637 return c->cols();
2638}
2639#endif
2640
2642{
2643 if ((r == NULL) || (s == NULL))
2644 return false;
2645 if (r == s)
2646 return true;
2647 if ((getCoeffType(r)==n_Z) && (getCoeffType(s)==n_Z))
2648 return true;
2649 if ((getCoeffType(r)==n_Zp) && (getCoeffType(s)==n_Zp))
2650 {
2651 if (r->ch == s->ch)
2652 return true;
2653 else
2654 return false;
2655 }
2656 // n_Zn stimmt wahrscheinlich noch nicht
2657 if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
2658 {
2659 if (r->ch == s->ch)
2660 return true;
2661 else
2662 return false;
2663 }
2664 if ((getCoeffType(r)==n_Q) && (getCoeffType(s)==n_Q))
2665 return true;
2666 // FALL n_Zn FEHLT NOCH!
2667 //if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
2668 return false;
2669}
2670
2672{
2673 coeffs r = basecoeffs();
2674 number g = get(1,1), h;
2675 int n=rows()*cols();
2676 for(int i=1; i<n && !n_IsOne(g, r); i++)
2677 {
2678 h = n_Gcd(g, view(i), r);
2679 n_Delete(&g, r);
2680 g=h;
2681 }
2682 return g;
2683}
2685{
2686 coeffs r = basecoeffs();
2687 number g = n_Copy(*d, r), h;
2688 int n=rows()*cols();
2689 for(int i=0; i<n && !n_IsOne(g, r); i++)
2690 {
2691 h = n_Gcd(g, view(i), r);
2692 n_Delete(&g, r);
2693 g=h;
2694 }
2695 *d = n_Div(*d, g, r);
2696 if (!n_IsOne(g, r))
2697 skaldiv(g);
2698}
All the auxiliary stuff.
#define TRUE
Definition auxiliary.h:100
#define FALSE
Definition auxiliary.h:96
static int si_min(const int a, const int b)
Definition auxiliary.h:125
static void reduce_mod_howell(bigintmat *A, bigintmat *b, bigintmat *eps, bigintmat *x)
int kernbase(bigintmat *a, bigintmat *c, number p, coeffs q)
a basis for the nullspace of a mod p: only used internally in Round2. Don't use it.
#define swap(_i, _j)
number solveAx(bigintmat *A, bigintmat *b, bigintmat *x)
solve Ax=b*d. x needs to be pre-allocated to the same number of columns as b. the minimal denominator...
static number solveAx_dixon(bigintmat *A, bigintmat *B, bigintmat *x, bigintmat *kern)
bigintmat * bimChangeCoeff(bigintmat *a, coeffs cnew)
Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.
#define MIN(a, b)
bigintmat * bimCopy(const bigintmat *b)
same as copy constructor - apart from it being able to accept NULL as input
Definition bigintmat.cc:405
static bigintmat * prependIdentity(bigintmat *A)
bool nCoeffs_are_equal(coeffs r, coeffs s)
intvec * bim2iv(bigintmat *b)
Definition bigintmat.cc:341
bigintmat * bimMult(bigintmat *a, bigintmat *b)
Definition bigintmat.cc:255
static int intArrSum(int *a, int length)
Definition bigintmat.cc:525
bigintmat * bimSub(bigintmat *a, bigintmat *b)
Definition bigintmat.cc:218
static number solveAx_howell(bigintmat *A, bigintmat *b, bigintmat *x, bigintmat *kern)
void diagonalForm(bigintmat *A, bigintmat **S, bigintmat **T)
bigintmat * iv2bim(intvec *b, const coeffs C)
Definition bigintmat.cc:349
static int findLongest(int *a, int length)
Definition bigintmat.cc:533
static int getShorter(int *a, int l, int j, int cols, int rows)
Definition bigintmat.cc:548
static number bimFarey(bigintmat *A, number N, bigintmat *L)
bool operator!=(const bigintmat &lhr, const bigintmat &rhr)
Definition bigintmat.cc:176
static coeffs numbercoeffs(number n, coeffs c)
create Z/nA of type n_Zn
Definition bigintmat.cc:21
bigintmat * bimAdd(bigintmat *a, bigintmat *b)
Matrix-Add/-Sub/-Mult so oder mit operator+/-/* ? @Note: NULL as a result means an error (non-compati...
Definition bigintmat.cc:182
bool operator==(const bigintmat &lhr, const bigintmat &rhr)
Definition bigintmat.cc:159
#define BIMATELEM(M, I, J)
Definition bigintmat.h:133
bigintmat * bimChangeCoeff(bigintmat *a, coeffs cnew)
Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.
bool nCoeffs_are_equal(coeffs r, coeffs s)
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition cf_gcd.cc:676
CanonicalForm den(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition cfEzgcd.cc:56
int l
Definition cfEzgcd.cc:100
int m
Definition cfEzgcd.cc:128
int i
Definition cfEzgcd.cc:132
int k
Definition cfEzgcd.cc:99
Variable x
Definition cfModGcd.cc:4090
int p
Definition cfModGcd.cc:4086
g
Definition cfModGcd.cc:4098
CanonicalForm cf
Definition cfModGcd.cc:4091
CanonicalForm b
Definition cfModGcd.cc:4111
FILE * f
Definition checklibs.c:9
CF_NO_INLINE bool isZero() const
Matrices of numbers.
Definition bigintmat.h:51
void Print()
IO: simply prints the matrix to the current output (screen?)
Definition bigintmat.cc:439
void colskaldiv(int j, number b)
Macht Ganzzahldivision aller j-ten Spalteneinträge mit b.
number det()
det (via LaPlace in general, hnf for euc. rings)
void splitrow(bigintmat *a, bigintmat *b)
Speichert in Matrix a den oberen, in b den unteren Teil der Matrix, vorausgesetzt die Dimensionen sti...
int isOne()
is matrix is identity
void zero()
Setzt alle Einträge auf 0.
void appendCol(bigintmat *a)
horizontally join the matrices, m <- m|a
number trace()
the trace ....
bool addrow(int i, int j, number a, coeffs c)
... Zeile ...
Definition bigintmat.cc:979
void swap(int i, int j)
swap columns i and j
Definition bigintmat.cc:681
void coltransform(int i, int j, number a, number b, number c, number d)
transforms cols (i,j) using the 2x2 matrix ((a,b)(c,d)) (hopefully)
bool addcol(int i, int j, number a, coeffs c)
addiert a-faches der j-ten Spalte zur i-ten dazu
Definition bigintmat.cc:955
number * v
Definition bigintmat.h:54
void hnf()
transforms INPLACE to HNF
char * StringAsPrinted()
Returns a string as it would have been printed in the interpreter.
Definition bigintmat.cc:447
void swapMatrix(bigintmat *a)
int cols() const
Definition bigintmat.h:144
int isZero()
number hnfdet()
det via HNF Primzahlen als long long int, müssen noch in number umgewandelt werden?
int findnonzero(int i)
find index of 1st non-zero entry in row i
Definition bigintmat.cc:719
bigintmat * modhnf(number p, coeffs c)
computes HNF(this | p*I)
void setcol(int j, bigintmat *m)
Setzt j-te Spalte gleich übergebenem Vektor (Matrix) m.
Definition bigintmat.cc:822
bool add(bigintmat *b)
Addiert zur Matrix die Matrix b dazu. Return false => an error occurred.
Definition bigintmat.cc:890
int * getwid(int maxwid)
Definition bigintmat.cc:576
void inpMult(number bintop, const coeffs C=NULL)
inplace version of skalar mult. CHANGES input.
Definition bigintmat.cc:145
bigintmat * transpose()
Definition bigintmat.cc:37
void setrow(int i, bigintmat *m)
Setzt i-te Zeile gleich übergebenem Vektor (Matrix) m.
Definition bigintmat.cc:856
void Write()
IO: writes the matrix into the current internal string buffer which must be started/ allocated before...
Definition bigintmat.cc:413
void rowskalmult(int i, number a, coeffs c)
... Zeile ...
bool skalmult(number b, coeffs c)
Multipliziert zur Matrix den Skalar b hinzu.
Definition bigintmat.cc:934
void simplifyContentDen(number *den)
ensures that Gcd(den, content)=1 enden hier wieder
void extendCols(int i)
append i zero-columns to the matrix
void splitcol(bigintmat *a, bigintmat *b)
... linken ... rechten ...
number content()
the content, the gcd of all entries. Only makes sense for Euclidean rings (or possibly constructive P...
void skaldiv(number b)
Macht Ganzzahldivision aller Matrixeinträge mit b.
void colskalmult(int i, number a, coeffs c)
Multipliziert zur i-ten Spalte den Skalar a hinzu.
int colIsZero(int i)
bool copy(bigintmat *b)
Kopiert Einträge von b auf Bigintmat.
int index(int r, int c) const
helper function to map from 2-dim coordinates, starting by 1 to 1-dim coordinate, starting by 0
Definition bigintmat.h:161
void getcol(int j, bigintmat *a)
copies the j-th column into the matrix a - which needs to be pre-allocated with the correct size.
Definition bigintmat.cc:743
number pseudoinv(bigintmat *a)
Speichert in Matrix a die Pseudoinverse, liefert den Nenner zurück.
int rows() const
Definition bigintmat.h:145
number get(int i, int j) const
get a copy of an entry. NOTE: starts at [1,1]
Definition bigintmat.cc:119
void pprint(int maxwid)
Definition bigintmat.cc:606
void getColRange(int j, int no, bigintmat *a)
copies the no-columns staring by j (so j...j+no-1) into the pre-allocated a
Definition bigintmat.cc:774
void concatcol(bigintmat *a, bigintmat *b)
int findcolnonzero(int j)
find index of 1st non-zero entry in column j
Definition bigintmat.cc:731
void copySubmatInto(bigintmat *, int sr, int sc, int nr, int nc, int tr, int tc)
copy the submatrix of b, staring at (a,b) having n rows, m cols into the given matrix at pos....
number view(int i, int j) const
view an entry an entry. NOTE: starts at [1,1]
Definition bigintmat.cc:127
void inpTranspose()
transpose in place
Definition bigintmat.cc:50
void one()
Macht Matrix (Falls quadratisch) zu Einheitsmatrix.
void howell()
dito, but Howell form (only different for zero-divsors)
void rawset(int i, number n, const coeffs C=NULL)
replace an entry with the given number n (only delete old). NOTE: starts at [0]. Should be named set_...
Definition bigintmat.h:196
void operator*=(int intop)
UEberladener *=-Operator (fuer int und bigint) Frage hier: *= verwenden oder lieber = und * einzeln?...
Definition bigintmat.cc:136
coeffs basecoeffs() const
Definition bigintmat.h:146
void getrow(int i, bigintmat *a)
Schreibt i-te Zeile in Vektor (Matrix) a.
Definition bigintmat.cc:787
void concatrow(bigintmat *a, bigintmat *b)
Fügt zwei Matrixen untereinander/nebeneinander in gegebene Matrix ein, bzw spaltet gegebenen Matrix a...
void set(int i, int j, number n, const coeffs C=NULL)
replace an entry with a copy (delete old + copy new!). NOTE: starts at [1,1]
Definition bigintmat.cc:95
bigintmat * elim(int i, int j)
Liefert Streichungsmatrix (i-te Zeile und j-te Spalte gestrichen) zurück.
int compare(const bigintmat *op) const
Definition bigintmat.cc:362
bool sub(bigintmat *b)
Subtrahiert ...
Definition bigintmat.cc:912
char * String()
IO: String returns a singular string containing the matrix, needs freeing afterwards.
Definition bigintmat.cc:432
void swaprow(int i, int j)
swap rows i and j
Definition bigintmat.cc:700
void mod(number p)
Reduziert komplette Matrix modulo p.
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of 'a' and 'b', i.e., a*b
Definition coeffs.h:640
static FORCE_INLINE void number2mpz(number n, coeffs c, mpz_t m)
Definition coeffs.h:991
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ....
Definition coeffs.h:551
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition coeffs.h:455
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of 'a' and 'b', i.e., a+b
Definition coeffs.h:654
static FORCE_INLINE number n_GetDenom(number &n, const coeffs r)
return the denominator of n (if elements of r are by nature not fractional, result is 1)
Definition coeffs.h:607
@ n_GF
\GF{p^n < 2^16}
Definition coeffs.h:32
@ n_Q
rational (GMP) numbers
Definition coeffs.h:30
@ n_Znm
only used if HAVE_RINGS is defined
Definition coeffs.h:45
@ n_algExt
used for all algebraic extensions, i.e., the top-most extension in an extension tower is algebraic
Definition coeffs.h:35
@ n_Zn
only used if HAVE_RINGS is defined
Definition coeffs.h:44
@ n_Z2m
only used if HAVE_RINGS is defined
Definition coeffs.h:46
@ n_Zp
\F{p < 2^31}
Definition coeffs.h:29
@ n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition coeffs.h:38
@ n_Z
only used if HAVE_RINGS is defined
Definition coeffs.h:43
static FORCE_INLINE number n_Gcd(number a, number b, const coeffs r)
in Z: return the gcd of 'a' and 'b' in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ,...
Definition coeffs.h:668
static FORCE_INLINE number n_QuotRem(number a, number b, number *q, const coeffs r)
Definition coeffs.h:685
static FORCE_INLINE char * nCoeffString(const coeffs cf)
TODO: make it a virtual method of coeffs, together with: Decompose & Compose, rParameter & rPar.
Definition coeffs.h:963
static FORCE_INLINE BOOLEAN n_GreaterZero(number n, const coeffs r)
ordered fields: TRUE iff 'n' is positive; in Z/pZ: TRUE iff 0 < m <= roundedBelow(p/2),...
Definition coeffs.h:498
static FORCE_INLINE number n_Ann(number a, const coeffs r)
if r is a ring with zero divisors, return an annihilator!=0 of b otherwise return NULL
Definition coeffs.h:683
void n_Print(number &a, const coeffs r)
print a number (BEWARE of string buffers!) mostly for debugging
Definition numbers.cc:673
static FORCE_INLINE nMapFunc n_SetMap(const coeffs src, const coeffs dst)
set the mapping function pointers for translating numbers from src to dst
Definition coeffs.h:704
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
Definition coeffs.h:561
static FORCE_INLINE number n_Farey(number a, number b, const coeffs r)
Definition coeffs.h:771
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition coeffs.h:619
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition numbers.cc:419
static FORCE_INLINE BOOLEAN n_Greater(number a, number b, const coeffs r)
ordered fields: TRUE iff 'a' is larger than 'b'; in Z/pZ: TRUE iff la > lb, where la and lb are the l...
Definition coeffs.h:515
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition coeffs.h:468
static FORCE_INLINE number n_GetUnit(number n, const coeffs r)
in Z: 1 in Z/kZ (where k is not a prime): largest divisor of n (taken in Z) that is co-prime with k i...
Definition coeffs.h:536
static FORCE_INLINE number n_Sub(number a, number b, const coeffs r)
return the difference of 'a' and 'b', i.e., a-b
Definition coeffs.h:659
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition coeffs.h:429
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition coeffs.h:459
static FORCE_INLINE void n_Write(number n, const coeffs r, const BOOLEAN bShortOut=TRUE)
Definition coeffs.h:595
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition coeffs.h:542
static FORCE_INLINE number n_IntMod(number a, number b, const coeffs r)
for r a field, return n_Init(0,r) always: n_Div(a,b,r)*b+n_IntMod(a,b,r)==a n_IntMod(a,...
Definition coeffs.h:632
static FORCE_INLINE void n_InpMult(number &a, number b, const coeffs r)
multiplication of 'a' and 'b'; replacement of 'a' by the product a*b
Definition coeffs.h:645
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff 'a' and 'b' represent the same number; they may have different representations.
Definition coeffs.h:464
static FORCE_INLINE number n_GetNumerator(number &n, const coeffs r)
return the numerator of n (if elements of r are by nature not fractional, result is n)
Definition coeffs.h:612
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition coeffs.h:80
void nKillChar(coeffs r)
undo all initialisations
Definition numbers.cc:574
static FORCE_INLINE void n_InpAdd(number &a, number b, const coeffs r)
addition of 'a' and 'b'; replacement of 'a' by the sum a+b
Definition coeffs.h:650
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition coeffs.h:472
static FORCE_INLINE number n_XExtGcd(number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
Definition coeffs.h:677
#define Print
Definition emacs.cc:80
#define WarnS
Definition emacs.cc:78
#define StringAppend
Definition emacs.cc:79
const CanonicalForm int s
Definition facAbsFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition facAbsFact.cc:53
b *CanonicalForm B
Definition facBivar.cc:52
CFList tmp1
Definition facFqBivar.cc:75
CFList tmp2
Definition facFqBivar.cc:75
int j
Definition facHensel.cc:110
fq_nmod_poly_t prod
Definition facHensel.cc:100
static int min(int a, int b)
Definition fast_mult.cc:268
void WerrorS(const char *s)
Definition feFopen.cc:24
std::pair< ideal, ring > flip(const ideal I, const ring r, const gfan::ZVector interiorPoint, const gfan::ZVector facetNormal, const gfan::ZVector adjustedInteriorPoint, const gfan::ZVector adjustedFacetNormal)
Definition flip.cc:17
static BOOLEAN length(leftv result, leftv arg)
Definition interval.cc:257
STATIC_VAR jList * T
Definition janet.cc:30
STATIC_VAR Poly * h
Definition janet.cc:971
int lcm(unsigned long *l, unsigned long *a, unsigned long *b, unsigned long p, int dega, int degb)
Definition minpoly.cc:709
#define assume(x)
Definition mod2.h:387
The main handler for Singular numbers which are suitable for Singular polynomials.
#define omAlloc(size)
#define omFree(addr)
#define omAlloc0(size)
#define NULL
Definition omList.c:12
static int index(p_Length length, p_Ord ord)
void StringSetS(const char *st)
Definition reporter.cc:128
void StringAppendS(const char *st)
Definition reporter.cc:107
void PrintS(const char *s)
Definition reporter.cc:284
char * StringEndS()
Definition reporter.cc:151
void PrintLn()
Definition reporter.cc:310
void Werror(const char *fmt,...)
Definition reporter.cc:189
#define R
Definition sirandom.c:27
#define A
Definition sirandom.c:24
#define Q
Definition sirandom.c:26
int gcd(int a, int b)