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