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