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