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