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