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