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