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
103number bigintmat::get(int i) const
104{
105 assume (i >= 0);
106 assume (i<rows()*cols());
107
108 return n_Copy(v[i], basecoeffs());
109}
110
111number bigintmat::view(int i) const
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{
138 number iop = n_Init(intop, basecoeffs());
139
140 inpMult(iop, basecoeffs());
141
142 n_Delete(&iop, basecoeffs());
143}
144
145void bigintmat::inpMult(number bintop, const coeffs C)
146{
147 assume (C == NULL || C == basecoeffs());
148
149 const int l = rows() * cols();
150
151 for (int i=0; i < l; i++)
152 n_InpMult(v[i], bintop, basecoeffs());
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
320bigintmat * bimMult(bigintmat * a, number b, const coeffs cf)
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 StringAppendS("[ ");
418 for(int i=1; i<= m; i++)
419 {
420 StringAppendS("[ ");
421 for(int j=1; j< n; j++)
422 {
423 n_Write(v[(i-1)*n+j-1], basecoeffs());
424 StringAppendS(", ");
425 }
426 if (n) n_Write(v[i*n-1], basecoeffs());
427 StringAppendS(" ]");
428 if (i<m)
429 {
430 StringAppendS(", ");
431 }
432 }
433 StringAppendS(" ] ");
434}
435
437{
438 StringSetS("");
439 Write();
440 return StringEndS();
441}
442
444{
445 char * s = String();
446 PrintS(s);
447 omFree(s);
448}
449
450
452{
453 if ((col==0) || (row==0))
454 return NULL;
455 else
456 {
457 int * colwid = getwid(80);
458 if (colwid == NULL)
459 {
460 WerrorS("not enough space to print bigintmat");
461 WerrorS("try string(...) for a unformatted output");
462 return NULL;
463 }
464 char * ps;
465 int slength = 0;
466 for (int j=0; j<col; j++)
467 slength += colwid[j]*row;
468 slength += col*row+row;
469 ps = (char*) omAlloc0(sizeof(char)*(slength));
470 int pos = 0;
471 for (int i=0; i<col*row; i++)
472 {
473 StringSetS("");
474 n_Write(v[i], basecoeffs());
475 char * ts = StringEndS();
476 const int _nl = strlen(ts);
477 int cj = i%col;
478 if (_nl > colwid[cj])
479 {
480 StringSetS("");
481 int ci = i/col;
482 StringAppend("[%d,%d]", ci+1, cj+1);
483 char * ph = StringEndS();
484 int phl = strlen(ph);
485 if (phl > colwid[cj])
486 {
487 for (int j=0; j<colwid[cj]-1; j++)
488 ps[pos+j] = ' ';
489 ps[pos+colwid[cj]-1] = '*';
490 }
491 else
492 {
493 for (int j=0; j<colwid[cj]-phl; j++)
494 ps[pos+j] = ' ';
495 for (int j=0; j<phl; j++)
496 ps[pos+colwid[cj]-phl+j] = ph[j];
497 }
498 omFree(ph);
499 }
500 else // Mit Leerzeichen auffüllen und zahl reinschreiben
501 {
502 for (int j=0; j<(colwid[cj]-_nl); j++)
503 ps[pos+j] = ' ';
504 for (int j=0; j<_nl; j++)
505 ps[pos+colwid[cj]-_nl+j] = ts[j];
506 }
507 // ", " und (evtl) "\n" einfügen
508 if ((i+1)%col == 0)
509 {
510 if (i != col*row-1)
511 {
512 ps[pos+colwid[cj]] = ',';
513 ps[pos+colwid[cj]+1] = '\n';
514 pos += colwid[cj]+2;
515 }
516 }
517 else
518 {
519 ps[pos+colwid[cj]] = ',';
520 pos += colwid[cj]+1;
521 }
522 omFree(ts); // Hier ts zerstören
523 }
524 return(ps);
525 // omFree(ps);
526}
527}
528
529static int intArrSum(int * a, int length)
530{
531 int sum = 0;
532 for (int i=0; i<length; i++)
533 sum += a[i];
534 return sum;
535}
536
537static int findLongest(int * a, int length)
538{
539 int l = 0;
540 int index;
541 for (int i=0; i<length; i++)
542 {
543 if (a[i] > l)
544 {
545 l = a[i];
546 index = i;
547 }
548 }
549 return index;
550}
551
552static int getShorter (int * a, int l, int j, int cols, int rows)
553{
554 int sndlong = 0;
555 int min;
556 for (int i=0; i<rows; i++)
557 {
558 int index = cols*i+j;
559 if ((a[index] > sndlong) && (a[index] < l))
560 {
561 min = floor(log10((double)cols))+floor(log10((double)rows))+5;
562 if ((a[index] < min) && (min < l))
563 sndlong = min;
564 else
565 sndlong = a[index];
566 }
567 }
568 if (sndlong == 0)
569 {
570 min = floor(log10((double)cols))+floor(log10((double)rows))+5;
571 if (min < l)
572 sndlong = min;
573 else
574 sndlong = 1;
575 }
576 return sndlong;
577}
578
579
580int * bigintmat::getwid(int maxwid)
581{
582 int const c = /*2**/(col-1)+1;
583 int * wv = (int*)omAlloc(sizeof(int)*col*row);
584 int * cwv = (int*)omAlloc(sizeof(int)*col);
585 for (int j=0; j<col; j++)
586 {
587 cwv[j] = 0;
588 for (int i=0; i<row; i++)
589 {
590 StringSetS("");
591 n_Write(v[col*i+j], basecoeffs());
592 char * tmp = StringEndS();
593 const int _nl = strlen(tmp);
594 wv[col*i+j] = _nl;
595 if (_nl > cwv[j]) cwv[j]=_nl;
596 omFree(tmp);
597 }
598 }
599
600 // Groesse verkleinern, bis < maxwid
601 if (intArrSum(cwv, col)+c > maxwid)
602 {
603 int j = findLongest(cwv, col);
604 cwv[j] = getShorter(wv, cwv[j], j, col, row);
605 }
606 omFree(wv);
607 return cwv;
608}
609
610void bigintmat::pprint(int maxwid)
611{
612 if ((col==0) || (row==0))
613 PrintS("");
614 else
615 {
616 int * colwid = getwid(maxwid);
617 char * ps;
618 int slength = 0;
619 for (int j=0; j<col; j++)
620 slength += colwid[j]*row;
621 slength += col*row+row;
622 ps = (char*) omAlloc0(sizeof(char)*(slength));
623 int pos = 0;
624 for (int i=0; i<col*row; i++)
625 {
626 StringSetS("");
627 n_Write(v[i], basecoeffs());
628 char * ts = StringEndS();
629 const int _nl = strlen(ts);
630 int cj = i%col;
631 if (_nl > colwid[cj])
632 {
633 StringSetS("");
634 int ci = i/col;
635 StringAppend("[%d,%d]", ci+1, cj+1);
636 char * ph = StringEndS();
637 int phl = strlen(ph);
638 if (phl > colwid[cj])
639 {
640 for (int j=0; j<colwid[cj]-1; j++)
641 ps[pos+j] = ' ';
642 ps[pos+colwid[cj]-1] = '*';
643 }
644 else
645 {
646 for (int j=0; j<colwid[cj]-phl; j++)
647 ps[pos+j] = ' ';
648 for (int j=0; j<phl; j++)
649 ps[pos+colwid[cj]-phl+j] = ph[j];
650 }
651 omFree(ph);
652 }
653 else // Mit Leerzeichen auffüllen und zahl reinschreiben
654 {
655 for (int j=0; j<colwid[cj]-_nl; j++)
656 ps[pos+j] = ' ';
657 for (int j=0; j<_nl; j++)
658 ps[pos+colwid[cj]-_nl+j] = ts[j];
659 }
660 // ", " und (evtl) "\n" einfügen
661 if ((i+1)%col == 0)
662 {
663 if (i != col*row-1)
664 {
665 ps[pos+colwid[cj]] = ',';
666 ps[pos+colwid[cj]+1] = '\n';
667 pos += colwid[cj]+2;
668 }
669 }
670 else
671 {
672 ps[pos+colwid[cj]] = ',';
673 pos += colwid[cj]+1;
674 }
675
676 omFree(ts); // Hier ts zerstören
677 }
678 PrintS(ps);
679 omFree(ps);
680 }
681}
682
683
684//swaps columns i and j
685void bigintmat::swap(int i, int j)
686{
687 if ((i <= col) && (j <= col) && (i>0) && (j>0))
688 {
689 number tmp;
690 number t;
691 for (int k=1; k<=row; k++)
692 {
693 tmp = get(k, i);
694 t = view(k, j);
695 set(k, i, t);
696 set(k, j, tmp);
697 n_Delete(&tmp, basecoeffs());
698 }
699 }
700 else
701 WerrorS("Error in swap");
702}
703
704void bigintmat::swaprow(int i, int j)
705{
706 if ((i <= row) && (j <= row) && (i>0) && (j>0))
707 {
708 number tmp;
709 number t;
710 for (int k=1; k<=col; k++)
711 {
712 tmp = get(i, k);
713 t = view(j, k);
714 set(i, k, t);
715 set(j, k, tmp);
716 n_Delete(&tmp, basecoeffs());
717 }
718 }
719 else
720 WerrorS("Error in swaprow");
721}
722
724{
725 for (int j=1; j<=col; j++)
726 {
727 if (!n_IsZero(view(i,j), basecoeffs()))
728 {
729 return j;
730 }
731 }
732 return 0;
733}
734
736{
737 for (int i=row; i>=1; i--)
738 {
739 if (!n_IsZero(view(i,j), basecoeffs()))
740 {
741 return i;
742 }
743 }
744 return 0;
745}
746
748{
749 assume((j<=col) && (j>=1));
750 if (((a->rows() != row) || (a->cols() != 1)) && ((a->rows() != 1) || (a->cols() != row)))
751 {
752 assume(0);
753 WerrorS("Error in getcol. Dimensions must agree!");
754 return;
755 }
757 {
759 number t1, t2;
760 for (int i=1; i<=row;i++)
761 {
762 t1 = get(i,j);
763 t2 = f(t1, basecoeffs(), a->basecoeffs());
764 a->set(i-1,t1);
765 n_Delete(&t1, basecoeffs());
766 n_Delete(&t2, a->basecoeffs());
767 }
768 return;
769 }
770 number t1;
771 for (int i=1; i<=row;i++)
772 {
773 t1 = view(i,j);
774 a->set(i-1,t1);
775 }
776}
777
778void bigintmat::getColRange(int j, int no, bigintmat *a)
779{
780 number t1;
781 for(int ii=0; ii< no; ii++)
782 {
783 for (int i=1; i<=row;i++)
784 {
785 t1 = view(i, ii+j);
786 a->set(i, ii+1, t1);
787 }
788 }
789}
790
792{
793 if ((i>row) || (i<1))
794 {
795 WerrorS("Error in getrow: Index out of range!");
796 return;
797 }
798 if (((a->rows() != 1) || (a->cols() != col)) && ((a->rows() != col) || (a->cols() != 1)))
799 {
800 WerrorS("Error in getrow. Dimensions must agree!");
801 return;
802 }
804 {
806 number t1, t2;
807 for (int j=1; j<=col;j++)
808 {
809 t1 = get(i,j);
810 t2 = f(t1, basecoeffs(), a->basecoeffs());
811 a->set(j-1,t2);
812 n_Delete(&t1, basecoeffs());
813 n_Delete(&t2, a->basecoeffs());
814 }
815 return;
816 }
817 number t1;
818 for (int j=1; j<=col;j++)
819 {
820 t1 = get(i,j);
821 a->set(j-1,t1);
822 n_Delete(&t1, basecoeffs());
823 }
824}
825
827{
828 if ((j>col) || (j<1))
829 {
830 WerrorS("Error in setcol: Index out of range!");
831 return;
832 }
833 if (((m->rows() != row) || (m->cols() != 1)) && ((m->rows() != 1) || (m->cols() != row)))
834 {
835 WerrorS("Error in setcol. Dimensions must agree!");
836 return;
837 }
838 if (!nCoeffs_are_equal(basecoeffs(), m->basecoeffs()))
839 {
840 nMapFunc f = n_SetMap(m->basecoeffs(), basecoeffs());
841 number t1,t2;
842 for (int i=1; i<=row; i++)
843 {
844 t1 = m->get(i-1);
845 t2 = f(t1, m->basecoeffs(),basecoeffs());
846 set(i, j, t2);
847 n_Delete(&t2, basecoeffs());
848 n_Delete(&t1, m->basecoeffs());
849 }
850 return;
851 }
852 number t1;
853 for (int i=1; i<=row; i++)
854 {
855 t1 = m->view(i-1);
856 set(i, j, t1);
857 }
858}
859
861{
862 if ((j>row) || (j<1))
863 {
864 WerrorS("Error in setrow: Index out of range!");
865 return;
866 }
867 if (((m->rows() != 1) || (m->cols() != col)) && ((m->rows() != col) || (m->cols() != 1)))
868 {
869 WerrorS("Error in setrow. Dimensions must agree!");
870 return;
871 }
872 if (!nCoeffs_are_equal(basecoeffs(), m->basecoeffs()))
873 {
874 nMapFunc f = n_SetMap(m->basecoeffs(), basecoeffs());
875 number tmp1,tmp2;
876 for (int i=1; i<=col; i++)
877 {
878 tmp1 = m->get(i-1);
879 tmp2 = f(tmp1, m->basecoeffs(),basecoeffs());
880 set(j, i, tmp2);
882 n_Delete(&tmp1, m->basecoeffs());
883 }
884 return;
885 }
886 number tmp;
887 for (int i=1; i<=col; i++)
888 {
889 tmp = m->view(i-1);
890 set(j, i, tmp);
891 }
892}
893
895{
896 if ((b->rows() != row) || (b->cols() != col))
897 {
898 WerrorS("Error in bigintmat::add. Dimensions do not agree!");
899 return false;
900 }
901 if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
902 {
903 WerrorS("Error in bigintmat::add. coeffs do not agree!");
904 return false;
905 }
906 for (int i=1; i<=row; i++)
907 {
908 for (int j=1; j<=col; j++)
909 {
910 rawset(i, j, n_Add(b->view(i,j), view(i,j), basecoeffs()));
911 }
912 }
913 return true;
914}
915
917{
918 if ((b->rows() != row) || (b->cols() != col))
919 {
920 WerrorS("Error in bigintmat::sub. Dimensions do not agree!");
921 return false;
922 }
923 if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
924 {
925 WerrorS("Error in bigintmat::sub. coeffs do not agree!");
926 return false;
927 }
928 for (int i=1; i<=row; i++)
929 {
930 for (int j=1; j<=col; j++)
931 {
932 rawset(i, j, n_Sub(view(i,j), b->view(i,j), basecoeffs()));
933 }
934 }
935 return true;
936}
937
939{
940 if (!nCoeffs_are_equal(c, basecoeffs()))
941 {
942 WerrorS("Wrong coeffs\n");
943 return false;
944 }
945 number t1, t2;
946 if ( n_IsOne(b,c)) return true;
947 for (int i=1; i<=row; i++)
948 {
949 for (int j=1; j<=col; j++)
950 {
951 t1 = view(i, j);
952 t2 = n_Mult(t1, b, basecoeffs());
953 rawset(i, j, t2);
954 }
955 }
956 return true;
957}
958
959bool bigintmat::addcol(int i, int j, number a, coeffs c)
960{
961 if ((i>col) || (j>col) || (i<1) || (j<1))
962 {
963 WerrorS("Error in addcol: Index out of range!");
964 return false;
965 }
966 if (!nCoeffs_are_equal(c, basecoeffs()))
967 {
968 WerrorS("Error in addcol: coeffs do not agree!");
969 return false;
970 }
971 number t1, t2, t3;
972 for (int k=1; k<=row; k++)
973 {
974 t1 = view(k, j);
975 t2 = view(k, i);
976 t3 = n_Mult(t1, a, basecoeffs());
977 n_InpAdd(t3, t2, basecoeffs());
978 rawset(k, i, t3);
979 }
980 return true;
981}
982
983bool bigintmat::addrow(int i, int j, number a, coeffs c)
984{
985 if ((i>row) || (j>row) || (i<1) || (j<1))
986 {
987 WerrorS("Error in addrow: Index out of range!");
988 return false;
989 }
990 if (!nCoeffs_are_equal(c, basecoeffs()))
991 {
992 WerrorS("Error in addrow: coeffs do not agree!");
993 return false;
994 }
995 number t1, t2, t3;
996 for (int k=1; k<=col; k++)
997 {
998 t1 = view(j, k);
999 t2 = view(i, k);
1000 t3 = n_Mult(t1, a, basecoeffs());
1001 n_InpAdd(t3, t2, basecoeffs());
1002 rawset(i, k, t3);
1003 }
1004 return true;
1005}
1006
1007void bigintmat::colskalmult(int i, number a, coeffs c)
1008{
1009 if ((i>=1) && (i<=col) && (nCoeffs_are_equal(c, basecoeffs())))
1010 {
1011 number t, tmult;
1012 for (int j=1; j<=row; j++)
1013 {
1014 t = view(j, i);
1015 tmult = n_Mult(a, t, basecoeffs());
1016 rawset(j, i, tmult);
1017 }
1018 }
1019 else
1020 WerrorS("Error in colskalmult");
1021}
1022
1023void bigintmat::rowskalmult(int i, number a, coeffs c)
1024{
1025 if ((i>=1) && (i<=row) && (nCoeffs_are_equal(c, basecoeffs())))
1026 {
1027 number t, tmult;
1028 for (int j=1; j<=col; j++)
1029 {
1030 t = view(i, j);
1031 tmult = n_Mult(a, t, basecoeffs());
1032 rawset(i, j, tmult);
1033 }
1034 }
1035 else
1036 WerrorS("Error in rowskalmult");
1037}
1038
1040{
1041 int ay = a->cols();
1042 int ax = a->rows();
1043 int by = b->cols();
1044 int bx = b->rows();
1045 number tmp;
1046 if (!((col == ay) && (col == by) && (ax+bx == row)))
1047 {
1048 WerrorS("Error in concatrow. Dimensions must agree!");
1049 return;
1050 }
1051 if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
1052 {
1053 WerrorS("Error in concatrow. coeffs do not agree!");
1054 return;
1055 }
1056 for (int i=1; i<=ax; i++)
1057 {
1058 for (int j=1; j<=ay; j++)
1059 {
1060 tmp = a->get(i,j);
1061 set(i, j, tmp);
1062 n_Delete(&tmp, basecoeffs());
1063 }
1064 }
1065 for (int i=1; i<=bx; i++)
1066 {
1067 for (int j=1; j<=by; j++)
1068 {
1069 tmp = b->get(i,j);
1070 set(i+ax, j, tmp);
1071 n_Delete(&tmp, basecoeffs());
1072 }
1073 }
1074}
1075
1077{
1078 bigintmat * tmp = new bigintmat(rows(), i, basecoeffs());
1079 appendCol(tmp);
1080 delete tmp;
1081}
1082
1084{
1085 coeffs R = basecoeffs();
1086 int ay = a->cols();
1087 int ax = a->rows();
1088 assume(row == ax);
1089
1091
1092 bigintmat * tmp = new bigintmat(rows(), cols() + ay, R);
1093 tmp->concatcol(this, a);
1094 this->swapMatrix(tmp);
1095 delete tmp;
1096}
1097
1099 int ay = a->cols();
1100 int ax = a->rows();
1101 int by = b->cols();
1102 int bx = b->rows();
1103 number tmp;
1104
1105 assume(row==ax && row == bx && ay+by ==col);
1106
1108
1109 for (int i=1; i<=ax; i++)
1110 {
1111 for (int j=1; j<=ay; j++)
1112 {
1113 tmp = a->view(i,j);
1114 set(i, j, tmp);
1115 }
1116 }
1117 for (int i=1; i<=bx; i++)
1118 {
1119 for (int j=1; j<=by; j++)
1120 {
1121 tmp = b->view(i,j);
1122 set(i, j+ay, tmp);
1123 }
1124 }
1125}
1126
1128{
1129 int ay = a->cols();
1130 int ax = a->rows();
1131 int by = b->cols();
1132 int bx = b->rows();
1133 number tmp;
1134 if (!(ax + bx == row))
1135 {
1136 WerrorS("Error in splitrow. Dimensions must agree!");
1137 }
1138 else if (!((col == ay) && (col == by)))
1139 {
1140 WerrorS("Error in splitrow. Dimensions must agree!");
1141 }
1142 else if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
1143 {
1144 WerrorS("Error in splitrow. coeffs do not agree!");
1145 }
1146 else
1147 {
1148 for(int i = 1; i<=ax; i++)
1149 {
1150 for(int j = 1; j<=ay;j++)
1151 {
1152 tmp = get(i,j);
1153 a->set(i,j,tmp);
1154 n_Delete(&tmp, basecoeffs());
1155 }
1156 }
1157 for (int i =1; i<=bx; i++)
1158 {
1159 for (int j=1;j<=col;j++)
1160 {
1161 tmp = get(i+ax, j);
1162 b->set(i,j,tmp);
1163 n_Delete(&tmp, basecoeffs());
1164 }
1165 }
1166 }
1167}
1168
1170{
1171 int ay = a->cols();
1172 int ax = a->rows();
1173 int by = b->cols();
1174 int bx = b->rows();
1175 number tmp;
1176 if (!((row == ax) && (row == bx)))
1177 {
1178 WerrorS("Error in splitcol. Dimensions must agree!");
1179 }
1180 else if (!(ay+by == col))
1181 {
1182 WerrorS("Error in splitcol. Dimensions must agree!");
1183 }
1184 else if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs()) && nCoeffs_are_equal(b->basecoeffs(), basecoeffs())))
1185 {
1186 WerrorS("Error in splitcol. coeffs do not agree!");
1187 }
1188 else
1189 {
1190 for (int i=1; i<=ax; i++)
1191 {
1192 for (int j=1; j<=ay; j++)
1193 {
1194 tmp = view(i,j);
1195 a->set(i,j,tmp);
1196 }
1197 }
1198 for (int i=1; i<=bx; i++)
1199 {
1200 for (int j=1; j<=by; j++)
1201 {
1202 tmp = view(i,j+ay);
1203 b->set(i,j,tmp);
1204 }
1205 }
1206 }
1207}
1208
1210{
1211 number tmp;
1212 if ((a->rows() != row) || (a->cols()+i-1 > col) || (i<1))
1213 {
1214 WerrorS("Error in splitcol. Dimensions must agree!");
1215 return;
1216 }
1217 if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs())))
1218 {
1219 WerrorS("Error in splitcol. coeffs do not agree!");
1220 return;
1221 }
1222 int width = a->cols();
1223 for (int j=1; j<=width; j++)
1224 {
1225 for (int k=1; k<=row; k++)
1226 {
1227 tmp = get(k, j+i-1);
1228 a->set(k, j, tmp);
1229 n_Delete(&tmp, basecoeffs());
1230 }
1231 }
1232}
1233
1235{
1236 number tmp;
1237 if ((a->cols() != col) || (a->rows()+i-1 > row) || (i<1))
1238 {
1239 WerrorS("Error in Marco-splitrow");
1240 return;
1241 }
1242
1243 if (!(nCoeffs_are_equal(a->basecoeffs(), basecoeffs())))
1244 {
1245 WerrorS("Error in splitrow. coeffs do not agree!");
1246 return;
1247 }
1248 int height = a->rows();
1249 for (int j=1; j<=height; j++)
1250 {
1251 for (int k=1; k<=col; k++)
1252 {
1253 tmp = view(j+i-1, k);
1254 a->set(j, k, tmp);
1255 }
1256 }
1257}
1258
1260{
1261 if ((b->rows() != row) || (b->cols() != col))
1262 {
1263 WerrorS("Error in bigintmat::copy. Dimensions do not agree!");
1264 return false;
1265 }
1266 if (!nCoeffs_are_equal(basecoeffs(), b->basecoeffs()))
1267 {
1268 WerrorS("Error in bigintmat::copy. coeffs do not agree!");
1269 return false;
1270 }
1271 number t1;
1272 for (int i=1; i<=row; i++)
1273 {
1274 for (int j=1; j<=col; j++)
1275 {
1276 t1 = b->view(i, j);
1277 set(i, j, t1);
1278 }
1279 }
1280 return true;
1281}
1282
1283/// copy the submatrix of b, staring at (a,b) having n rows, m cols into
1284/// the given matrix at pos. (c,d)
1285/// needs c+n, d+m <= rows, cols
1286/// a+n, b+m <= b.rows(), b.cols()
1287void bigintmat::copySubmatInto(bigintmat *B, int a, int b, int n, int m, int c, int d)
1288{
1289 number t1;
1290 for (int i=1; i<=n; i++)
1291 {
1292 for (int j=1; j<=m; j++)
1293 {
1294 t1 = B->view(a+i-1, b+j-1);
1295 set(c+i-1, d+j-1, t1);
1296 }
1297 }
1298}
1299
1301{
1302 coeffs r = basecoeffs();
1303 if (row==col)
1304 {
1305 for (int i=1; i<=row; i++)
1306 {
1307 for (int j=1; j<=col; j++)
1308 {
1309 if (i==j)
1310 {
1311 if (!n_IsOne(view(i, j), r))
1312 return 0;
1313 }
1314 else
1315 {
1316 if (!n_IsZero(view(i,j), r))
1317 return 0;
1318 }
1319 }
1320 }
1321 }
1322 return 1;
1323}
1324
1326{
1327 if (row==col)
1328 {
1329 number one = n_Init(1, basecoeffs()),
1330 zero = n_Init(0, basecoeffs());
1331 for (int i=1; i<=row; i++)
1332 {
1333 for (int j=1; j<=col; j++)
1334 {
1335 if (i==j)
1336 {
1337 set(i, j, one);
1338 }
1339 else
1340 {
1341 set(i, j, zero);
1342 }
1343 }
1344 }
1345 n_Delete(&one, basecoeffs());
1347 }
1348}
1349
1351{
1352 number tmp = n_Init(0, basecoeffs());
1353 for (int i=1; i<=row; i++)
1354 {
1355 for (int j=1; j<=col; j++)
1356 {
1357 set(i, j, tmp);
1358 }
1359 }
1360 n_Delete(&tmp,basecoeffs());
1361}
1362
1364{
1365 for (int i=1; i<=row; i++) {
1366 for (int j=1; j<=col; j++) {
1367 if (!n_IsZero(view(i,j), basecoeffs()))
1368 return FALSE;
1369 }
1370 }
1371 return TRUE;
1372}
1373//****************************************************************************
1374//
1375//****************************************************************************
1376
1377//used in the det function. No idea what it does.
1378//looks like it return the submatrix where the i-th row
1379//and j-th column has been removed in the LaPlace generic
1380//determinant algorithm
1382{
1383 if ((i<=0) || (i>row) || (j<=0) || (j>col))
1384 return NULL;
1385 int cx, cy;
1386 cx=1;
1387 cy=1;
1388 number t;
1389 bigintmat *b = new bigintmat(row-1, col-1, basecoeffs());
1390 for (int k=1; k<=row; k++) {
1391 if (k!=i)
1392 {
1393 cy=1;
1394 for (int l=1; l<=col; l++)
1395 {
1396 if (l!=j)
1397 {
1398 t = get(k, l);
1399 b->set(cx, cy, t);
1400 n_Delete(&t, basecoeffs());
1401 cy++;
1402 }
1403 }
1404 cx++;
1405 }
1406 }
1407 return b;
1408}
1409
1410
1411//returns d such that a/d is the inverse of the input
1412//TODO: make work for p not prime using the euc stuff.
1413//long term: rewrite for Z using p-adic lifting
1414//and Dixon. Possibly even the sparse recent Storjohann stuff
1416
1417 // Falls Matrix über reellen Zahlen nicht invertierbar, breche ab
1418 assume((a->rows() == row) && (a->rows() == a->cols()) && (row == col));
1419
1420 number det = this->det(); //computes the HNF, so should e reused.
1421 if ((n_IsZero(det, basecoeffs())))
1422 return det;
1423
1424 // Hänge Einheitsmatrix über Matrix und wendet HNF an. An Stelle der Einheitsmatrix steht im Ergebnis die Transformationsmatrix dazu
1425 a->one();
1426 bigintmat *m = new bigintmat(2*row, col, basecoeffs());
1427 m->concatrow(a,this);
1428 m->hnf();
1429 // Arbeite weiterhin mit der zusammengehängten Matrix
1430 // 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
1431 number diag;
1432 number temp, ttemp;
1433 for (int i=1; i<=col; i++) {
1434 diag = m->get(row+i, i);
1435 for (int j=i+1; j<=col; j++) {
1436 temp = m->get(row+i, j);
1437 m->colskalmult(j, diag, basecoeffs());
1438 temp = n_InpNeg(temp, basecoeffs());
1439 m->addcol(j, i, temp, basecoeffs());
1440 n_Delete(&temp, basecoeffs());
1441 }
1442 n_Delete(&diag, basecoeffs());
1443 }
1444 // Falls wir nicht modulo n arbeiten, können wir die Spalten durch den ggT teilen, um die Einträge kleiner zu bekommen
1445 // Bei Z/n sparen wir uns das, da es hier sinnlos ist
1446 number g;
1447 number gcd;
1448 for (int j=1; j<=col; j++) {
1449 g = n_Init(0, basecoeffs());
1450 for (int i=1; i<=2*row; i++) {
1451 temp = m->get(i,j);
1452 gcd = n_Gcd(g, temp, basecoeffs());
1453 n_Delete(&g, basecoeffs());
1454 n_Delete(&temp, basecoeffs());
1455 g = n_Copy(gcd, basecoeffs());
1456 n_Delete(&gcd, basecoeffs());
1457 }
1458 if (!(n_IsOne(g, basecoeffs())))
1459 m->colskaldiv(j, g);
1460 n_Delete(&g, basecoeffs());
1461 }
1462
1463 // 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)
1464
1465 g = n_Init(0, basecoeffs());
1466 number prod = n_Init(1, basecoeffs());
1467 for (int i=1; i<=col; i++) {
1468 gcd = n_Gcd(g, m->get(row+i, i), basecoeffs());
1469 n_Delete(&g, basecoeffs());
1470 g = n_Copy(gcd, basecoeffs());
1471 n_Delete(&gcd, basecoeffs());
1472 ttemp = n_Copy(prod, basecoeffs());
1473 temp = m->get(row+i, i);
1475 prod = n_Mult(ttemp, temp, basecoeffs());
1476 n_Delete(&ttemp, basecoeffs());
1477 n_Delete(&temp, basecoeffs());
1478 }
1479 number lcm;
1480 lcm = n_Div(prod, g, basecoeffs());
1481 for (int j=1; j<=col; j++) {
1482 ttemp = m->get(row+j,j);
1483 temp = n_QuotRem(lcm, ttemp, NULL, basecoeffs());
1484 m->colskalmult(j, temp, basecoeffs());
1485 n_Delete(&ttemp, basecoeffs());
1486 n_Delete(&temp, basecoeffs());
1487 }
1488 n_Delete(&lcm, basecoeffs());
1490
1491 number divisor = m->get(row+1, 1);
1492 m->splitrow(a, 1);
1493 delete m;
1494 n_Delete(&det, basecoeffs());
1495 return divisor;
1496}
1497
1499{
1500 assume (col == row);
1501 number t = get(1,1),
1502 h;
1503 coeffs r = basecoeffs();
1504 for(int i=2; i<= col; i++) {
1505 h = n_Add(t, view(i,i), r);
1506 n_Delete(&t, r);
1507 t = h;
1508 }
1509 return t;
1510}
1511
1513{
1514 assume (row==col);
1515
1516 if (col == 1)
1517 return get(1, 1);
1518 // should work as well in Z/pZ of type n_Zp?
1519 // relies on XExtGcd and the other euc. functinos.
1521 return hnfdet();
1522 }
1523 number sum = n_Init(0, basecoeffs());
1524 number t1, t2, t3, t4;
1525 bigintmat *b;
1526 for (int i=1; i<=row; i++) {
1527 b = elim(i, 1);
1528 t1 = get(i, 1);
1529 t2 = b->det();
1530 t3 = n_Mult(t1, t2, basecoeffs());
1531 t4 = n_Copy(sum, basecoeffs());
1532 n_Delete(&sum, basecoeffs());
1533 if ((i+1)>>1<<1==(i+1))
1534 sum = n_Add(t4, t3, basecoeffs());
1535 else
1536 sum = n_Sub(t4, t3, basecoeffs());
1537 n_Delete(&t1, basecoeffs());
1538 n_Delete(&t2, basecoeffs());
1539 n_Delete(&t3, basecoeffs());
1540 n_Delete(&t4, basecoeffs());
1541 }
1542 return sum;
1543}
1544
1546{
1547 assume (col == row);
1548
1549 if (col == 1)
1550 return get(1, 1);
1551 bigintmat *m = new bigintmat(this);
1552 m->hnf();
1553 number prod = n_Init(1, basecoeffs());
1554 number temp, temp2;
1555 for (int i=1; i<=col; i++) {
1556 temp = m->get(i, i);
1557 temp2 = n_Mult(temp, prod, basecoeffs());
1559 prod = temp2;
1560 n_Delete(&temp, basecoeffs());
1561 }
1562 delete m;
1563 return prod;
1564}
1565
1567{
1568 int n = rows(), m = cols();
1569 row = a->rows();
1570 col = a->cols();
1571 number * V = v;
1572 v = a->v;
1573 a->v = V;
1574 a->row = n;
1575 a->col = m;
1576}
1578{
1579 coeffs R = basecoeffs();
1580 for(int i=1; i<=rows(); i++)
1581 if (!n_IsZero(view(i, j), R)) return FALSE;
1582 return TRUE;
1583}
1584
1586{
1587 coeffs R = basecoeffs();
1588 hnf(); // as a starting point...
1589 if (getCoeffType(R)== n_Z) return; //wrong, need to prune!
1590
1591 int n = cols(), m = rows(), i, j, k;
1592
1593 //make sure, the matrix has enough space. We need no rows+1 columns.
1594 //The resulting Howell form will be pruned to be at most square.
1595 bigintmat * t = new bigintmat(m, m+1, R);
1596 t->copySubmatInto(this, 1, n>m ? n-m+1 : 1, m, n>m ? m : n, 1, n>m ? 2 : m+2-n );
1597 swapMatrix(t);
1598 delete t;
1599 for(i=1; i<= cols(); i++) {
1600 if (!colIsZero(i)) break;
1601 }
1602 assume (i>1);
1603 if (i>cols()) {
1604 t = new bigintmat(rows(), 0, R);
1605 swapMatrix(t);
1606 delete t;
1607 return; // zero matrix found, clearly normal.
1608 }
1609
1610 int last_zero_col = i-1;
1611 for (int c = cols(); c>0; c--) {
1612 for(i=rows(); i>0; i--) {
1613 if (!n_IsZero(view(i, c), R)) break;
1614 }
1615 if (i==0) break; // matrix SHOULD be zero from here on
1616 number a = n_Ann(view(i, c), R);
1617 addcol(last_zero_col, c, a, R);
1618 n_Delete(&a, R);
1619 for(j = c-1; j>last_zero_col; j--) {
1620 for(k=rows(); k>0; k--) {
1621 if (!n_IsZero(view(k, j), R)) break;
1622 if (!n_IsZero(view(k, last_zero_col), R)) break;
1623 }
1624 if (k==0) break;
1625 if (!n_IsZero(view(k, last_zero_col), R)) {
1626 number gcd, co1, co2, co3, co4;
1627 gcd = n_XExtGcd(view(k, last_zero_col), view(k, j), &co1, &co2, &co3, &co4, R);
1628 if (n_Equal(gcd, view(k, j), R)) {
1629 number q = n_Div(view(k, last_zero_col), gcd, R);
1630 q = n_InpNeg(q, R);
1631 addcol(last_zero_col, j, q, R);
1632 n_Delete(&q, R);
1633 } else if (n_Equal(gcd, view(k, last_zero_col), R)) {
1634 swap(last_zero_col, k);
1635 number q = n_Div(view(k, last_zero_col), gcd, R);
1636 q = n_InpNeg(q, R);
1637 addcol(last_zero_col, j, q, R);
1638 n_Delete(&q, R);
1639 } else {
1640 coltransform(last_zero_col, j, co3, co4, co1, co2);
1641 }
1642 n_Delete(&gcd, R);
1643 n_Delete(&co1, R);
1644 n_Delete(&co2, R);
1645 n_Delete(&co3, R);
1646 n_Delete(&co4, R);
1647 }
1648 }
1649 for(k=rows(); k>0; k--) {
1650 if (!n_IsZero(view(k, last_zero_col), R)) break;
1651 }
1652 if (k) last_zero_col--;
1653 }
1654 t = new bigintmat(rows(), cols()-last_zero_col, R);
1655 t->copySubmatInto(this, 1, last_zero_col+1, rows(), cols()-last_zero_col, 1, 1);
1656 swapMatrix(t);
1657 delete t;
1658}
1659
1661{
1662 // Laufen von unten nach oben und von links nach rechts
1663 // CF: TODO: for n_Z: write a recursive version. This one will
1664 // have exponential blow-up. Look at Michianchio
1665 // Alternatively, do p-adic det and modular method
1666
1667#if 0
1668 char * s;
1669 ::PrintS("mat over Z is \n");
1670 ::Print("%s\n", s = nCoeffString(basecoeffs()));
1671 omFree(s);
1672 Print();
1673 ::Print("\n(%d x %d)\n", rows(), cols());
1674#endif
1675
1676 int i = rows();
1677 int j = cols();
1678 number q = n_Init(0, basecoeffs());
1679 number one = n_Init(1, basecoeffs());
1680 number minusone = n_Init(-1, basecoeffs());
1681 number tmp1 = n_Init(0, basecoeffs());
1682 number tmp2 = n_Init(0, basecoeffs());
1683 number co1, co2, co3, co4;
1684 number ggt = n_Init(0, basecoeffs());
1685
1686 while ((i>0) && (j>0))
1687 {
1688 // Falls erstes Nicht-Null-Element in Zeile i nicht existiert, oder hinter Spalte j vorkommt, gehe in nächste Zeile
1689 if ((findnonzero(i)==0) || (findnonzero(i)>j))
1690 {
1691 i--;
1692 }
1693 else
1694 {
1695 // Laufe von links nach rechts durch die Zeile:
1696 for (int l=1; l<=j-1; l++)
1697 {
1699 tmp1 = get(i, l);
1700 // Falls Eintrag (im folgenden x genannt) gleich 0, gehe eine Spalte weiter. Ansonsten...
1701 if (!n_IsZero(tmp1, basecoeffs()))
1702 {
1704 tmp2 = get(i, l+1);
1705 // Falls Eintrag (i.f. y g.) rechts daneben gleich 0, tausche beide Spalten, sonst...
1706 if (!n_IsZero(tmp2, basecoeffs()))
1707 {
1708 n_Delete(&ggt, basecoeffs());
1709 ggt = n_XExtGcd(tmp1, tmp2, &co1, &co2, &co3, &co4, basecoeffs());
1710 // 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:
1711 if (n_Equal(tmp1, ggt, basecoeffs()))
1712 {
1713 swap(l, l+1);
1714 n_Delete(&q, basecoeffs());
1715 q = n_Div(tmp2, ggt, basecoeffs());
1716 q = n_InpNeg(q, basecoeffs());
1717 // Dann addiere das -q-fache der (neuen) rechten Spalte zur linken dazu. Damit erhalten wir die gewünschte 0
1718
1719 addcol(l, l+1, q, basecoeffs());
1720 n_Delete(&q, basecoeffs());
1721 }
1722 else if (n_Equal(tmp1, minusone, basecoeffs()))
1723 {
1724 // 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
1725 // Die Berechnung von q (=y/ggt) entfällt, da ggt=1
1726 swap(l, l+1);
1727 colskalmult(l+1, minusone, basecoeffs());
1729 addcol(l, l+1, tmp2, basecoeffs());
1730 }
1731 else
1732 {
1733 // CF: use the 2x2 matrix (co1, co2)(co3, co4) to
1734 // get the gcd in position and the 0 in the other:
1735#ifdef CF_DEB
1736 ::PrintS("applying trafo\n");
1737 StringSetS("");
1738 n_Write(co1, basecoeffs()); StringAppendS("\t");
1739 n_Write(co2, basecoeffs()); StringAppendS("\t");
1740 n_Write(co3, basecoeffs()); StringAppendS("\t");
1741 n_Write(co4, basecoeffs()); StringAppendS("\t");
1742 ::Print("%s\nfor l=%d\n", StringEndS(), l);
1743 {char * s = String();
1744 ::Print("to %s\n", s);omFree(s);};
1745#endif
1746 coltransform(l, l+1, co3, co4, co1, co2);
1747#ifdef CF_DEB
1748 {char * s = String();
1749 ::Print("gives %s\n", s);}
1750#endif
1751 }
1752 n_Delete(&co1, basecoeffs());
1753 n_Delete(&co2, basecoeffs());
1754 n_Delete(&co3, basecoeffs());
1755 n_Delete(&co4, basecoeffs());
1756 }
1757 else
1758 {
1759 swap(l, l+1);
1760 }
1761 // Dann betrachte die vormals rechte Spalte als neue linke, und die rechts daneben als neue rechte.
1762 }
1763 }
1764
1765 #ifdef HAVE_RINGS
1766 // normalize by units:
1767 if (!n_IsZero(view(i, j), basecoeffs()))
1768 {
1769 number u = n_GetUnit(view(i, j), basecoeffs());
1770 if (!n_IsOne(u, basecoeffs()))
1771 {
1772 colskaldiv(j, u);
1773 }
1774 n_Delete(&u, basecoeffs());
1775 }
1776 #endif
1777 // Zum Schluss mache alle Einträge rechts vom Diagonalelement betragsmäßig kleiner als dieses
1778 for (int l=j+1; l<=col; l++)
1779 {
1780 n_Delete(&q, basecoeffs());
1781 q = n_QuotRem(view(i, l), view(i, j), NULL, basecoeffs());
1782 q = n_InpNeg(q, basecoeffs());
1783 addcol(l, j, q, basecoeffs());
1784 }
1785 i--;
1786 j--;
1787 // Dann betrachte die Zeile darüber und gehe dort wie vorher vor
1788 }
1789 }
1790 n_Delete(&q, basecoeffs());
1793 n_Delete(&ggt, basecoeffs());
1794 n_Delete(&one, basecoeffs());
1795 n_Delete(&minusone, basecoeffs());
1796
1797#if 0
1798 ::PrintS("hnf over Z is \n");
1799 Print();
1800 ::Print("\n(%d x %d)\n", rows(), cols());
1801#endif
1802}
1803
1805{
1806 coeffs cold = a->basecoeffs();
1807 bigintmat *b = new bigintmat(a->rows(), a->cols(), cnew);
1808 // Erzeugt Karte von alten coeffs nach neuen
1809 nMapFunc f = n_SetMap(cold, cnew);
1810 number t1;
1811 number t2;
1812 // apply map to all entries.
1813 for (int i=1; i<=a->rows(); i++)
1814 {
1815 for (int j=1; j<=a->cols(); j++)
1816 {
1817 t1 = a->get(i, j);
1818 t2 = f(t1, cold, cnew);
1819 b->set(i, j, t2);
1820 n_Delete(&t1, cold);
1821 n_Delete(&t2, cnew);
1822 }
1823 }
1824 return b;
1825}
1826
1827#ifdef HAVE_RINGS
1828//OK: a HNF of (this | p*I)
1829//so the result will always have FULL rank!!!!
1830//(This is different form a lift of the HNF mod p: consider the matrix (p)
1831//to see the difference. It CAN be computed as HNF mod p^2 usually..)
1833{
1834 coeffs Rp = numbercoeffs(p, R); // R/pR
1835 bigintmat *m = bimChangeCoeff(this, Rp);
1836 m->howell();
1837 bigintmat *a = bimChangeCoeff(m, R);
1838 delete m;
1839 bigintmat *C = new bigintmat(rows(), rows(), R);
1840 int piv = rows(), i = a->cols();
1841 while (piv)
1842 {
1843 if (!i || n_IsZero(a->view(piv, i), R))
1844 {
1845 C->set(piv, piv, p, R);
1846 }
1847 else
1848 {
1849 C->copySubmatInto(a, 1, i, rows(), 1, 1, piv);
1850 i--;
1851 }
1852 piv--;
1853 }
1854 delete a;
1855 return C;
1856}
1857#endif
1858
1859
1860//exactly divide matrix by b
1862{
1863 number tmp1, tmp2;
1864 for (int i=1; i<=row; i++)
1865 {
1866 for (int j=1; j<=col; j++)
1867 {
1868 tmp1 = view(i, j);
1869 tmp2 = n_Div(tmp1, b, basecoeffs());
1870 rawset(i, j, tmp2);
1871 }
1872 }
1873}
1874
1875//exactly divide col j by b
1876void bigintmat::colskaldiv(int j, number b)
1877{
1878 number tmp1, tmp2;
1879 for (int i=1; i<=row; i++)
1880 {
1881 tmp1 = view(i, j);
1882 tmp2 = n_Div(tmp1, b, basecoeffs());
1883 rawset(i, j, tmp2);
1884 }
1885}
1886
1887// col(j, k) <- col(j,k)*matrix((a, c)(b, d))
1888// mostly used internally in the hnf and Howell stuff
1889void bigintmat::coltransform(int j, int k, number a, number b, number c, number d)
1890{
1891 number tmp1, tmp2, tmp3, tmp4;
1892 for (int i=1; i<=row; i++)
1893 {
1894 tmp1 = get(i, j);
1895 tmp2 = get(i, k);
1896 tmp3 = n_Mult(tmp1, a, basecoeffs());
1897 tmp4 = n_Mult(tmp2, b, basecoeffs());
1898 n_InpAdd(tmp3, tmp4, basecoeffs());
1899 n_Delete(&tmp4, basecoeffs());
1900
1901 n_InpMult(tmp1, c, basecoeffs());
1902 n_InpMult(tmp2, d, basecoeffs());
1905
1906 set(i, j, tmp3);
1907 set(i, k, tmp1);
1909 n_Delete(&tmp3, basecoeffs());
1910 }
1911}
1912
1913
1914
1915//reduce all entries mod p. Does NOT change the coeffs type
1916void bigintmat::mod(number p)
1917{
1918 // produce the matrix in Z/pZ
1919 number tmp1, tmp2;
1920 for (int i=1; i<=row; i++)
1921 {
1922 for (int j=1; j<=col; j++)
1923 {
1924 tmp1 = get(i, j);
1925 tmp2 = n_IntMod(tmp1, p, basecoeffs());
1927 set(i, j, tmp2);
1928 }
1929 }
1930}
1931
1933{
1934 if (!nCoeffs_are_equal(a->basecoeffs(), b->basecoeffs()))
1935 {
1936 WerrorS("Error in bimMult. Coeffs do not agree!");
1937 return;
1938 }
1939 if ((a->rows() != c->rows()) || (b->cols() != c->cols()) || (a->cols() != b->rows()))
1940 {
1941 WerrorS("Error in bimMult. Dimensions do not agree!");
1942 return;
1943 }
1944 bigintmat *tmp = bimMult(a, b);
1945 c->copy(tmp);
1946
1947 delete tmp;
1948}
1949
1951{
1952 //write b = Ax + eps where eps is "small" in the sense of bounded by the
1953 //pivot entries in H. H does not need to be Howell (or HNF) but need
1954 //to be triagonal in the same direction.
1955 //b can have multiple columns.
1956#if 0
1957 PrintS("reduce_mod_howell: A:\n");
1958 A->Print();
1959 PrintS("\nb:\n");
1960 b->Print();
1961#endif
1962
1963 coeffs R = A->basecoeffs();
1964 assume(x->basecoeffs() == R);
1965 assume(b->basecoeffs() == R);
1966 assume(eps->basecoeffs() == R);
1967 if (!A->cols())
1968 {
1969 x->zero();
1970 eps->copy(b);
1971
1972#if 0
1973 PrintS("\nx:\n");
1974 x->Print();
1975 PrintS("\neps:\n");
1976 eps->Print();
1977 PrintS("\n****************************************\n");
1978#endif
1979 return;
1980 }
1981
1982 bigintmat * B = new bigintmat(b->rows(), 1, R);
1983 for(int i=1; i<= b->cols(); i++)
1984 {
1985 int A_col = A->cols();
1986 b->getcol(i, B);
1987 for(int j = B->rows(); j>0; j--)
1988 {
1989 number Ai = A->view(A->rows() - B->rows() + j, A_col);
1990 if (n_IsZero(Ai, R) &&
1991 n_IsZero(B->view(j, 1), R))
1992 {
1993 continue; //all is fine: 0*x = 0
1994 }
1995 else if (n_IsZero(B->view(j, 1), R))
1996 {
1997 x->rawset(x->rows() - B->rows() + j, i, n_Init(0, R));
1998 A_col--;
1999 }
2000 else if (n_IsZero(Ai, R))
2001 {
2002 A_col--;
2003 }
2004 else
2005 {
2006 // "solve" ax=b, possibly enlarging d
2007 number Bj = B->view(j, 1);
2008 number q = n_Div(Bj, Ai, R);
2009 x->rawset(x->rows() - B->rows() + j, i, q);
2010 for(int k=j; k>B->rows() - A->rows(); k--)
2011 {
2012 //B[k] = B[k] - x[k]A[k][j]
2013 number s = n_Mult(q, A->view(A->rows() - B->rows() + k, A_col), R);
2014 B->rawset(k, 1, n_Sub(B->view(k, 1), s, R));
2015 n_Delete(&s, R);
2016 }
2017 A_col--;
2018 }
2019 if (!A_col)
2020 {
2021 break;
2022 }
2023 }
2024 eps->setcol(i, B);
2025 }
2026 delete B;
2027#if 0
2028 PrintS("\nx:\n");
2029 x->Print();
2030 PrintS("\neps:\n");
2031 eps->Print();
2032 PrintS("\n****************************************\n");
2033#endif
2034}
2035
2037{
2038 coeffs R = A->basecoeffs();
2039 bigintmat *m = new bigintmat(A->rows()+A->cols(), A->cols(), R);
2040 m->copySubmatInto(A, 1, 1, A->rows(), A->cols(), A->cols()+1, 1);
2041 number one = n_Init(1, R);
2042 for(int i=1; i<= A->cols(); i++)
2043 m->set(i,i,one);
2044 n_Delete(&one, R);
2045 return m;
2046}
2047
2048static number bimFarey(bigintmat *A, number N, bigintmat *L)
2049{
2050 coeffs Z = A->basecoeffs(),
2051 Q = nInitChar(n_Q, 0);
2052 number den = n_Init(1, Z);
2053 nMapFunc f = n_SetMap(Q, Z);
2054
2055 for(int i=1; i<= A->rows(); i++)
2056 {
2057 for(int j=1; j<= A->cols(); j++)
2058 {
2059 number ad = n_Mult(den, A->view(i, j), Z);
2060 number re = n_IntMod(ad, N, Z);
2061 n_Delete(&ad, Z);
2062 number q = n_Farey(re, N, Z);
2063 n_Delete(&re, Z);
2064 if (!q)
2065 {
2066 n_Delete(&ad, Z);
2067 n_Delete(&den, Z);
2068 return NULL;
2069 }
2070
2071 number d = n_GetDenom(q, Q),
2072 n = n_GetNumerator(q, Q);
2073
2074 n_Delete(&q, Q);
2075 n_Delete(&ad, Z);
2076 number dz = f(d, Q, Z),
2077 nz = f(n, Q, Z);
2078 n_Delete(&d, Q);
2079 n_Delete(&n, Q);
2080
2081 if (!n_IsOne(dz, Z))
2082 {
2083 L->skalmult(dz, Z);
2084 n_InpMult(den, dz, Z);
2085#if 0
2086 PrintS("den increasing to ");
2087 n_Print(den, Z);
2088 PrintLn();
2089#endif
2090 }
2091 n_Delete(&dz, Z);
2092 L->rawset(i, j, nz);
2093 }
2094 }
2095
2096 nKillChar(Q);
2097 PrintS("bimFarey worked\n");
2098#if 0
2099 L->Print();
2100 PrintS("\n * 1/");
2101 n_Print(den, Z);
2102 PrintLn();
2103#endif
2104 return den;
2105}
2106
2107#ifdef HAVE_RINGS
2109 coeffs R = A->basecoeffs();
2110
2111 assume(getCoeffType(R) == n_Z);
2112
2113 number p = n_Init(536870909, R); // PreviousPrime(2^29); not clever
2114 coeffs Rp = numbercoeffs(p, R); // R/pR
2115 bigintmat *Ap = bimChangeCoeff(A, Rp),
2116 *m = prependIdentity(Ap),
2117 *Tp, *Hp;
2118 delete Ap;
2119
2120 m->howell();
2121 Hp = new bigintmat(A->rows(), A->cols(), Rp);
2122 Hp->copySubmatInto(m, A->cols()+1, 1, A->rows(), A->cols(), 1, 1);
2123 Tp = new bigintmat(A->cols(), A->cols(), Rp);
2124 Tp->copySubmatInto(m, 1, 1, A->cols(), A->cols(), 1, 1);
2125
2126 int i, j;
2127
2128 for(i=1; i<= A->cols(); i++)
2129 {
2130 for(j=m->rows(); j>A->cols(); j--)
2131 {
2132 if (!n_IsZero(m->view(j, i), Rp)) break;
2133 }
2134 if (j>A->cols()) break;
2135 }
2136// Print("Found nullity (kern dim) of %d\n", i-1);
2137 bigintmat * kp = new bigintmat(A->cols(), i-1, Rp);
2138 kp->copySubmatInto(Tp, 1, 1, A->cols(), i-1, 1, 1);
2139 kp->howell();
2140
2141 delete m;
2142
2143 //Hp is the mod-p howell form
2144 //Tp the transformation, mod p
2145 //kp a basis for the kernel, in howell form, mod p
2146
2147 bigintmat * eps_p = new bigintmat(B->rows(), B->cols(), Rp),
2148 * x_p = new bigintmat(A->cols(), B->cols(), Rp),
2149 * fps_p = new bigintmat(kp->cols(), B->cols(), Rp);
2150
2151 //initial solution
2152
2153 number zero = n_Init(0, R);
2154 x->skalmult(zero, R);
2155 n_Delete(&zero, R);
2156
2157 bigintmat * b = new bigintmat(B);
2158 number pp = n_Init(1, R);
2159 i = 1;
2160 do
2161 {
2162 bigintmat * b_p = bimChangeCoeff(b, Rp), * s;
2163 bigintmat * t1, *t2;
2164 reduce_mod_howell(Hp, b_p, eps_p, x_p);
2165 delete b_p;
2166 if (!eps_p->isZero())
2167 {
2168 PrintS("no solution, since no modular solution\n");
2169
2170 delete eps_p;
2171 delete x_p;
2172 delete Hp;
2173 delete kp;
2174 delete Tp;
2175 delete b;
2176 n_Delete(&pp, R);
2177 n_Delete(&p, R);
2178 nKillChar(Rp);
2179
2180 return NULL;
2181 }
2182 t1 = bimMult(Tp, x_p);
2183 delete x_p;
2184 x_p = t1;
2185 reduce_mod_howell(kp, x_p, x_p, fps_p); //we're not all interested in fps_p
2186 s = bimChangeCoeff(x_p, R);
2187 t1 = bimMult(A, s);
2188 t2 = bimSub(b, t1);
2189 t2->skaldiv(p);
2190 delete b;
2191 delete t1;
2192 b = t2;
2193 s->skalmult(pp, R);
2194 t1 = bimAdd(x, s);
2195 delete s;
2196 x->swapMatrix(t1);
2197 delete t1;
2198
2199 if(kern && i==1)
2200 {
2201 bigintmat * ker = bimChangeCoeff(kp, R);
2202 t1 = bimMult(A, ker);
2203 t1->skaldiv(p);
2204 t1->skalmult(n_Init(-1, R), R);
2205 b->appendCol(t1);
2206 delete t1;
2207 x->appendCol(ker);
2208 delete ker;
2209 x_p->extendCols(kp->cols());
2210 eps_p->extendCols(kp->cols());
2211 fps_p->extendCols(kp->cols());
2212 }
2213
2214 n_InpMult(pp, p, R);
2215
2216 if (b->isZero())
2217 {
2218 //exact solution found, stop
2219 delete eps_p;
2220 delete fps_p;
2221 delete x_p;
2222 delete Hp;
2223 delete kp;
2224 delete Tp;
2225 delete b;
2226 n_Delete(&pp, R);
2227 n_Delete(&p, R);
2228 nKillChar(Rp);
2229
2230 return n_Init(1, R);
2231 }
2232 else
2233 {
2234 bigintmat *y = new bigintmat(x->rows(), x->cols(), R);
2235 number d = bimFarey(x, pp, y);
2236 if (d)
2237 {
2238 bigintmat *c = bimMult(A, y);
2239 bigintmat *bd = new bigintmat(B);
2240 bd->skalmult(d, R);
2241 if (kern)
2242 {
2243 bd->extendCols(kp->cols());
2244 }
2245 if (*c == *bd)
2246 {
2247 x->swapMatrix(y);
2248 delete y;
2249 delete c;
2250 if (kern)
2251 {
2252 y = new bigintmat(x->rows(), B->cols(), R);
2253 c = new bigintmat(x->rows(), kp->cols(), R);
2254 x->splitcol(y, c);
2255 x->swapMatrix(y);
2256 delete y;
2257 kern->swapMatrix(c);
2258 delete c;
2259 }
2260
2261 delete bd;
2262
2263 delete eps_p;
2264 delete fps_p;
2265 delete x_p;
2266 delete Hp;
2267 delete kp;
2268 delete Tp;
2269 delete b;
2270 n_Delete(&pp, R);
2271 n_Delete(&p, R);
2272 nKillChar(Rp);
2273
2274 return d;
2275 }
2276 delete c;
2277 delete bd;
2278 n_Delete(&d, R);
2279 }
2280 delete y;
2281 }
2282 i++;
2283 } while (1);
2284 delete eps_p;
2285 delete fps_p;
2286 delete x_p;
2287 delete Hp;
2288 delete kp;
2289 delete Tp;
2290 n_Delete(&pp, R);
2291 n_Delete(&p, R);
2292 nKillChar(Rp);
2293 return NULL;
2294}
2295#endif
2296
2297//TODO: re-write using reduce_mod_howell
2299{
2300 // try to solve Ax=b, more precisely, find
2301 // number d
2302 // bigintmat x
2303 // sth. Ax=db
2304 // where d is small-ish (divides the determinant of A if this makes sense)
2305 // return 0 if there is no solution.
2306 //
2307 // if kern is non-NULL, return a basis for the kernel
2308
2309 //Algo: we do row-howell (triangular matrix). The idea is
2310 // Ax = b <=> AT T^-1x = b
2311 // y := T^-1 x, solve AT y = b
2312 // and return Ty.
2313 //Howell does not compute the trafo, hence we need to cheat:
2314 //B := (I_n | A^t)^t, then the top part of the Howell form of
2315 //B will give a useful trafo
2316 //Then we can find x by back-substitution and lcm/gcd to find the denominator
2317 //The defining property of Howell makes this work.
2318
2319 coeffs R = A->basecoeffs();
2321 m->howell(); // since m contains the identity, we'll have A->cols()
2322 // many cols.
2323 number den = n_Init(1, R);
2324
2325 bigintmat * B = new bigintmat(A->rows(), 1, R);
2326 for(int i=1; i<= b->cols(); i++)
2327 {
2328 int A_col = A->cols();
2329 b->getcol(i, B);
2330 B->skalmult(den, R);
2331 for(int j = B->rows(); j>0; j--)
2332 {
2333 number Ai = m->view(m->rows()-B->rows() + j, A_col);
2334 if (n_IsZero(Ai, R) &&
2335 n_IsZero(B->view(j, 1), R))
2336 {
2337 continue; //all is fine: 0*x = 0
2338 }
2339 else if (n_IsZero(B->view(j, 1), R))
2340 {
2341 x->rawset(x->rows() - B->rows() + j, i, n_Init(0, R));
2342 A_col--;
2343 }
2344 else if (n_IsZero(Ai, R))
2345 {
2346 delete m;
2347 delete B;
2348 n_Delete(&den, R);
2349 return 0;
2350 }
2351 else
2352 {
2353 // solve ax=db, possibly enlarging d
2354 // so x = db/a
2355 number Bj = B->view(j, 1);
2356 number g = n_Gcd(Bj, Ai, R);
2357 number xi;
2358 if (n_Equal(Ai, g, R))
2359 { //good: den stable!
2360 xi = n_Div(Bj, Ai, R);
2361 }
2362 else
2363 { //den <- den * (a/g), so old sol. needs to be adjusted
2364 number inc_d = n_Div(Ai, g, R);
2365 n_InpMult(den, inc_d, R);
2366 x->skalmult(inc_d, R);
2367 B->skalmult(inc_d, R);
2368 xi = n_Div(Bj, g, R);
2369 n_Delete(&inc_d, R);
2370 } //now for the back-substitution:
2371 x->rawset(x->rows() - B->rows() + j, i, xi);
2372 for(int k=j; k>0; k--)
2373 {
2374 //B[k] = B[k] - x[k]A[k][j]
2375 number s = n_Mult(xi, m->view(m->rows()-B->rows() + k, A_col), R);
2376 B->rawset(k, 1, n_Sub(B->view(k, 1), s, R));
2377 n_Delete(&s, R);
2378 }
2379 n_Delete(&g, R);
2380 A_col--;
2381 }
2382 if (!A_col)
2383 {
2384 if (B->isZero()) break;
2385 else
2386 {
2387 delete m;
2388 delete B;
2389 n_Delete(&den, R);
2390 return 0;
2391 }
2392 }
2393 }
2394 }
2395 delete B;
2396 bigintmat *T = new bigintmat(A->cols(), A->cols(), R);
2397 T->copySubmatInto(m, 1, 1, A->cols(), A->cols(), 1, 1);
2398 if (kern)
2399 {
2400 int i, j;
2401 for(i=1; i<= A->cols(); i++)
2402 {
2403 for(j=m->rows(); j>A->cols(); j--)
2404 {
2405 if (!n_IsZero(m->view(j, i), R)) break;
2406 }
2407 if (j>A->cols()) break;
2408 }
2409 Print("Found nullity (kern dim) of %d\n", i-1);
2410 bigintmat * ker = new bigintmat(A->rows(), i-1, R);
2411 ker->copySubmatInto(T, 1, 1, A->rows(), i-1, 1, 1);
2412 kern->swapMatrix(ker);
2413 delete ker;
2414 }
2415 delete m;
2416 bigintmat * y = bimMult(T, x);
2417 x->swapMatrix(y);
2418 delete y;
2419 x->simplifyContentDen(&den);
2420#if 0
2421 PrintS("sol = 1/");
2422 n_Print(den, R);
2423 PrintS(" *\n");
2424 x->Print();
2425 PrintLn();
2426#endif
2427 return den;
2428}
2429
2431{
2432#if 0
2433 PrintS("Solve Ax=b for A=\n");
2434 A->Print();
2435 PrintS("\nb = \n");
2436 b->Print();
2437 PrintS("\nx = \n");
2438 x->Print();
2439 PrintLn();
2440#endif
2441
2442 coeffs R = A->basecoeffs();
2443 assume (R == b->basecoeffs());
2444 assume (R == x->basecoeffs());
2445 assume ((x->cols() == b->cols()) && (x->rows() == A->cols()) && (A->rows() == b->rows()));
2446
2447 switch (getCoeffType(R))
2448 {
2449 #ifdef HAVE_RINGS
2450 case n_Z:
2451 return solveAx_dixon(A, b, x, NULL);
2452 case n_Zn:
2453 case n_Znm:
2454 case n_Z2m:
2455 return solveAx_howell(A, b, x, NULL);
2456 #endif
2457 case n_Zp:
2458 case n_Q:
2459 case n_GF:
2460 case n_algExt:
2461 case n_transExt:
2462 WarnS("have field, should use Gauss or better");
2463 break;
2464 default:
2465 if (R->cfXExtGcd && R->cfAnn)
2466 { //assume it's Euclidean
2467 return solveAx_howell(A, b, x, NULL);
2468 }
2469 WerrorS("have no solve algorithm");
2470 break;
2471 }
2472 return NULL;
2473}
2474
2476{
2477 bigintmat * t, *s, *a=A;
2478 coeffs R = a->basecoeffs();
2479 if (T)
2480 {
2481 *T = new bigintmat(a->cols(), a->cols(), R),
2482 (*T)->one();
2483 t = new bigintmat(*T);
2484 }
2485 else
2486 {
2487 t = *T;
2488 }
2489
2490 if (S)
2491 {
2492 *S = new bigintmat(a->rows(), a->rows(), R);
2493 (*S)->one();
2494 s = new bigintmat(*S);
2495 }
2496 else
2497 {
2498 s = *S;
2499 }
2500
2501 int flip=0;
2502 do
2503 {
2504 bigintmat * x, *X;
2505 if (flip)
2506 {
2507 x = s;
2508 X = *S;
2509 }
2510 else
2511 {
2512 x = t;
2513 X = *T;
2514 }
2515
2516 if (x)
2517 {
2518 x->one();
2519 bigintmat * r = new bigintmat(a->rows()+a->cols(), a->cols(), R);
2520 bigintmat * rw = new bigintmat(1, a->cols(), R);
2521 for(int i=0; i<a->cols(); i++)
2522 {
2523 x->getrow(i+1, rw);
2524 r->setrow(i+1, rw);
2525 }
2526 for (int i=0; i<a->rows(); i++)
2527 {
2528 a->getrow(i+1, rw);
2529 r->setrow(i+a->cols()+1, rw);
2530 }
2531 r->hnf();
2532 for(int i=0; i<a->cols(); i++)
2533 {
2534 r->getrow(i+1, rw);
2535 x->setrow(i+1, rw);
2536 }
2537 for(int i=0; i<a->rows(); i++)
2538 {
2539 r->getrow(i+a->cols()+1, rw);
2540 a->setrow(i+1, rw);
2541 }
2542 delete rw;
2543 delete r;
2544
2545#if 0
2546 Print("X: %ld\n", X);
2547 X->Print();
2548 Print("\nx: %ld\n", x);
2549 x->Print();
2550#endif
2551 bimMult(X, x, X);
2552#if 0
2553 Print("\n2:X: %ld %ld %ld\n", X, *S, *T);
2554 X->Print();
2555 Print("\n2:x: %ld\n", x);
2556 x->Print();
2557 PrintLn();
2558#endif
2559 }
2560 else
2561 {
2562 a->hnf();
2563 }
2564
2565 int diag = 1;
2566 for(int i=a->rows(); diag && i>0; i--)
2567 {
2568 for(int j=a->cols(); j>0; j--)
2569 {
2570 if ((a->rows()-i)!=(a->cols()-j) && !n_IsZero(a->view(i, j), R))
2571 {
2572 diag = 0;
2573 break;
2574 }
2575 }
2576 }
2577#if 0
2578 PrintS("Diag ? %d\n", diag);
2579 a->Print();
2580 PrintLn();
2581#endif
2582 if (diag) break;
2583
2584 a = a->transpose(); // leaks - I need to write inpTranspose
2585 flip = 1-flip;
2586 } while (1);
2587 if (flip)
2588 a = a->transpose();
2589
2590 if (S) *S = (*S)->transpose();
2591 if (s) delete s;
2592 if (t) delete t;
2593 A->copy(a);
2594}
2595
2596#ifdef HAVE_RINGS
2597//a "q-base" for the kernel of a.
2598//Should be re-done to use Howell rather than smith.
2599//can be done using solveAx now.
2600int kernbase (bigintmat *a, bigintmat *c, number p, coeffs q)
2601{
2602#if 0
2603 PrintS("Kernel of ");
2604 a->Print();
2605 PrintS(" modulo ");
2606 n_Print(p, q);
2607 PrintLn();
2608#endif
2609
2610 coeffs coe = numbercoeffs(p, q);
2611 bigintmat *m = bimChangeCoeff(a, coe), *U, *V;
2612 diagonalForm(m, &U, &V);
2613#if 0
2614 PrintS("\ndiag form: ");
2615 m->Print();
2616 PrintS("\nU:\n");
2617 U->Print();
2618 PrintS("\nV:\n");
2619 V->Print();
2620 PrintLn();
2621#endif
2622
2623 int rg = 0;
2624#undef MIN
2625#define MIN(a,b) (a < b ? a : b)
2626 for(rg=0; rg<MIN(m->rows(), m->cols()) && !n_IsZero(m->view(m->rows()-rg,m->cols()-rg), coe); rg++);
2627
2628 bigintmat * k = new bigintmat(m->cols(), m->rows(), coe);
2629 for(int i=0; i<rg; i++)
2630 {
2631 number A = n_Ann(m->view(m->rows()-i, m->cols()-i), coe);
2632 k->set(m->cols()-i, i+1, A);
2633 n_Delete(&A, coe);
2634 }
2635 for(int i=rg; i<m->cols(); i++)
2636 {
2637 k->set(m->cols()-i, i+1-rg, n_Init(1, coe));
2638 }
2639 bimMult(V, k, k);
2640 c->copy(bimChangeCoeff(k, q));
2641 return c->cols();
2642}
2643#endif
2644
2646{
2647 if ((r == NULL) || (s == NULL))
2648 return false;
2649 if (r == s)
2650 return true;
2651 if ((getCoeffType(r)==n_Z) && (getCoeffType(s)==n_Z))
2652 return true;
2653 if ((getCoeffType(r)==n_Zp) && (getCoeffType(s)==n_Zp))
2654 {
2655 if (r->ch == s->ch)
2656 return true;
2657 else
2658 return false;
2659 }
2660 // n_Zn stimmt wahrscheinlich noch nicht
2661 if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
2662 {
2663 if (r->ch == s->ch)
2664 return true;
2665 else
2666 return false;
2667 }
2668 if ((getCoeffType(r)==n_Q) && (getCoeffType(s)==n_Q))
2669 return true;
2670 // FALL n_Zn FEHLT NOCH!
2671 //if ((getCoeffType(r)==n_Zn) && (getCoeffType(s)==n_Zn))
2672 return false;
2673}
2674
2676{
2677 coeffs r = basecoeffs();
2678 number g = get(1,1), h;
2679 int n=rows()*cols();
2680 for(int i=1; i<n && !n_IsOne(g, r); i++)
2681 {
2682 h = n_Gcd(g, view(i), r);
2683 n_Delete(&g, r);
2684 g=h;
2685 }
2686 return g;
2687}
2689{
2690 coeffs r = basecoeffs();
2691 number g = n_Copy(*d, r), h;
2692 int n=rows()*cols();
2693 for(int i=0; i<n && !n_IsOne(g, r); i++)
2694 {
2695 h = n_Gcd(g, view(i), r);
2696 n_Delete(&g, r);
2697 g=h;
2698 }
2699 *d = n_Div(*d, g, r);
2700 if (!n_IsOne(g, r))
2701 skaldiv(g);
2702}
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)
Definition: bigintmat.cc:1950
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.
Definition: bigintmat.cc:2600
#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...
Definition: bigintmat.cc:2430
static number solveAx_dixon(bigintmat *A, bigintmat *B, bigintmat *x, bigintmat *kern)
Definition: bigintmat.cc:2108
bigintmat * bimChangeCoeff(bigintmat *a, coeffs cnew)
Liefert Kopier von Matrix a zurück, mit coeffs cnew statt den ursprünglichen.
Definition: bigintmat.cc:1804
#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)
Definition: bigintmat.cc:2036
bool nCoeffs_are_equal(coeffs r, coeffs s)
Definition: bigintmat.cc:2645
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:529
bigintmat * bimSub(bigintmat *a, bigintmat *b)
Definition: bigintmat.cc:218
static number solveAx_howell(bigintmat *A, bigintmat *b, bigintmat *x, bigintmat *kern)
Definition: bigintmat.cc:2298
void diagonalForm(bigintmat *A, bigintmat **S, bigintmat **T)
Definition: bigintmat.cc:2475
bigintmat * iv2bim(intvec *b, const coeffs C)
Definition: bigintmat.cc:349
static int findLongest(int *a, int length)
Definition: bigintmat.cc:537
static int getShorter(int *a, int l, int j, int cols, int rows)
Definition: bigintmat.cc:552
static number bimFarey(bigintmat *A, number N, bigintmat *L)
Definition: bigintmat.cc:2048
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.
Definition: bigintmat.cc:1804
bool nCoeffs_are_equal(coeffs r, coeffs s)
Definition: bigintmat.cc:2645
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:4082
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
CanonicalForm cf
Definition: cfModGcd.cc:4083
CanonicalForm b
Definition: cfModGcd.cc:4103
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:443
void colskaldiv(int j, number b)
Macht Ganzzahldivision aller j-ten Spalteneinträge mit b.
Definition: bigintmat.cc:1876
number det()
det (via LaPlace in general, hnf for euc. rings)
Definition: bigintmat.cc:1512
void splitrow(bigintmat *a, bigintmat *b)
Speichert in Matrix a den oberen, in b den unteren Teil der Matrix, vorausgesetzt die Dimensionen sti...
Definition: bigintmat.cc:1127
int isOne()
is matrix is identity
Definition: bigintmat.cc:1300
void zero()
Setzt alle Einträge auf 0.
Definition: bigintmat.cc:1350
void appendCol(bigintmat *a)
horizontally join the matrices, m <- m|a
Definition: bigintmat.cc:1083
number trace()
the trace ....
Definition: bigintmat.cc:1498
bool addrow(int i, int j, number a, coeffs c)
... Zeile ...
Definition: bigintmat.cc:983
void swap(int i, int j)
swap columns i and j
Definition: bigintmat.cc:685
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)
Definition: bigintmat.cc:1889
bool addcol(int i, int j, number a, coeffs c)
addiert a-faches der j-ten Spalte zur i-ten dazu
Definition: bigintmat.cc:959
number * v
Definition: bigintmat.h:54
void hnf()
transforms INPLACE to HNF
Definition: bigintmat.cc:1660
char * StringAsPrinted()
Returns a string as it would have been printed in the interpreter.
Definition: bigintmat.cc:451
void swapMatrix(bigintmat *a)
Definition: bigintmat.cc:1566
int cols() const
Definition: bigintmat.h:144
int isZero()
Definition: bigintmat.cc:1363
number hnfdet()
det via HNF Primzahlen als long long int, müssen noch in number umgewandelt werden?
Definition: bigintmat.cc:1545
int findnonzero(int i)
find index of 1st non-zero entry in row i
Definition: bigintmat.cc:723
bigintmat * modhnf(number p, coeffs c)
computes HNF(this | p*I)
Definition: bigintmat.cc:1832
void setcol(int j, bigintmat *m)
Setzt j-te Spalte gleich übergebenem Vektor (Matrix) m.
Definition: bigintmat.cc:826
bool add(bigintmat *b)
Addiert zur Matrix die Matrix b dazu. Return false => an error occurred.
Definition: bigintmat.cc:894
int * getwid(int maxwid)
Definition: bigintmat.cc:580
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:860
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 ...
Definition: bigintmat.cc:1023
bool skalmult(number b, coeffs c)
Multipliziert zur Matrix den Skalar b hinzu.
Definition: bigintmat.cc:938
void simplifyContentDen(number *den)
ensures that Gcd(den, content)=1 enden hier wieder
Definition: bigintmat.cc:2688
void extendCols(int i)
append i zero-columns to the matrix
Definition: bigintmat.cc:1076
void splitcol(bigintmat *a, bigintmat *b)
... linken ... rechten ...
Definition: bigintmat.cc:1169
number content()
the content, the gcd of all entries. Only makes sense for Euclidean rings (or possibly constructive P...
Definition: bigintmat.cc:2675
void skaldiv(number b)
Macht Ganzzahldivision aller Matrixeinträge mit b.
Definition: bigintmat.cc:1861
void colskalmult(int i, number a, coeffs c)
Multipliziert zur i-ten Spalte den Skalar a hinzu.
Definition: bigintmat.cc:1007
int colIsZero(int i)
Definition: bigintmat.cc:1577
bool copy(bigintmat *b)
Kopiert Einträge von b auf Bigintmat.
Definition: bigintmat.cc:1259
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:747
number pseudoinv(bigintmat *a)
Speichert in Matrix a die Pseudoinverse, liefert den Nenner zurück.
Definition: bigintmat.cc:1415
int col
Definition: bigintmat.h:56
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:610
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:778
void concatcol(bigintmat *a, bigintmat *b)
Definition: bigintmat.cc:1098
int findcolnonzero(int j)
find index of 1st non-zero entry in column j
Definition: bigintmat.cc:735
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....
Definition: bigintmat.cc:1287
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.
Definition: bigintmat.cc:1325
void howell()
dito, but Howell form (only different for zero-divsors)
Definition: bigintmat.cc:1585
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:791
void concatrow(bigintmat *a, bigintmat *b)
Fügt zwei Matrixen untereinander/nebeneinander in gegebene Matrix ein, bzw spaltet gegebenen Matrix a...
Definition: bigintmat.cc:1039
bigintmat()
Definition: bigintmat.h:59
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.
Definition: bigintmat.cc:1381
int compare(const bigintmat *op) const
Definition: bigintmat.cc:362
bool sub(bigintmat *b)
Subtrahiert ...
Definition: bigintmat.cc:916
int row
Definition: bigintmat.h:55
char * String()
IO: String returns a singular string containing the matrix, needs freeing afterwards.
Definition: bigintmat.cc:436
void swaprow(int i, int j)
swap rows i and j
Definition: bigintmat.cc:704
void mod(number p)
Reduziert komplette Matrix modulo p.
Definition: bigintmat.cc:1916
Definition: intvec.h:23
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:633
static FORCE_INLINE void number2mpz(number n, coeffs c, mpz_t m)
Definition: coeffs.h:984
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:544
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:448
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:647
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:600
@ 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:661
static FORCE_INLINE number n_QuotRem(number a, number b, number *q, const coeffs r)
Definition: coeffs.h:678
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:956
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:491
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:676
void n_Print(number &a, const coeffs r)
print a number (BEWARE of string buffers!) mostly for debugging
Definition: numbers.cc:667
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:697
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:554
static FORCE_INLINE number n_Farey(number a, number b, const coeffs r)
Definition: coeffs.h:764
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:612
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:413
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:508
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff 'n' represents the zero element.
Definition: coeffs.h:461
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:529
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:652
static FORCE_INLINE n_coeffType getCoeffType(const coeffs r)
Returns the type of coeffs domain.
Definition: coeffs.h:422
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:452
static FORCE_INLINE void n_Write(number n, const coeffs r, const BOOLEAN bShortOut=TRUE)
Definition: coeffs.h:588
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:535
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:625
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:638
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:457
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:605
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
void nKillChar(coeffs r)
undo all initialisations
Definition: numbers.cc:568
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:643
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff 'n' represents the one element.
Definition: coeffs.h:465
static FORCE_INLINE number n_XExtGcd(number a, number b, number *s, number *t, number *u, number *v, const coeffs r)
Definition: coeffs.h:670
#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:72
CFList tmp2
Definition: facFqBivar.cc:72
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:389
The main handler for Singular numbers which are suitable for Singular polynomials.
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omFree(addr)
Definition: omAllocDecl.h:261
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define NULL
Definition: omList.c:12
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
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)
Definition: walkSupport.cc:836