linalg_from_matpol.cc
Go to the documentation of this file.
1 typedef int perm[100];
2 static void mpReplace(int j, int n, int &sign, int *perm);
3 static int mpNextperm(perm * z, int max);
4 static poly mp_Leibnitz(matrix a, const ring);
5 static poly minuscopy (poly p, const ring);
6 static poly p_Insert(poly p1, poly p2, const ring);
7 
8 static void mp_PartClean(matrix, int, int, const ring);
9 static int mp_PrepareRow (matrix, int, int, const ring);
10 static int mp_PreparePiv (matrix, int, int, const ring);
11 static int mp_PivBar(matrix, int, int, const ring);
12 static int mp_PivRow(matrix, int, int, const ring);
13 static float mp_PolyWeight(poly, const ring);
14 static void mp_SwapRow(matrix, int, int, int, const ring);
15 static void mp_SwapCol(matrix, int, int, int, const ring);
16 static void mp_ElimBar(matrix, matrix, poly, int, int, const ring);
17 
18 /*2
19 * prepare one step of 'Bareiss' algorithm
20 * for application in minor
21 */
22 static int mp_PrepareRow (matrix a, int lr, int lc, const ring R)
23 {
24  int r;
25 
26  r = mp_PivBar(a,lr,lc, R);
27  if(r==0) return 0;
28  if(r<lr) mp_SwapRow(a, r, lr, lc, R);
29  return 1;
30 }
31 
32 /*2
33 * prepare one step of 'Bareiss' algorithm
34 * for application in minor
35 */
36 static int mp_PreparePiv (matrix a, int lr, int lc, const ring R)
37 {
38  int c;
39 
40  c = mp_PivRow(a, lr, lc, R);
41  if(c==0) return 0;
42  if(c<lc) mp_SwapCol(a, c, lr, lc, R);
43  return 1;
44 }
45 
46 /*
47 * find best row
48 */
49 static int mp_PivBar(matrix a, int lr, int lc, const ring R)
50 {
51  float f1, f2;
52  poly *q1;
53  int i,j,io;
54 
55  io = -1;
56  f1 = 1.0e30;
57  for (i=lr-1;i>=0;i--)
58  {
59  q1 = &(a->m)[i*a->ncols];
60  f2 = 0.0;
61  for (j=lc-1;j>=0;j--)
62  {
63  if (q1[j]!=NULL)
64  f2 += mp_PolyWeight(q1[j], R);
65  }
66  if ((f2!=0.0) && (f2<f1))
67  {
68  f1 = f2;
69  io = i;
70  }
71  }
72  if (io<0) return 0;
73  else return io+1;
74 }
75 
76 /*
77 * find pivot in the last row
78 */
79 static int mp_PivRow(matrix a, int lr, int lc, const ring R)
80 {
81  float f1, f2;
82  poly *q1;
83  int j,jo;
84 
85  jo = -1;
86  f1 = 1.0e30;
87  q1 = &(a->m)[(lr-1)*a->ncols];
88  for (j=lc-1;j>=0;j--)
89  {
90  if (q1[j]!=NULL)
91  {
92  f2 = mp_PolyWeight(q1[j], R);
93  if (f2<f1)
94  {
95  f1 = f2;
96  jo = j;
97  }
98  }
99  }
100  if (jo<0) return 0;
101  else return jo+1;
102 }
103 
104 /*
105 * weigth of a polynomial, for pivot strategy
106 */
107 static float mp_PolyWeight(poly p, const ring R)
108 {
109  int i;
110  float res;
111 
112  if (pNext(p) == NULL)
113  {
114  res = (float)n_Size(p_GetCoeff(p, R), R);
115  for (i=rVar(R);i>0;i--)
116  {
117  if(p_GetExp(p,i, R)!=0)
118  {
119  res += 2.0;
120  break;
121  }
122  }
123  }
124  else
125  {
126  res = 0.0;
127  do
128  {
129  res += (float)n_Size(p_GetCoeff(p, R), R) + 2.0;
130  pIter(p);
131  }
132  while (p);
133  }
134  return res;
135 }
136 
137 static void mpSwapRow(matrix a, int pos, int lr, int lc)
138 {
139  poly sw;
140  int j;
141  poly* a2 = a->m;
142  poly* a1 = &a2[a->ncols*(pos-1)];
143 
144  a2 = &a2[a->ncols*(lr-1)];
145  for (j=lc-1; j>=0; j--)
146  {
147  sw = a1[j];
148  a1[j] = a2[j];
149  a2[j] = sw;
150  }
151 }
152 
153 static void mpSwapCol(matrix a, int pos, int lr, int lc)
154 {
155  poly sw;
156  int j;
157  poly* a2 = a->m;
158  poly* a1 = &a2[pos-1];
159 
160  a2 = &a2[lc-1];
161  for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
162  {
163  sw = a1[j];
164  a1[j] = a2[j];
165  a2[j] = sw;
166  }
167 }
168 
169 /*
170 * C++ classes for Bareiss algorithm
171 */
173 {
174  private:
175  int ym, yn;
176  public:
177  float *wrow, *wcol;
178  row_col_weight() : ym(0) {}
179  row_col_weight(int, int);
180  ~row_col_weight();
181 };
182 
183 /*2
184 * a submatrix M of a matrix X[m,n]:
185 * 0 <= i < s_m <= a_m
186 * 0 <= j < s_n <= a_n
187 * M = ( Xarray[qrow[i],qcol[j]] )
188 * if a_m = a_n and s_m = s_n
189 * det(X) = sign*div^(s_m-1)*det(M)
190 * resticted pivot for elimination
191 * 0 <= j < piv_s
192 */
194 {
195  private:
196  int a_m, a_n, s_m, s_n, sign, piv_s;
197  int *qrow, *qcol;
199  ring R;
200 
201  void mpInitMat();
202  poly * mpRowAdr(int);
203  poly * mpColAdr(int);
204  void mpRowWeight(float *);
205  void mpColWeight(float *);
206  void mpRowSwap(int, int);
207  void mpColSwap(int, int);
208  public:
209  mp_permmatrix() : a_m(0), R(NULL) {}
210  mp_permmatrix(matrix, const ring);
212  ~mp_permmatrix();
213  int mpGetRow();
214  int mpGetCol();
215  int mpGetRdim();
216  int mpGetCdim();
217  int mpGetSign();
218  void mpSetSearch(int s);
219  void mpSaveArray();
220  poly mpGetElem(int, int);
221  void mpSetElem(poly, int, int);
222  void mpDelElem(int, int);
223  void mpElimBareiss(poly);
224  int mpPivotBareiss(row_col_weight *);
225  int mpPivotRow(row_col_weight *, int);
226  void mpToIntvec(intvec *);
227  void mpRowReorder();
228  void mpColReorder();
229 };
230 
231 
232 static poly minuscopy (poly p, const ring R)
233 {
234  poly w;
235  number e;
236  e = n_Init(-1, R);
237  w = p_Copy(p, R);
238  p_Mult_nn(w, e, R);
239  n_Delete(&e, R);
240  return w;
241 }
242 
243 
244 /*2
245 *returns the determinant of the matrix m;
246 *uses Newtons formulea for symmetric functions
247 */
248 poly mp_Det (matrix m, const ring R)
249 {
250  int i,j,k,n;
251  poly p,q;
252  matrix a, s;
253  matrix ma[100];
254  number c=NULL, d=NULL, ONE=NULL;
255 
256  n = MATROWS(m);
257  if (n != MATCOLS(m))
258  {
259  Werror("det of %d x %d matrix",n,MATCOLS(m));
260  return NULL;
261  }
262  k=rChar(R);
263  if ((k > 0) && (k <= n))
264  return mp_Leibnitz(m, R);
265  ONE = n_Init(1, R);
266  ma[1]=mp_Copy(m, R);
267  k = (n+1) / 2;
268  s = mpNew(1, n);
269  MATELEM(s,1,1) = mp_Trace(m, R);
270  for (i=2; i<=k; i++)
271  {
272  //ma[i] = mpNew(n,n);
273  ma[i]=mp_Mult(ma[i-1], ma[1], R);
274  MATELEM(s,1,i) = mp_Trace(ma[i], R);
275  p_Test(MATELEM(s,1,i), R);
276  }
277  for (i=k+1; i<=n; i++)
278  {
279  MATELEM(s,1,i) = TraceOfProd(ma[i / 2], ma[(i+1) / 2], n, R);
280  p_Test(MATELEM(s,1,i), R);
281  }
282  for (i=1; i<=k; i++)
283  id_Delete((ideal *)&(ma[i]), R);
284 /* the array s contains the traces of the powers of the matrix m,
285 * these are the power sums of the eigenvalues of m */
286  a = mpNew(1,n);
287  MATELEM(a,1,1) = minuscopy(MATELEM(s,1,1), R);
288  for (i=2; i<=n; i++)
289  {
290  p = p_Copy(MATELEM(s,1,i), R);
291  for (j=i-1; j>=1; j--)
292  {
293  q = pp_Mult_qq(MATELEM(s,1,j), MATELEM(a,1,i-j), R);
294  p_Test(q, R);
295  p = p_Add_q(p,q, R);
296  }
297  // c= -1/i
298  d = n_Init(-(int)i, R);
299  c = n_Div(ONE, d, R);
300  n_Delete(&d, R);
301 
302  p_Mult_nn(p, c, R);
303  p_Test(p, R);
304  MATELEM(a,1,i) = p;
305  n_Delete(&c, R);
306  }
307 /* the array a contains the elementary symmetric functions of the
308 * eigenvalues of m */
309  for (i=1; i<=n-1; i++)
310  {
311  //p_Delete(&(MATELEM(a,1,i)), R);
312  p_Delete(&(MATELEM(s,1,i)), R);
313  }
314  p_Delete(&(MATELEM(s,1,n)), R);
315 /* up to a sign, the determinant is the n-th elementary symmetric function */
316  if ((n/2)*2 < n)
317  {
318  d = n_Init(-1, R);
319  p_Mult_nn(MATELEM(a,1,n), d, R);
320  n_Delete(&d, R);
321  }
322  n_Delete(&ONE, R);
323  id_Delete((ideal *)&s, R);
324  poly result=MATELEM(a,1,n);
325  MATELEM(a,1,n)=NULL;
326  id_Delete((ideal *)&a, R);
327  return result;
328 }
329 
330 
331 ///*2
332 //*homogenize all elements of matrix (not the matrix itself)
333 //*/
334 //matrix mpHomogen(matrix a, int v)
335 //{
336 // int i,j;
337 // poly p;
338 //
339 // for (i=1;i<=MATROWS(a);i++)
340 // {
341 // for (j=1;j<=MATCOLS(a);j++)
342 // {
343 // p=pHomogen(MATELEM(a,i,j),v);
344 // p_Delete(&(MATELEM(a,i,j)), ?);
345 // MATELEM(a,i,j)=p;
346 // }
347 // }
348 // return a;
349 //}
350 
351 
352 
353 /* --------------- internal stuff ------------------- */
354 
356 {
357  ym = i;
358  yn = j;
359  wrow = (float *)omAlloc(i*sizeof(float));
360  wcol = (float *)omAlloc(j*sizeof(float));
361 }
362 
364 {
365  if (ym!=0)
366  {
367  omFreeSize((ADDRESS)wcol, yn*sizeof(float));
368  omFreeSize((ADDRESS)wrow, ym*sizeof(float));
369  }
370 }
371 
373 {
374  a_m = A->nrows;
375  a_n = A->ncols;
376  this->mpInitMat();
377  Xarray = A->m;
378 }
379 
381 {
382  poly p, *athis, *aM;
383  int i, j;
384 
385  a_m = M->s_m;
386  a_n = M->s_n;
387  sign = M->sign;
388  R = M->R;
389 
390  this->mpInitMat();
391  Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
392  for (i=a_m-1; i>=0; i--)
393  {
394  athis = this->mpRowAdr(i);
395  aM = M->mpRowAdr(i);
396  for (j=a_n-1; j>=0; j--)
397  {
398  p = aM[M->qcol[j]];
399  if (p)
400  {
401  athis[j] = p_Copy(p, R);
402  }
403  }
404  }
405 }
406 
408 {
409  int k;
410 
411  if (a_m != 0)
412  {
413  omFreeSize((ADDRESS)qrow,a_m*sizeof(int));
414  omFreeSize((ADDRESS)qcol,a_n*sizeof(int));
415  if (Xarray != NULL)
416  {
417  for (k=a_m*a_n-1; k>=0; k--)
418  p_Delete(&Xarray[k], R);
419  omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
420  }
421  }
422 }
423 
424 int mp_permmatrix::mpGetRdim() { return s_m; }
425 
426 int mp_permmatrix::mpGetCdim() { return s_n; }
427 
429 
431 
433 
435 {
436  return Xarray[a_n*qrow[r]+qcol[c]];
437 }
438 
439 void mp_permmatrix::mpSetElem(poly p, int r, int c)
440 {
441  Xarray[a_n*qrow[r]+qcol[c]] = p;
442 }
443 
444 void mp_permmatrix::mpDelElem(int r, int c)
445 {
446  p_Delete(&Xarray[a_n*qrow[r]+qcol[c]], R);
447 }
448 
449 /*
450 * the Bareiss-type elimination with division by div (div != NULL)
451 */
453 {
454  poly piv, elim, q1, q2, *ap, *a;
455  int i, j, jj;
456 
457  ap = this->mpRowAdr(s_m);
458  piv = ap[qcol[s_n]];
459  for(i=s_m-1; i>=0; i--)
460  {
461  a = this->mpRowAdr(i);
462  elim = a[qcol[s_n]];
463  if (elim != NULL)
464  {
465  elim = p_Neg(elim, R);
466  for (j=s_n-1; j>=0; j--)
467  {
468  q2 = NULL;
469  jj = qcol[j];
470  if (ap[jj] != NULL)
471  {
472  q2 = SM_MULT(ap[jj], elim, div, R);
473  if (a[jj] != NULL)
474  {
475  q1 = SM_MULT(a[jj], piv, div, R);
476  p_Delete(&a[jj], R);
477  q2 = p_Add_q(q2, q1, R);
478  }
479  }
480  else if (a[jj] != NULL)
481  {
482  q2 = SM_MULT(a[jj], piv, div, R);
483  }
484  if ((q2!=NULL) && div)
485  SM_DIV(q2, div, R);
486  a[jj] = q2;
487  }
488  p_Delete(&a[qcol[s_n]], R);
489  }
490  else
491  {
492  for (j=s_n-1; j>=0; j--)
493  {
494  jj = qcol[j];
495  if (a[jj] != NULL)
496  {
497  q2 = SM_MULT(a[jj], piv, div, R);
498  p_Delete(&a[jj], R);
499  if (div)
500  SM_DIV(q2, div, R);
501  a[jj] = q2;
502  }
503  }
504  }
505  }
506 }
507 
508 /*2
509 * pivot strategy for Bareiss algorithm
510 */
512 {
513  poly p, *a;
514  int i, j, iopt, jopt;
515  float sum, f1, f2, fo, r, ro, lp;
516  float *dr = C->wrow, *dc = C->wcol;
517 
518  fo = 1.0e20;
519  ro = 0.0;
520  iopt = jopt = -1;
521 
522  s_n--;
523  s_m--;
524  if (s_m == 0)
525  return 0;
526  if (s_n == 0)
527  {
528  for(i=s_m; i>=0; i--)
529  {
530  p = this->mpRowAdr(i)[qcol[0]];
531  if (p)
532  {
533  f1 = mp_PolyWeight(p, R);
534  if (f1 < fo)
535  {
536  fo = f1;
537  if (iopt >= 0)
538  p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), R);
539  iopt = i;
540  }
541  else
542  p_Delete(&(this->mpRowAdr(i)[qcol[0]]), R);
543  }
544  }
545  if (iopt >= 0)
546  mpReplace(iopt, s_m, sign, qrow);
547  return 0;
548  }
549  this->mpRowWeight(dr);
550  this->mpColWeight(dc);
551  sum = 0.0;
552  for(i=s_m; i>=0; i--)
553  sum += dr[i];
554  for(i=s_m; i>=0; i--)
555  {
556  r = dr[i];
557  a = this->mpRowAdr(i);
558  for(j=s_n; j>=0; j--)
559  {
560  p = a[qcol[j]];
561  if (p)
562  {
563  lp = mp_PolyWeight(p, R);
564  ro = r - lp;
565  f1 = ro * (dc[j]-lp);
566  if (f1 != 0.0)
567  {
568  f2 = lp * (sum - ro - dc[j]);
569  f2 += f1;
570  }
571  else
572  f2 = lp-r-dc[j];
573  if (f2 < fo)
574  {
575  fo = f2;
576  iopt = i;
577  jopt = j;
578  }
579  }
580  }
581  }
582  if (iopt < 0)
583  return 0;
584  mpReplace(iopt, s_m, sign, qrow);
585  mpReplace(jopt, s_n, sign, qcol);
586  return 1;
587 }
588 
589 /*2
590 * pivot strategy for Bareiss algorithm with defined row
591 */
593 {
594  poly p, *a;
595  int j, iopt, jopt;
596  float sum, f1, f2, fo, r, ro, lp;
597  float *dr = C->wrow, *dc = C->wcol;
598 
599  fo = 1.0e20;
600  ro = 0.0;
601  iopt = jopt = -1;
602 
603  s_n--;
604  s_m--;
605  if (s_m == 0)
606  return 0;
607  if (s_n == 0)
608  {
609  p = this->mpRowAdr(row)[qcol[0]];
610  if (p)
611  {
612  f1 = mp_PolyWeight(p, R);
613  if (f1 < fo)
614  {
615  fo = f1;
616  if (iopt >= 0)
617  p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]), R);
618  iopt = row;
619  }
620  else
621  p_Delete(&(this->mpRowAdr(row)[qcol[0]]), R);
622  }
623  if (iopt >= 0)
624  mpReplace(iopt, s_m, sign, qrow);
625  return 0;
626  }
627  this->mpRowWeight(dr);
628  this->mpColWeight(dc);
629  sum = 0.0;
630  for(j=s_m; j>=0; j--)
631  sum += dr[j];
632  r = dr[row];
633  a = this->mpRowAdr(row);
634  for(j=s_n; j>=0; j--)
635  {
636  p = a[qcol[j]];
637  if (p)
638  {
639  lp = mp_PolyWeight(p, R);
640  ro = r - lp;
641  f1 = ro * (dc[j]-lp);
642  if (f1 != 0.0)
643  {
644  f2 = lp * (sum - ro - dc[j]);
645  f2 += f1;
646  }
647  else
648  f2 = lp-r-dc[j];
649  if (f2 < fo)
650  {
651  fo = f2;
652  iopt = row;
653  jopt = j;
654  }
655  }
656  }
657  if (iopt < 0)
658  return 0;
659  mpReplace(iopt, s_m, sign, qrow);
660  mpReplace(jopt, s_n, sign, qcol);
661  return 1;
662 }
663 
665 {
666  int i;
667 
668  for (i=v->rows()-1; i>=0; i--)
669  (*v)[i] = qcol[i]+1;
670 }
671 
673 {
674  int k, i, i1, i2;
675 
676  if (a_m > a_n)
677  k = a_m - a_n;
678  else
679  k = 0;
680  for (i=a_m-1; i>=k; i--)
681  {
682  i1 = qrow[i];
683  if (i1 != i)
684  {
685  this->mpRowSwap(i1, i);
686  i2 = 0;
687  while (qrow[i2] != i) i2++;
688  qrow[i2] = i1;
689  }
690  }
691 }
692 
694 {
695  int k, j, j1, j2;
696 
697  if (a_n > a_m)
698  k = a_n - a_m;
699  else
700  k = 0;
701  for (j=a_n-1; j>=k; j--)
702  {
703  j1 = qcol[j];
704  if (j1 != j)
705  {
706  this->mpColSwap(j1, j);
707  j2 = 0;
708  while (qcol[j2] != j) j2++;
709  qcol[j2] = j1;
710  }
711  }
712 }
713 
714 // private
716 {
717  int k;
718 
719  s_m = a_m;
720  s_n = a_n;
721  piv_s = 0;
722  qrow = (int *)omAlloc(a_m*sizeof(int));
723  qcol = (int *)omAlloc(a_n*sizeof(int));
724  for (k=a_m-1; k>=0; k--) qrow[k] = k;
725  for (k=a_n-1; k>=0; k--) qcol[k] = k;
726 }
727 
729 {
730  return &(Xarray[a_n*qrow[r]]);
731 }
732 
734 {
735  return &(Xarray[qcol[c]]);
736 }
737 
738 void mp_permmatrix::mpRowWeight(float *wrow)
739 {
740  poly p, *a;
741  int i, j;
742  float count;
743 
744  for (i=s_m; i>=0; i--)
745  {
746  a = this->mpRowAdr(i);
747  count = 0.0;
748  for(j=s_n; j>=0; j--)
749  {
750  p = a[qcol[j]];
751  if (p)
752  count += mp_PolyWeight(p, R);
753  }
754  wrow[i] = count;
755  }
756 }
757 
758 void mp_permmatrix::mpColWeight(float *wcol)
759 {
760  poly p, *a;
761  int i, j;
762  float count;
763 
764  for (j=s_n; j>=0; j--)
765  {
766  a = this->mpColAdr(j);
767  count = 0.0;
768  for(i=s_m; i>=0; i--)
769  {
770  p = a[a_n*qrow[i]];
771  if (p)
772  count += mp_PolyWeight(p, R);
773  }
774  wcol[j] = count;
775  }
776 }
777 
778 void mp_permmatrix::mpRowSwap(int i1, int i2)
779 {
780  poly p, *a1, *a2;
781  int j;
782 
783  a1 = &(Xarray[a_n*i1]);
784  a2 = &(Xarray[a_n*i2]);
785  for (j=a_n-1; j>= 0; j--)
786  {
787  p = a1[j];
788  a1[j] = a2[j];
789  a2[j] = p;
790  }
791 }
792 
793 void mp_permmatrix::mpColSwap(int j1, int j2)
794 {
795  poly p, *a1, *a2;
796  int i, k = a_n*a_m;
797 
798  a1 = &(Xarray[j1]);
799  a2 = &(Xarray[j2]);
800  for (i=0; i< k; i+=a_n)
801  {
802  p = a1[i];
803  a1[i] = a2[i];
804  a2[i] = p;
805  }
806 }
807 
809 {
810  return qrow[s_m];
811 }
812 
814 {
815  return qcol[s_n];
816 }
817 
818 /// perform replacement for pivot strategy in Bareiss algorithm
819 /// change sign of determinant
820 static void mpReplace(int j, int n, int &sign, int *perm)
821 {
822  int k;
823 
824  if (j != n)
825  {
826  k = perm[n];
827  perm[n] = perm[j];
828  perm[j] = k;
829  sign = -sign;
830  }
831 }
832 
833 static int mpNextperm(perm * z, int max)
834 {
835  int s, i, k, t;
836  s = max;
837  do
838  {
839  s--;
840  }
841  while ((s > 0) && ((*z)[s] >= (*z)[s+1]));
842  if (s==0)
843  return 0;
844  do
845  {
846  (*z)[s]++;
847  k = 0;
848  do
849  {
850  k++;
851  }
852  while (((*z)[k] != (*z)[s]) && (k!=s));
853  }
854  while (k < s);
855  for (i=s+1; i <= max; i++)
856  {
857  (*z)[i]=0;
858  do
859  {
860  (*z)[i]++;
861  k=0;
862  do
863  {
864  k++;
865  }
866  while (((*z)[k] != (*z)[i]) && (k != i));
867  }
868  while (k < i);
869  }
870  s = max+1;
871  do
872  {
873  s--;
874  }
875  while ((s > 0) && ((*z)[s] > (*z)[s+1]));
876  t = 1;
877  for (i=1; i<max; i++)
878  for (k=i+1; k<=max; k++)
879  if ((*z)[k] < (*z)[i])
880  t = -t;
881  (*z)[0] = t;
882  return s;
883 }
884 
885 static poly mp_Leibnitz(matrix a, const ring R)
886 {
887  int i, e, n;
888  poly p, d;
889  perm z;
890 
891  n = MATROWS(a);
892  memset(&z,0,(n+2)*sizeof(int));
893  p = p_One(R);
894  for (i=1; i <= n; i++)
895  p = p_Mult_q(p, p_Copy(MATELEM(a, i, i), R), R);
896  d = p;
897  for (i=1; i<= n; i++)
898  z[i] = i;
899  z[0]=1;
900  e = 1;
901  if (n!=1)
902  {
903  while (e)
904  {
905  e = mpNextperm((perm *)&z, n);
906  p = p_One(R);
907  for (i = 1; i <= n; i++)
908  p = p_Mult_q(p, p_Copy(MATELEM(a, i, z[i]), R), R);
909  if (z[0] > 0)
910  d = p_Add_q(d, p, R);
911  else
912  d = p_Sub(d, p, R);
913  }
914  }
915  return d;
916 }
917 
918 static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
919 {
920  int r=lr-1, c=lc-1;
921  poly *b = a0->m, *x = re->m;
922  poly piv, elim, q1, q2, *ap, *a, *q;
923  int i, j;
924 
925  ap = &b[r*a0->ncols];
926  piv = ap[c];
927  for(j=c-1; j>=0; j--)
928  if (ap[j] != NULL) ap[j] = p_Neg(ap[j], R);
929  for(i=r-1; i>=0; i--)
930  {
931  a = &b[i*a0->ncols];
932  q = &x[i*re->ncols];
933  if (a[c] != NULL)
934  {
935  elim = a[c];
936  for (j=c-1; j>=0; j--)
937  {
938  q1 = NULL;
939  if (a[j] != NULL)
940  {
941  q1 = SM_MULT(a[j], piv, div, R);
942  if (ap[j] != NULL)
943  {
944  q2 = SM_MULT(ap[j], elim, div, R);
945  q1 = p_Add_q(q1,q2, R);
946  }
947  }
948  else if (ap[j] != NULL)
949  q1 = SM_MULT(ap[j], elim, div, R);
950  if (q1 != NULL)
951  {
952  if (div)
953  SM_DIV(q1, div, R);
954  q[j] = q1;
955  }
956  }
957  }
958  else
959  {
960  for (j=c-1; j>=0; j--)
961  {
962  if (a[j] != NULL)
963  {
964  q1 = SM_MULT(a[j], piv, div, R);
965  if (div)
966  SM_DIV(q1, div, R);
967  q[j] = q1;
968  }
969  }
970  }
971  }
972 }
973 
int status int void size_t count
Definition: si_signals.h:59
static void mp_PartClean(matrix, int, int, const ring)
const CanonicalForm int s
Definition: facAbsFact.cc:55
void mpElimBareiss(poly)
const poly a
Definition: syzextra.cc:212
static int mpNextperm(perm *z, int max)
void mpRowSwap(int, int)
static void mp_SwapCol(matrix, int, int, int, const ring)
int ncols
Definition: matpol.h:22
static void mp_SwapRow(matrix, int, int, int, const ring)
return P p
Definition: myNF.cc:203
poly mp_Det(matrix m, const ring R)
poly mp_Trace(matrix a, const ring R)
Definition: matpol.cc:288
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 rows() const
Definition: intvec.h:88
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
int rChar(ring r)
Definition: ring.cc:686
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:583
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
static int mp_PivRow(matrix, int, int, const ring)
void mpColWeight(float *)
void * ADDRESS
Definition: auxiliary.h:115
int k
Definition: cfEzgcd.cc:93
void mpRowWeight(float *)
static void mpSwapRow(matrix a, int pos, int lr, int lc)
void mpSetElem(poly, int, int)
int mpPivotBareiss(row_col_weight *)
#define omAlloc(size)
Definition: omAllocDecl.h:210
poly p_Sub(poly p1, poly p2, const ring r)
Definition: p_polys.cc:1943
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:804
CanonicalForm lc(const CanonicalForm &f)
#define pIter(p)
Definition: monomials.h:44
poly res
Definition: myNF.cc:322
#define M
Definition: sirandom.c:24
void mpSetSearch(int s)
void mpToIntvec(intvec *)
poly * m
Definition: matpol.h:19
const ring r
Definition: syzextra.cc:208
Definition: intvec.h:14
poly p_One(const ring r)
Definition: p_polys.cc:1312
for(int i=0;i< R->ExpL_Size;i++) Print("%09lx "
Definition: cfEzgcd.cc:66
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent : the integer VarOffset encodes:
Definition: p_polys.h:464
int nrows
Definition: matpol.h:21
int j
Definition: myNF.cc:70
static int max(int a, int b)
Definition: fast_mult.cc:264
static poly p_Insert(poly p1, poly p2, const ring)
static int mp_PrepareRow(matrix, int, int, const ring)
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1070
int mpPivotRow(row_col_weight *, int)
#define A
Definition: sirandom.c:23
const ring R
Definition: DebugPrint.cc:36
static poly minuscopy(poly p, const ring)
poly * mpColAdr(int)
int m
Definition: cfEzgcd.cc:119
int i
Definition: cfEzgcd.cc:123
static poly p_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:895
static int mp_PivBar(matrix, int, int, const ring)
static float mp_PolyWeight(poly, const ring)
#define p_Test(p, r)
Definition: p_polys.h:160
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:47
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:843
static void mpSwapCol(matrix a, int pos, int lr, int lc)
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
void mpDelElem(int, int)
matrix mp_Mult(matrix a, matrix b, const ring R)
Definition: matpol.cc:223
CF_NO_INLINE CanonicalForm div(const CanonicalForm &, const CanonicalForm &)
CF_INLINE CanonicalForm div, mod ( const CanonicalForm & lhs, const CanonicalForm & rhs ) ...
Definition: cf_inline.cc:553
#define MATCOLS(i)
Definition: matpol.h:28
#define NULL
Definition: omList.c:10
static void mp_ElimBar(matrix, matrix, poly, int, int, const ring)
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
const CanonicalForm & w
Definition: facAbsFact.cc:55
#define SM_MULT
Definition: sparsmat.h:23
Variable x
Definition: cfModGcd.cc:4023
#define pNext(p)
Definition: monomials.h:43
#define p_GetCoeff(p, r)
Definition: monomials.h:57
#define SM_DIV
Definition: sparsmat.h:24
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:74
poly * mpRowAdr(int)
static void mpReplace(int j, int n, int &sign, int *perm)
perform replacement for pivot strategy in Bareiss algorithm change sign of determinant ...
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1013
END_NAMESPACE const void * p2
Definition: syzextra.cc:202
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete &#39;p&#39;
Definition: coeffs.h:459
poly TraceOfProd(matrix a, matrix b, int n, const ring R)
Definition: matpol.cc:302
#define MATROWS(i)
Definition: matpol.h:27
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
int perm[100]
const poly b
Definition: syzextra.cc:213
static FORCE_INLINE int n_Size(number n, const coeffs r)
return a non-negative measure for the complexity of n; return 0 only when n represents zero; (used fo...
Definition: coeffs.h:574
poly mpGetElem(int, int)
static poly mp_Leibnitz(matrix a, const ring)
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1020
static int sign(int x)
Definition: ring.cc:3333
void Werror(const char *fmt,...)
Definition: reporter.cc:189
static int mp_PreparePiv(matrix, int, int, const ring)
#define omAlloc0(size)
Definition: omAllocDecl.h:211
return result
Definition: facAbsBiFact.cc:76
#define MATELEM(mat, i, j)
Definition: matpol.h:29
void mpColSwap(int, int)