My Project  debian-1:4.1.1-p2+ds-4build1
matpol.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 
5 /*
6 * ABSTRACT:
7 */
8 
9 #include "misc/auxiliary.h"
10 
11 #include "omalloc/omalloc.h"
12 #include "misc/mylimits.h"
13 
14 #include "misc/intvec.h"
15 #include "coeffs/numbers.h"
16 
17 #include "reporter/reporter.h"
18 
19 
20 #include "monomials/ring.h"
21 #include "monomials/p_polys.h"
22 
23 #include "simpleideals.h"
24 #include "matpol.h"
25 #include "prCopy.h"
26 
27 #include "sparsmat.h"
28 
29 //omBin sip_sideal_bin = omGetSpecBin(sizeof(ip_smatrix));
30 /*0 implementation*/
31 
32 static poly mp_Exdiv ( poly m, poly d, poly vars, const ring);
33 static poly mp_Select (poly fro, poly what, const ring);
34 
35 /// create a r x c zero-matrix
36 matrix mpNew(int r, int c)
37 {
38  int rr=r;
39  if (rr<=0) rr=1;
40  //if ( (((int)(MAX_INT_VAL/sizeof(poly))) / rr) <= c)
41  //{
42  // Werror("internal error: creating matrix[%d][%d]",r,c);
43  // return NULL;
44  //}
46  rc->nrows = r;
47  rc->ncols = c;
48  rc->rank = r;
49  if ((c != 0)&&(r!=0))
50  {
51  size_t s=((size_t)r)*((size_t)c)*sizeof(poly);
52  rc->m = (poly*)omAlloc0(s);
53  //if (rc->m==NULL)
54  //{
55  // Werror("internal error: creating matrix[%d][%d]",r,c);
56  // return NULL;
57  //}
58  }
59  return rc;
60 }
61 
62 /// copies matrix a (from ring r to r)
63 matrix mp_Copy (matrix a, const ring r)
64 {
65  id_Test((ideal)a, r);
66  poly t;
67  int i, m=MATROWS(a), n=MATCOLS(a);
68  matrix b = mpNew(m, n);
69 
70  for (i=m*n-1; i>=0; i--)
71  {
72  t = a->m[i];
73  if (t!=NULL)
74  {
75  p_Normalize(t, r);
76  b->m[i] = p_Copy(t, r);
77  }
78  }
79  b->rank=a->rank;
80  return b;
81 }
82 
83 /// copies matrix a from rSrc into rDst
84 matrix mp_Copy(const matrix a, const ring rSrc, const ring rDst)
85 {
86  id_Test((ideal)a, rSrc);
87 
88  poly t;
89  int i, m=MATROWS(a), n=MATCOLS(a);
90 
91  matrix b = mpNew(m, n);
92 
93  for (i=m*n-1; i>=0; i--)
94  {
95  t = a->m[i];
96  if (t!=NULL)
97  {
98  b->m[i] = prCopyR_NoSort(t, rSrc, rDst);
99  p_Normalize(b->m[i], rDst);
100  }
101  }
102  b->rank=a->rank;
103 
104  id_Test((ideal)b, rDst);
105 
106  return b;
107 }
108 
109 
110 
111 /// make it a p * unit matrix
112 matrix mp_InitP(int r, int c, poly p, const ring R)
113 {
114  matrix rc = mpNew(r,c);
115  int i=si_min(r,c), n = c*(i-1)+i-1, inc = c+1;
116 
117  p_Normalize(p, R);
118  while (n>0)
119  {
120  rc->m[n] = p_Copy(p, R);
121  n -= inc;
122  }
123  rc->m[0]=p;
124  return rc;
125 }
126 
127 /// make it a v * unit matrix
128 matrix mp_InitI(int r, int c, int v, const ring R)
129 {
130  return mp_InitP(r, c, p_ISet(v, R), R);
131 }
132 
133 /// c = f*a
134 matrix mp_MultI(matrix a, int f, const ring R)
135 {
136  int k, n = a->nrows, m = a->ncols;
137  poly p = p_ISet(f, R);
138  matrix c = mpNew(n,m);
139 
140  for (k=m*n-1; k>0; k--)
141  c->m[k] = pp_Mult_qq(a->m[k], p, R);
142  c->m[0] = p_Mult_q(p_Copy(a->m[0], R), p, R);
143  return c;
144 }
145 
146 /// multiply a matrix 'a' by a poly 'p', destroy the args
147 matrix mp_MultP(matrix a, poly p, const ring R)
148 {
149  int k, n = a->nrows, m = a->ncols;
150 
151  p_Normalize(p, R);
152  for (k=m*n-1; k>0; k--)
153  {
154  if (a->m[k]!=NULL)
155  a->m[k] = p_Mult_q(a->m[k], p_Copy(p, R), R);
156  }
157  a->m[0] = p_Mult_q(a->m[0], p, R);
158  return a;
159 }
160 
161 /*2
162 * multiply a poly 'p' by a matrix 'a', destroy the args
163 */
164 matrix pMultMp(poly p, matrix a, const ring R)
165 {
166  int k, n = a->nrows, m = a->ncols;
167 
168  p_Normalize(p, R);
169  for (k=m*n-1; k>0; k--)
170  {
171  if (a->m[k]!=NULL)
172  a->m[k] = p_Mult_q(p_Copy(p, R), a->m[k], R);
173  }
174  a->m[0] = p_Mult_q(p, a->m[0], R);
175  return a;
176 }
177 
178 matrix mp_Add(matrix a, matrix b, const ring R)
179 {
180  int k, n = a->nrows, m = a->ncols;
181  if ((n != b->nrows) || (m != b->ncols))
182  {
183 /*
184 * Werror("cannot add %dx%d matrix and %dx%d matrix",
185 * m,n,b->cols(),b->rows());
186 */
187  return NULL;
188  }
189  matrix c = mpNew(n,m);
190  for (k=m*n-1; k>=0; k--)
191  c->m[k] = p_Add_q(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
192  return c;
193 }
194 
195 matrix mp_Sub(matrix a, matrix b, const ring R)
196 {
197  int k, n = a->nrows, m = a->ncols;
198  if ((n != b->nrows) || (m != b->ncols))
199  {
200 /*
201 * Werror("cannot sub %dx%d matrix and %dx%d matrix",
202 * m,n,b->cols(),b->rows());
203 */
204  return NULL;
205  }
206  matrix c = mpNew(n,m);
207  for (k=m*n-1; k>=0; k--)
208  c->m[k] = p_Sub(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
209  return c;
210 }
211 
212 matrix mp_Mult(matrix a, matrix b, const ring R)
213 {
214  int i, j, k;
215  int m = MATROWS(a);
216  int p = MATCOLS(a);
217  int q = MATCOLS(b);
218 
219  if (p!=MATROWS(b))
220  {
221 /*
222 * Werror("cannot multiply %dx%d matrix and %dx%d matrix",
223 * m,p,b->rows(),q);
224 */
225  return NULL;
226  }
227  matrix c = mpNew(m,q);
228 
229  for (i=1; i<=m; i++)
230  {
231  for (k=1; k<=p; k++)
232  {
233  poly aik;
234  if ((aik=MATELEM(a,i,k))!=NULL)
235  {
236  for (j=1; j<=q; j++)
237  {
238  poly bkj;
239  if ((bkj=MATELEM(b,k,j))!=NULL)
240  {
241  poly *cij=&(MATELEM(c,i,j));
242  poly s = pp_Mult_qq(aik /*MATELEM(a,i,k)*/, bkj/*MATELEM(b,k,j)*/, R);
243  if (/*MATELEM(c,i,j)*/ (*cij)==NULL) (*cij)=s;
244  else (*cij) = p_Add_q((*cij) /*MATELEM(c,i,j)*/ ,s, R);
245  }
246  }
247  }
248  // pNormalize(t);
249  // MATELEM(c,i,j) = t;
250  }
251  }
252  for(i=m*q-1;i>=0;i--) p_Normalize(c->m[i], R);
253  return c;
254 }
255 
256 matrix mp_Transp(matrix a, const ring R)
257 {
258  int i, j, r = MATROWS(a), c = MATCOLS(a);
259  poly *p;
260  matrix b = mpNew(c,r);
261 
262  p = b->m;
263  for (i=0; i<c; i++)
264  {
265  for (j=0; j<r; j++)
266  {
267  if (a->m[j*c+i]!=NULL) *p = p_Copy(a->m[j*c+i], R);
268  p++;
269  }
270  }
271  return b;
272 }
273 
274 /*2
275 *returns the trace of matrix a
276 */
277 poly mp_Trace ( matrix a, const ring R)
278 {
279  int i;
280  int n = (MATCOLS(a)<MATROWS(a)) ? MATCOLS(a) : MATROWS(a);
281  poly t = NULL;
282 
283  for (i=1; i<=n; i++)
284  t = p_Add_q(t, p_Copy(MATELEM(a,i,i), R), R);
285  return t;
286 }
287 
288 /*2
289 *returns the trace of the product of a and b
290 */
291 poly TraceOfProd ( matrix a, matrix b, int n, const ring R)
292 {
293  int i, j;
294  poly p, t = NULL;
295 
296  for (i=1; i<=n; i++)
297  {
298  for (j=1; j<=n; j++)
299  {
300  p = pp_Mult_qq(MATELEM(a,i,j), MATELEM(b,j,i), R);
301  t = p_Add_q(t, p, R);
302  }
303  }
304  return t;
305 }
306 
307 // #ifndef SIZE_OF_SYSTEM_PAGE
308 // #define SIZE_OF_SYSTEM_PAGE 4096
309 // #endif
310 
311 /*2
312 * corresponds to Maple's coeffs:
313 * var has to be the number of a variable
314 */
315 matrix mp_Coeffs (ideal I, int var, const ring R)
316 {
317  poly h,f;
318  int l, i, c, m=0;
319  /* look for maximal power m of x_var in I */
320  for (i=IDELEMS(I)-1; i>=0; i--)
321  {
322  f=I->m[i];
323  while (f!=NULL)
324  {
325  l=p_GetExp(f,var, R);
326  if (l>m) m=l;
327  pIter(f);
328  }
329  }
330  matrix co=mpNew((m+1)*I->rank,IDELEMS(I));
331  /* divide each monomial by a power of x_var,
332  * remember the power in l and the component in c*/
333  for (i=IDELEMS(I)-1; i>=0; i--)
334  {
335  f=I->m[i];
336  I->m[i]=NULL;
337  while (f!=NULL)
338  {
339  l=p_GetExp(f,var, R);
340  p_SetExp(f,var,0, R);
341  c=si_max((int)p_GetComp(f, R),1);
342  p_SetComp(f,0, R);
343  p_Setm(f, R);
344  /* now add the resulting monomial to co*/
345  h=pNext(f);
346  pNext(f)=NULL;
347  //MATELEM(co,c*(m+1)-l,i+1)
348  // =p_Add_q(MATELEM(co,c*(m+1)-l,i+1),f, R);
349  MATELEM(co,(c-1)*(m+1)+l+1,i+1)
350  =p_Add_q(MATELEM(co,(c-1)*(m+1)+l+1,i+1),f, R);
351  /* iterate f*/
352  f=h;
353  }
354  }
355  id_Delete(&I, R);
356  return co;
357 }
358 
359 /*2
360 * given the result c of mpCoeffs(ideal/module i, var)
361 * i of rank r
362 * build the matrix of the corresponding monomials in m
363 */
364 void mp_Monomials(matrix c, int r, int var, matrix m, const ring R)
365 {
366  /* clear contents of m*/
367  int k,l;
368  for (k=MATROWS(m);k>0;k--)
369  {
370  for(l=MATCOLS(m);l>0;l--)
371  {
372  p_Delete(&MATELEM(m,k,l), R);
373  }
374  }
375  omfreeSize((ADDRESS)m->m,MATROWS(m)*MATCOLS(m)*sizeof(poly));
376  /* allocate monoms in the right size r x MATROWS(c)*/
377  m->m=(poly*)omAlloc0(r*MATROWS(c)*sizeof(poly));
378  MATROWS(m)=r;
379  MATCOLS(m)=MATROWS(c);
380  m->rank=r;
381  /* the maximal power p of x_var: MATCOLS(m)=r*(p+1) */
382  int p=MATCOLS(m)/r-1;
383  /* fill in the powers of x_var=h*/
384  poly h=p_One(R);
385  for(k=r;k>0; k--)
386  {
387  MATELEM(m,k,k*(p+1))=p_One(R);
388  }
389  for(l=p;l>=0; l--)
390  {
391  p_SetExp(h,var,p-l, R);
392  p_Setm(h, R);
393  for(k=r;k>0; k--)
394  {
395  MATELEM(m,k,k*(p+1)-l)=p_Copy(h, R);
396  }
397  }
398  p_Delete(&h, R);
399 }
400 
401 matrix mp_CoeffProc (poly f, poly vars, const ring R)
402 {
403  assume(vars!=NULL);
404  poly sel, h;
405  int l, i;
406  int pos_of_1 = -1;
407  matrix co;
408 
409  if (f==NULL)
410  {
411  co = mpNew(2, 1);
412  MATELEM(co,1,1) = p_One(R);
413  MATELEM(co,2,1) = NULL;
414  return co;
415  }
416  sel = mp_Select(f, vars, R);
417  l = pLength(sel);
418  co = mpNew(2, l);
419 
421  {
422  for (i=l; i>=1; i--)
423  {
424  h = sel;
425  pIter(sel);
426  pNext(h)=NULL;
427  MATELEM(co,1,i) = h;
428  MATELEM(co,2,i) = NULL;
429  if (p_IsConstant(h, R)) pos_of_1 = i;
430  }
431  }
432  else
433  {
434  for (i=1; i<=l; i++)
435  {
436  h = sel;
437  pIter(sel);
438  pNext(h)=NULL;
439  MATELEM(co,1,i) = h;
440  MATELEM(co,2,i) = NULL;
441  if (p_IsConstant(h, R)) pos_of_1 = i;
442  }
443  }
444  while (f!=NULL)
445  {
446  i = 1;
447  loop
448  {
449  if (i!=pos_of_1)
450  {
451  h = mp_Exdiv(f, MATELEM(co,1,i),vars, R);
452  if (h!=NULL)
453  {
454  MATELEM(co,2,i) = p_Add_q(MATELEM(co,2,i), h, R);
455  break;
456  }
457  }
458  if (i == l)
459  {
460  // check monom 1 last:
461  if (pos_of_1 != -1)
462  {
463  h = mp_Exdiv(f, MATELEM(co,1,pos_of_1),vars, R);
464  if (h!=NULL)
465  {
466  MATELEM(co,2,pos_of_1) = p_Add_q(MATELEM(co,2,pos_of_1), h, R);
467  }
468  }
469  break;
470  }
471  i ++;
472  }
473  pIter(f);
474  }
475  return co;
476 }
477 
478 /*2
479 *exact divisor: let d == x^i*y^j, m is thought to have only one term;
480 * return m/d iff d divides m, and no x^k*y^l (k>i or l>j) divides m
481 * consider all variables in vars
482 */
483 static poly mp_Exdiv ( poly m, poly d, poly vars, const ring R)
484 {
485  int i;
486  poly h = p_Head(m, R);
487  for (i=1; i<=rVar(R); i++)
488  {
489  if (p_GetExp(vars,i, R) > 0)
490  {
491  if (p_GetExp(d,i, R) != p_GetExp(h,i, R))
492  {
493  p_Delete(&h, R);
494  return NULL;
495  }
496  p_SetExp(h,i,0, R);
497  }
498  }
499  p_Setm(h, R);
500  return h;
501 }
502 
503 void mp_Coef2(poly v, poly mon, matrix *c, matrix *m, const ring R)
504 {
505  poly* s;
506  poly p;
507  int sl,i,j;
508  int l=0;
509  poly sel=mp_Select(v,mon, R);
510 
511  p_Vec2Polys(sel,&s,&sl, R);
512  for (i=0; i<sl; i++)
513  l=si_max(l,pLength(s[i]));
514  *c=mpNew(sl,l);
515  *m=mpNew(sl,l);
516  poly h;
517  int isConst;
518  for (j=1; j<=sl;j++)
519  {
520  p=s[j-1];
521  if (p_IsConstant(p, R)) /*p != NULL */
522  {
523  isConst=-1;
524  i=l;
525  }
526  else
527  {
528  isConst=1;
529  i=1;
530  }
531  while(p!=NULL)
532  {
533  h = p_Head(p, R);
534  MATELEM(*m,j,i) = h;
535  i+=isConst;
536  p = p->next;
537  }
538  }
539  while (v!=NULL)
540  {
541  i = 1;
542  j = __p_GetComp(v, R);
543  loop
544  {
545  poly mp=MATELEM(*m,j,i);
546  if (mp!=NULL)
547  {
548  h = mp_Exdiv(v, mp /*MATELEM(*m,j,i)*/, mp, R);
549  if (h!=NULL)
550  {
551  p_SetComp(h,0, R);
552  MATELEM(*c,j,i) = p_Add_q(MATELEM(*c,j,i), h, R);
553  break;
554  }
555  }
556  if (i < l)
557  i++;
558  else
559  break;
560  }
561  v = v->next;
562  }
563 }
564 
565 int mp_Compare(matrix a, matrix b, const ring R)
566 {
567  if (MATCOLS(a)<MATCOLS(b)) return -1;
568  else if (MATCOLS(a)>MATCOLS(b)) return 1;
569  if (MATROWS(a)<MATROWS(b)) return -1;
570  else if (MATROWS(a)<MATROWS(b)) return 1;
571 
572  unsigned ii=MATCOLS(a)*MATROWS(a)-1;
573  unsigned j=0;
574  int r=0;
575  while (j<=ii)
576  {
577  r=p_Compare(a->m[j],b->m[j],R);
578  if (r!=0) return r;
579  j++;
580  }
581  return r;
582 }
583 
584 BOOLEAN mp_Equal(matrix a, matrix b, const ring R)
585 {
586  if ((MATCOLS(a)!=MATCOLS(b)) || (MATROWS(a)!=MATROWS(b)))
587  return FALSE;
588  int i=MATCOLS(a)*MATROWS(a)-1;
589  while (i>=0)
590  {
591  if (a->m[i]==NULL)
592  {
593  if (b->m[i]!=NULL) return FALSE;
594  }
595  else if (b->m[i]==NULL) return FALSE;
596  else if (p_Cmp(a->m[i],b->m[i], R)!=0) return FALSE;
597  i--;
598  }
599  i=MATCOLS(a)*MATROWS(a)-1;
600  while (i>=0)
601  {
602  if(!p_EqualPolys(a->m[i],b->m[i], R)) return FALSE;
603  i--;
604  }
605  return TRUE;
606 }
607 
608 /*2
609 * insert a monomial into a list, avoid duplicates
610 * arguments are destroyed
611 */
612 static poly p_Insert(poly p1, poly p2, const ring R)
613 {
614  poly a1, p, a2, a;
615  int c;
616 
617  if (p1==NULL) return p2;
618  if (p2==NULL) return p1;
619  a1 = p1;
620  a2 = p2;
621  a = p = p_One(R);
622  loop
623  {
624  c = p_Cmp(a1, a2, R);
625  if (c == 1)
626  {
627  a = pNext(a) = a1;
628  pIter(a1);
629  if (a1==NULL)
630  {
631  pNext(a) = a2;
632  break;
633  }
634  }
635  else if (c == -1)
636  {
637  a = pNext(a) = a2;
638  pIter(a2);
639  if (a2==NULL)
640  {
641  pNext(a) = a1;
642  break;
643  }
644  }
645  else
646  {
647  p_LmDelete(&a2, R);
648  a = pNext(a) = a1;
649  pIter(a1);
650  if (a1==NULL)
651  {
652  pNext(a) = a2;
653  break;
654  }
655  else if (a2==NULL)
656  {
657  pNext(a) = a1;
658  break;
659  }
660  }
661  }
662  p_LmDelete(&p, R);
663  return p;
664 }
665 
666 /*2
667 *if what == xy the result is the list of all different power products
668 * x^i*y^j (i, j >= 0) that appear in fro
669 */
670 static poly mp_Select (poly fro, poly what, const ring R)
671 {
672  int i;
673  poly h, res;
674  res = NULL;
675  while (fro!=NULL)
676  {
677  h = p_One(R);
678  for (i=1; i<=rVar(R); i++)
679  p_SetExp(h,i, p_GetExp(fro,i, R) * p_GetExp(what, i, R), R);
680  p_SetComp(h, p_GetComp(fro, R), R);
681  p_Setm(h, R);
682  res = p_Insert(h, res, R);
683  fro = fro->next;
684  }
685  return res;
686 }
687 
688 /*
689 *static void ppp(matrix a)
690 *{
691 * int j,i,r=a->nrows,c=a->ncols;
692 * for(j=1;j<=r;j++)
693 * {
694 * for(i=1;i<=c;i++)
695 * {
696 * if(MATELEM(a,j,i)!=NULL) PrintS("X");
697 * else PrintS("0");
698 * }
699 * PrintLn();
700 * }
701 *}
702 */
703 
704 static void mp_PartClean(matrix a, int lr, int lc, const ring R)
705 {
706  poly *q1;
707  int i,j;
708 
709  for (i=lr-1;i>=0;i--)
710  {
711  q1 = &(a->m)[i*a->ncols];
712  for (j=lc-1;j>=0;j--) if(q1[j]) p_Delete(&q1[j], R);
713  }
714 }
715 
716 BOOLEAN mp_IsDiagUnit(matrix U, const ring R)
717 {
718  if(MATROWS(U)!=MATCOLS(U))
719  return FALSE;
720  for(int i=MATCOLS(U);i>=1;i--)
721  {
722  for(int j=MATCOLS(U); j>=1; j--)
723  {
724  if (i==j)
725  {
726  if (!p_IsUnit(MATELEM(U,i,i), R)) return FALSE;
727  }
728  else if (MATELEM(U,i,j)!=NULL) return FALSE;
729  }
730  }
731  return TRUE;
732 }
733 
734 void iiWriteMatrix(matrix im, const char *n, int dim, const ring r, int spaces)
735 {
736  int i,ii = MATROWS(im)-1;
737  int j,jj = MATCOLS(im)-1;
738  poly *pp = im->m;
739 
740  for (i=0; i<=ii; i++)
741  {
742  for (j=0; j<=jj; j++)
743  {
744  if (spaces>0)
745  Print("%-*.*s",spaces,spaces," ");
746  if (dim == 2) Print("%s[%u,%u]=",n,i+1,j+1);
747  else if (dim == 1) Print("%s[%u]=",n,j+1);
748  else if (dim == 0) Print("%s=",n);
749  if ((i<ii)||(j<jj)) p_Write(*pp++, r);
750  else p_Write0(*pp, r);
751  }
752  }
753 }
754 
755 char * iiStringMatrix(matrix im, int dim, const ring r, char ch)
756 {
757  int i,ii = MATROWS(im);
758  int j,jj = MATCOLS(im);
759  poly *pp = im->m;
760  char ch_s[2];
761  ch_s[0]=ch;
762  ch_s[1]='\0';
763 
764  StringSetS("");
765 
766  for (i=0; i<ii; i++)
767  {
768  for (j=0; j<jj; j++)
769  {
770  p_String0(*pp++, r);
771  StringAppendS(ch_s);
772  if (dim > 1) StringAppendS("\n");
773  }
774  }
775  char *s=StringEndS();
776  s[strlen(s)- (dim > 1 ? 2 : 1)]='\0';
777  return s;
778 }
779 
780 void mp_Delete(matrix* a, const ring r)
781 {
782  id_Delete((ideal *) a, r);
783 }
784 
785 /*
786 * C++ classes for Bareiss algorithm
787 */
788 class row_col_weight
789 {
790  private:
791  int ym, yn;
792  public:
793  float *wrow, *wcol;
794  row_col_weight() : ym(0) {}
795  row_col_weight(int, int);
796  ~row_col_weight();
797 };
798 
800 {
801  ym = i;
802  yn = j;
803  wrow = (float *)omAlloc(i*sizeof(float));
804  wcol = (float *)omAlloc(j*sizeof(float));
805 }
807 {
808  if (ym!=0)
809  {
810  omFreeSize((ADDRESS)wcol, yn*sizeof(float));
811  omFreeSize((ADDRESS)wrow, ym*sizeof(float));
812  }
813 }
814 
815 /*2
816 * a submatrix M of a matrix X[m,n]:
817 * 0 <= i < s_m <= a_m
818 * 0 <= j < s_n <= a_n
819 * M = ( Xarray[qrow[i],qcol[j]] )
820 * if a_m = a_n and s_m = s_n
821 * det(X) = sign*div^(s_m-1)*det(M)
822 * resticted pivot for elimination
823 * 0 <= j < piv_s
824 */
825 class mp_permmatrix
826 {
827  private:
828  int a_m, a_n, s_m, s_n, sign, piv_s;
829  int *qrow, *qcol;
830  poly *Xarray;
831  ring _R;
832  void mpInitMat();
833  poly * mpRowAdr(int r)
834  { return &(Xarray[a_n*qrow[r]]); }
835  poly * mpColAdr(int c)
836  { return &(Xarray[qcol[c]]); }
837  void mpRowWeight(float *);
838  void mpColWeight(float *);
839  void mpRowSwap(int, int);
840  void mpColSwap(int, int);
841  public:
842  mp_permmatrix() : a_m(0) {}
843  mp_permmatrix(matrix, ring);
845  ~mp_permmatrix();
846  int mpGetRow();
847  int mpGetCol();
848  int mpGetRdim() { return s_m; }
849  int mpGetCdim() { return s_n; }
850  int mpGetSign() { return sign; }
851  void mpSetSearch(int s);
852  void mpSaveArray() { Xarray = NULL; }
853  poly mpGetElem(int, int);
854  void mpSetElem(poly, int, int);
855  void mpDelElem(int, int);
856  void mpElimBareiss(poly);
858  int mpPivotRow(row_col_weight *, int);
859  void mpToIntvec(intvec *);
860  void mpRowReorder();
861  void mpColReorder();
862 };
864 {
865  a_m = A->nrows;
866  a_n = A->ncols;
867  this->mpInitMat();
868  Xarray = A->m;
869  _R=R;
870 }
871 
873 {
874  poly p, *athis, *aM;
875  int i, j;
876 
877  _R=M->_R;
878  a_m = M->s_m;
879  a_n = M->s_n;
880  sign = M->sign;
881  this->mpInitMat();
882  Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
883  for (i=a_m-1; i>=0; i--)
884  {
885  athis = this->mpRowAdr(i);
886  aM = M->mpRowAdr(i);
887  for (j=a_n-1; j>=0; j--)
888  {
889  p = aM[M->qcol[j]];
890  if (p)
891  {
892  athis[j] = p_Copy(p,_R);
893  }
894  }
895  }
896 }
897 
899 {
900  int k;
901 
902  if (a_m != 0)
903  {
904  omFreeSize((ADDRESS)qrow,a_m*sizeof(int));
905  omFreeSize((ADDRESS)qcol,a_n*sizeof(int));
906  if (Xarray != NULL)
907  {
908  for (k=a_m*a_n-1; k>=0; k--)
909  p_Delete(&Xarray[k],_R);
910  omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
911  }
912  }
913 }
914 
915 
916 static float mp_PolyWeight(poly p, const ring r);
917 void mp_permmatrix::mpColWeight(float *wcol)
918 {
919  poly p, *a;
920  int i, j;
921  float count;
922 
923  for (j=s_n; j>=0; j--)
924  {
925  a = this->mpColAdr(j);
926  count = 0.0;
927  for(i=s_m; i>=0; i--)
928  {
929  p = a[a_n*qrow[i]];
930  if (p)
931  count += mp_PolyWeight(p,_R);
932  }
933  wcol[j] = count;
934  }
935 }
936 void mp_permmatrix::mpRowWeight(float *wrow)
937 {
938  poly p, *a;
939  int i, j;
940  float count;
941 
942  for (i=s_m; i>=0; i--)
943  {
944  a = this->mpRowAdr(i);
945  count = 0.0;
946  for(j=s_n; j>=0; j--)
947  {
948  p = a[qcol[j]];
949  if (p)
950  count += mp_PolyWeight(p,_R);
951  }
952  wrow[i] = count;
953  }
954 }
955 
956 void mp_permmatrix::mpRowSwap(int i1, int i2)
957 {
958  poly p, *a1, *a2;
959  int j;
960 
961  a1 = &(Xarray[a_n*i1]);
962  a2 = &(Xarray[a_n*i2]);
963  for (j=a_n-1; j>= 0; j--)
964  {
965  p = a1[j];
966  a1[j] = a2[j];
967  a2[j] = p;
968  }
969 }
970 
971 void mp_permmatrix::mpColSwap(int j1, int j2)
972 {
973  poly p, *a1, *a2;
974  int i, k = a_n*a_m;
975 
976  a1 = &(Xarray[j1]);
977  a2 = &(Xarray[j2]);
978  for (i=0; i< k; i+=a_n)
979  {
980  p = a1[i];
981  a1[i] = a2[i];
982  a2[i] = p;
983  }
984 }
986 {
987  int k;
988 
989  s_m = a_m;
990  s_n = a_n;
991  piv_s = 0;
992  qrow = (int *)omAlloc(a_m*sizeof(int));
993  qcol = (int *)omAlloc(a_n*sizeof(int));
994  for (k=a_m-1; k>=0; k--) qrow[k] = k;
995  for (k=a_n-1; k>=0; k--) qcol[k] = k;
996 }
997 
999 {
1000  int k, j, j1, j2;
1001 
1002  if (a_n > a_m)
1003  k = a_n - a_m;
1004  else
1005  k = 0;
1006  for (j=a_n-1; j>=k; j--)
1007  {
1008  j1 = qcol[j];
1009  if (j1 != j)
1010  {
1011  this->mpColSwap(j1, j);
1012  j2 = 0;
1013  while (qcol[j2] != j) j2++;
1014  qcol[j2] = j1;
1015  }
1016  }
1017 }
1020 {
1021  int k, i, i1, i2;
1022 
1023  if (a_m > a_n)
1024  k = a_m - a_n;
1025  else
1026  k = 0;
1027  for (i=a_m-1; i>=k; i--)
1028  {
1029  i1 = qrow[i];
1030  if (i1 != i)
1031  {
1032  this->mpRowSwap(i1, i);
1033  i2 = 0;
1034  while (qrow[i2] != i) i2++;
1035  qrow[i2] = i1;
1036  }
1037  }
1038 }
1039 
1040 /*
1041 * perform replacement for pivot strategy in Bareiss algorithm
1042 * change sign of determinant
1043 */
1044 static void mpReplace(int j, int n, int &sign, int *perm)
1045 {
1046  int k;
1047 
1048  if (j != n)
1049  {
1050  k = perm[n];
1051  perm[n] = perm[j];
1052  perm[j] = k;
1053  sign = -sign;
1054  }
1055 }
1056 /*2
1057 * pivot strategy for Bareiss algorithm
1058 */
1060 {
1061  poly p, *a;
1062  int i, j, iopt, jopt;
1063  float sum, f1, f2, fo, r, ro, lp;
1064  float *dr = C->wrow, *dc = C->wcol;
1065 
1066  fo = 1.0e20;
1067  ro = 0.0;
1068  iopt = jopt = -1;
1069 
1070  s_n--;
1071  s_m--;
1072  if (s_m == 0)
1073  return 0;
1074  if (s_n == 0)
1075  {
1076  for(i=s_m; i>=0; i--)
1077  {
1078  p = this->mpRowAdr(i)[qcol[0]];
1079  if (p)
1080  {
1081  f1 = mp_PolyWeight(p,_R);
1082  if (f1 < fo)
1083  {
1084  fo = f1;
1085  if (iopt >= 0)
1086  p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]),_R);
1087  iopt = i;
1088  }
1089  else
1090  p_Delete(&(this->mpRowAdr(i)[qcol[0]]),_R);
1091  }
1092  }
1093  if (iopt >= 0)
1094  mpReplace(iopt, s_m, sign, qrow);
1095  return 0;
1096  }
1097  this->mpRowWeight(dr);
1098  this->mpColWeight(dc);
1099  sum = 0.0;
1100  for(i=s_m; i>=0; i--)
1101  sum += dr[i];
1102  for(i=s_m; i>=0; i--)
1103  {
1104  r = dr[i];
1105  a = this->mpRowAdr(i);
1106  for(j=s_n; j>=0; j--)
1107  {
1108  p = a[qcol[j]];
1109  if (p)
1110  {
1111  lp = mp_PolyWeight(p,_R);
1112  ro = r - lp;
1113  f1 = ro * (dc[j]-lp);
1114  if (f1 != 0.0)
1115  {
1116  f2 = lp * (sum - ro - dc[j]);
1117  f2 += f1;
1118  }
1119  else
1120  f2 = lp-r-dc[j];
1121  if (f2 < fo)
1122  {
1123  fo = f2;
1124  iopt = i;
1125  jopt = j;
1126  }
1127  }
1128  }
1129  }
1130  if (iopt < 0)
1131  return 0;
1132  mpReplace(iopt, s_m, sign, qrow);
1133  mpReplace(jopt, s_n, sign, qcol);
1134  return 1;
1136 poly mp_permmatrix::mpGetElem(int r, int c)
1137 {
1138  return Xarray[a_n*qrow[r]+qcol[c]];
1139 }
1140 
1141 /*
1142 * the Bareiss-type elimination with division by div (div != NULL)
1143 */
1145 {
1146  poly piv, elim, q1, q2, *ap, *a;
1147  int i, j, jj;
1148 
1149  ap = this->mpRowAdr(s_m);
1150  piv = ap[qcol[s_n]];
1151  for(i=s_m-1; i>=0; i--)
1152  {
1153  a = this->mpRowAdr(i);
1154  elim = a[qcol[s_n]];
1155  if (elim != NULL)
1156  {
1157  elim = p_Neg(elim,_R);
1158  for (j=s_n-1; j>=0; j--)
1159  {
1160  q2 = NULL;
1161  jj = qcol[j];
1162  if (ap[jj] != NULL)
1163  {
1164  q2 = SM_MULT(ap[jj], elim, div,_R);
1165  if (a[jj] != NULL)
1166  {
1167  q1 = SM_MULT(a[jj], piv, div,_R);
1168  p_Delete(&a[jj],_R);
1169  q2 = p_Add_q(q2, q1, _R);
1170  }
1171  }
1172  else if (a[jj] != NULL)
1173  {
1174  q2 = SM_MULT(a[jj], piv, div, _R);
1175  }
1176  if ((q2!=NULL) && div)
1177  SM_DIV(q2, div, _R);
1178  a[jj] = q2;
1179  }
1180  p_Delete(&a[qcol[s_n]], _R);
1181  }
1182  else
1183  {
1184  for (j=s_n-1; j>=0; j--)
1185  {
1186  jj = qcol[j];
1187  if (a[jj] != NULL)
1188  {
1189  q2 = SM_MULT(a[jj], piv, div, _R);
1190  p_Delete(&a[jj], _R);
1191  if (div)
1192  SM_DIV(q2, div, _R);
1193  a[jj] = q2;
1194  }
1195  }
1196  }
1197  }
1198 }
1199 /*
1200 * weigth of a polynomial, for pivot strategy
1201 */
1202 static float mp_PolyWeight(poly p, const ring r)
1203 {
1204  int i;
1205  float res;
1206 
1207  if (pNext(p) == NULL)
1208  {
1209  res = (float)n_Size(pGetCoeff(p),r->cf);
1210  for (i=r->N;i>0;i--)
1211  {
1212  if(p_GetExp(p,i,r)!=0)
1213  {
1214  res += 2.0;
1215  break;
1216  }
1217  }
1218  }
1219  else
1220  {
1221  res = 0.0;
1222  do
1223  {
1224  res += (float)n_Size(pGetCoeff(p),r->cf)+2.0;
1225  pIter(p);
1226  }
1227  while (p);
1228  }
1229  return res;
1230 }
1231 /*
1232 * find best row
1233 */
1234 static int mp_PivBar(matrix a, int lr, int lc, const ring r)
1235 {
1236  float f1, f2;
1237  poly *q1;
1238  int i,j,io;
1239 
1240  io = -1;
1241  f1 = 1.0e30;
1242  for (i=lr-1;i>=0;i--)
1243  {
1244  q1 = &(a->m)[i*a->ncols];
1245  f2 = 0.0;
1246  for (j=lc-1;j>=0;j--)
1247  {
1248  if (q1[j]!=NULL)
1249  f2 += mp_PolyWeight(q1[j],r);
1250  }
1251  if ((f2!=0.0) && (f2<f1))
1252  {
1253  f1 = f2;
1254  io = i;
1255  }
1256  }
1257  if (io<0) return 0;
1258  else return io+1;
1259 }
1261 static void mpSwapRow(matrix a, int pos, int lr, int lc)
1262 {
1263  poly sw;
1264  int j;
1265  poly* a2 = a->m;
1266  poly* a1 = &a2[a->ncols*(pos-1)];
1267 
1268  a2 = &a2[a->ncols*(lr-1)];
1269  for (j=lc-1; j>=0; j--)
1270  {
1271  sw = a1[j];
1272  a1[j] = a2[j];
1273  a2[j] = sw;
1274  }
1275 }
1276 
1277 /*2
1278 * prepare one step of 'Bareiss' algorithm
1279 * for application in minor
1280 */
1281 static int mp_PrepareRow (matrix a, int lr, int lc, const ring R)
1282 {
1283  int r;
1284 
1285  r = mp_PivBar(a,lr,lc,R);
1286  if(r==0) return 0;
1287  if(r<lr) mpSwapRow(a, r, lr, lc);
1288  return 1;
1289 }
1290 
1291 /*
1292 * find pivot in the last row
1293 */
1294 static int mp_PivRow(matrix a, int lr, int lc, const ring r)
1295 {
1296  float f1, f2;
1297  poly *q1;
1298  int j,jo;
1299 
1300  jo = -1;
1301  f1 = 1.0e30;
1302  q1 = &(a->m)[(lr-1)*a->ncols];
1303  for (j=lc-1;j>=0;j--)
1304  {
1305  if (q1[j]!=NULL)
1306  {
1307  f2 = mp_PolyWeight(q1[j],r);
1308  if (f2<f1)
1309  {
1310  f1 = f2;
1311  jo = j;
1312  }
1313  }
1314  }
1315  if (jo<0) return 0;
1316  else return jo+1;
1317 }
1319 static void mpSwapCol(matrix a, int pos, int lr, int lc)
1320 {
1321  poly sw;
1322  int j;
1323  poly* a2 = a->m;
1324  poly* a1 = &a2[pos-1];
1325 
1326  a2 = &a2[lc-1];
1327  for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
1328  {
1329  sw = a1[j];
1330  a1[j] = a2[j];
1331  a2[j] = sw;
1332  }
1333 }
1334 
1335 /*2
1336 * prepare one step of 'Bareiss' algorithm
1337 * for application in minor
1338 */
1339 static int mp_PreparePiv (matrix a, int lr, int lc,const ring r)
1340 {
1341  int c;
1342 
1343  c = mp_PivRow(a, lr, lc,r);
1344  if(c==0) return 0;
1345  if(c<lc) mpSwapCol(a, c, lr, lc);
1346  return 1;
1347 }
1349 static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
1350 {
1351  int r=lr-1, c=lc-1;
1352  poly *b = a0->m, *x = re->m;
1353  poly piv, elim, q1, *ap, *a, *q;
1354  int i, j;
1355 
1356  ap = &b[r*a0->ncols];
1357  piv = ap[c];
1358  for(j=c-1; j>=0; j--)
1359  if (ap[j] != NULL) ap[j] = p_Neg(ap[j],R);
1360  for(i=r-1; i>=0; i--)
1361  {
1362  a = &b[i*a0->ncols];
1363  q = &x[i*re->ncols];
1364  if (a[c] != NULL)
1365  {
1366  elim = a[c];
1367  for (j=c-1; j>=0; j--)
1368  {
1369  q1 = NULL;
1370  if (a[j] != NULL)
1371  {
1372  q1 = sm_MultDiv(a[j], piv, div,R);
1373  if (ap[j] != NULL)
1374  {
1375  poly q2 = sm_MultDiv(ap[j], elim, div, R);
1376  q1 = p_Add_q(q1,q2,R);
1377  }
1378  }
1379  else if (ap[j] != NULL)
1380  q1 = sm_MultDiv(ap[j], elim, div, R);
1381  if (q1 != NULL)
1382  {
1383  if (div)
1384  sm_SpecialPolyDiv(q1, div,R);
1385  q[j] = q1;
1386  }
1387  }
1388  }
1389  else
1390  {
1391  for (j=c-1; j>=0; j--)
1392  {
1393  if (a[j] != NULL)
1394  {
1395  q1 = sm_MultDiv(a[j], piv, div, R);
1396  if (div)
1397  sm_SpecialPolyDiv(q1, div, R);
1398  q[j] = q1;
1399  }
1400  }
1401  }
1402  }
1403 }
1404 
1405 /*2*/
1406 /// entries of a are minors and go to result (only if not in R)
1407 void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
1408  ideal R, const ring)
1409 {
1410  poly *q1;
1411  int e=IDELEMS(result);
1412  int i,j;
1413 
1414  if (R != NULL)
1415  {
1416  for (i=r-1;i>=0;i--)
1417  {
1418  q1 = &(a->m)[i*a->ncols];
1419  //for (j=c-1;j>=0;j--)
1420  //{
1421  // if (q1[j]!=NULL) q1[j] = kNF(R,currRing->qideal,q1[j]);
1422  //}
1423  }
1424  }
1425  for (i=r-1;i>=0;i--)
1426  {
1427  q1 = &(a->m)[i*a->ncols];
1428  for (j=c-1;j>=0;j--)
1429  {
1430  if (q1[j]!=NULL)
1431  {
1432  if (elems>=e)
1433  {
1434  pEnlargeSet(&(result->m),e,e);
1435  e += e;
1436  IDELEMS(result) =e;
1437  }
1438  result->m[elems] = q1[j];
1439  q1[j] = NULL;
1440  elems++;
1441  }
1442  }
1443  }
1444 }
1445 /*
1446 // from linalg_from_matpol.cc: TODO: compare with above & remove...
1447 void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
1448  ideal R, const ring R)
1449 {
1450  poly *q1;
1451  int e=IDELEMS(result);
1452  int i,j;
1453 
1454  if (R != NULL)
1455  {
1456  for (i=r-1;i>=0;i--)
1457  {
1458  q1 = &(a->m)[i*a->ncols];
1459  for (j=c-1;j>=0;j--)
1460  {
1461  if (q1[j]!=NULL) q1[j] = kNF(R,currRing->qideal,q1[j]);
1462  }
1463  }
1464  }
1465  for (i=r-1;i>=0;i--)
1466  {
1467  q1 = &(a->m)[i*a->ncols];
1468  for (j=c-1;j>=0;j--)
1469  {
1470  if (q1[j]!=NULL)
1471  {
1472  if (elems>=e)
1473  {
1474  if(e<SIZE_OF_SYSTEM_PAGE)
1475  {
1476  pEnlargeSet(&(result->m),e,e);
1477  e += e;
1478  }
1479  else
1480  {
1481  pEnlargeSet(&(result->m),e,SIZE_OF_SYSTEM_PAGE);
1482  e += SIZE_OF_SYSTEM_PAGE;
1483  }
1484  IDELEMS(result) =e;
1485  }
1486  result->m[elems] = q1[j];
1487  q1[j] = NULL;
1488  elems++;
1489  }
1490  }
1491  }
1492 }
1493 */
1495 static void mpFinalClean(matrix a)
1496 {
1497  omFreeSize((ADDRESS)a->m,a->nrows*a->ncols*sizeof(poly));
1499 }
1500 
1501 /*2*/
1502 /// produces recursively the ideal of all arxar-minors of a
1503 void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
1504  poly barDiv, ideal R, const ring r)
1505 {
1506  int k;
1507  int kr=lr-1,kc=lc-1;
1508  matrix nextLevel=mpNew(kr,kc);
1509 
1510  loop
1511  {
1512 /*--- look for an optimal row and bring it to last position ------------*/
1513  if(mp_PrepareRow(a,lr,lc,r)==0) break;
1514 /*--- now take all pivots from the last row ------------*/
1515  k = lc;
1516  loop
1517  {
1518  if(mp_PreparePiv(a,lr,k,r)==0) break;
1519  mp_ElimBar(a,nextLevel,barDiv,lr,k,r);
1520  k--;
1521  if (ar>1)
1522  {
1523  mp_RecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R,r);
1524  mp_PartClean(nextLevel,kr,k, r);
1525  }
1526  else mp_MinorToResult(result,elems,nextLevel,kr,k,R,r);
1527  if (ar>k-1) break;
1528  }
1529  if (ar>=kr) break;
1530 /*--- now we have to take out the last row...------------*/
1531  lr = kr;
1532  kr--;
1533  }
1534  mpFinalClean(nextLevel);
1535 }
1536 /*
1537 // from linalg_from_matpol.cc: TODO: compare with above & remove...
1538 void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
1539  poly barDiv, ideal R, const ring R)
1540 {
1541  int k;
1542  int kr=lr-1,kc=lc-1;
1543  matrix nextLevel=mpNew(kr,kc);
1544 
1545  loop
1546  {
1547 // --- look for an optimal row and bring it to last position ------------
1548  if(mpPrepareRow(a,lr,lc)==0) break;
1549 // --- now take all pivots from the last row ------------
1550  k = lc;
1551  loop
1552  {
1553  if(mpPreparePiv(a,lr,k)==0) break;
1554  mpElimBar(a,nextLevel,barDiv,lr,k);
1555  k--;
1556  if (ar>1)
1557  {
1558  mpRecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R);
1559  mpPartClean(nextLevel,kr,k);
1560  }
1561  else mpMinorToResult(result,elems,nextLevel,kr,k,R);
1562  if (ar>k-1) break;
1563  }
1564  if (ar>=kr) break;
1565 // --- now we have to take out the last row...------------
1566  lr = kr;
1567  kr--;
1568  }
1569  mpFinalClean(nextLevel);
1570 }
1571 */
1572 
1573 /*2*/
1574 /// returns the determinant of the matrix m;
1575 /// uses Bareiss algorithm
1576 poly mp_DetBareiss (matrix a, const ring r)
1577 {
1578  int s;
1579  poly div, res;
1580  if (MATROWS(a) != MATCOLS(a))
1581  {
1582  Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
1583  return NULL;
1584  }
1585  matrix c = mp_Copy(a,r);
1586  mp_permmatrix *Bareiss = new mp_permmatrix(c,r);
1587  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1588 
1589  /* Bareiss */
1590  div = NULL;
1591  while(Bareiss->mpPivotBareiss(&w))
1592  {
1593  Bareiss->mpElimBareiss(div);
1594  div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1595  }
1596  Bareiss->mpRowReorder();
1597  Bareiss->mpColReorder();
1598  Bareiss->mpSaveArray();
1599  s = Bareiss->mpGetSign();
1600  delete Bareiss;
1601 
1602  /* result */
1603  res = MATELEM(c,1,1);
1604  MATELEM(c,1,1) = NULL;
1605  id_Delete((ideal *)&c,r);
1606  if (s < 0)
1607  res = p_Neg(res,r);
1608  return res;
1609 }
1610 /*
1611 // from linalg_from_matpol.cc: TODO: compare with above & remove...
1612 poly mp_DetBareiss (matrix a, const ring R)
1613 {
1614  int s;
1615  poly div, res;
1616  if (MATROWS(a) != MATCOLS(a))
1617  {
1618  Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
1619  return NULL;
1620  }
1621  matrix c = mp_Copy(a, R);
1622  mp_permmatrix *Bareiss = new mp_permmatrix(c, R);
1623  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1624 
1625  // Bareiss
1626  div = NULL;
1627  while(Bareiss->mpPivotBareiss(&w))
1628  {
1629  Bareiss->mpElimBareiss(div);
1630  div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1631  }
1632  Bareiss->mpRowReorder();
1633  Bareiss->mpColReorder();
1634  Bareiss->mpSaveArray();
1635  s = Bareiss->mpGetSign();
1636  delete Bareiss;
1637 
1638  // result
1639  res = MATELEM(c,1,1);
1640  MATELEM(c,1,1) = NULL;
1641  id_Delete((ideal *)&c, R);
1642  if (s < 0)
1643  res = p_Neg(res, R);
1644  return res;
1645 }
1646 */
1647 
1648 /*2
1649 * compute all ar-minors of the matrix a
1650 */
1651 matrix mp_Wedge(matrix a, int ar, const ring R)
1652 {
1653  int i,j,k,l;
1654  int *rowchoise,*colchoise;
1655  BOOLEAN rowch,colch;
1656  matrix result;
1657  matrix tmp;
1658  poly p;
1659 
1660  i = binom(a->nrows,ar);
1661  j = binom(a->ncols,ar);
1662 
1663  rowchoise=(int *)omAlloc(ar*sizeof(int));
1664  colchoise=(int *)omAlloc(ar*sizeof(int));
1665  result = mpNew(i,j);
1666  tmp = mpNew(ar,ar);
1667  l = 1; /* k,l:the index in result*/
1668  idInitChoise(ar,1,a->nrows,&rowch,rowchoise);
1669  while (!rowch)
1670  {
1671  k=1;
1672  idInitChoise(ar,1,a->ncols,&colch,colchoise);
1673  while (!colch)
1674  {
1675  for (i=1; i<=ar; i++)
1676  {
1677  for (j=1; j<=ar; j++)
1678  {
1679  MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1680  }
1681  }
1682  p = mp_DetBareiss(tmp, R);
1683  if ((k+l) & 1) p=p_Neg(p, R);
1684  MATELEM(result,l,k) = p;
1685  k++;
1686  idGetNextChoise(ar,a->ncols,&colch,colchoise);
1687  }
1688  idGetNextChoise(ar,a->nrows,&rowch,rowchoise);
1689  l++;
1690  }
1691 
1692  /*delete the matrix tmp*/
1693  for (i=1; i<=ar; i++)
1694  {
1695  for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1696  }
1697  id_Delete((ideal *) &tmp, R);
1698  return (result);
1699 }
1701 static void p_DecomposeComp(poly p, poly *a, int l, const ring r)
1702 {
1703  poly h=p;
1704  while(h!=NULL)
1705  {
1706  poly hh=pNext(h);
1707  pNext(h)=a[__p_GetComp(h,r)-1];
1708  a[__p_GetComp(h,r)-1]=h;
1709  p_SetComp(h,0,r);
1710  p_SetmComp(h,r);
1711  h=hh;
1712  }
1713  for(int i=0;i<l;i++)
1714  {
1715  if(a[i]!=NULL) a[i]=pReverse(a[i]);
1716  }
1717 }
1718 // helper for mp_Tensor
1719 // destroyes f, keeps B
1720 static ideal mp_MultAndShift(poly f, ideal B, int s, const ring r)
1721 {
1722  assume(f!=NULL);
1723  ideal res=idInit(IDELEMS(B),B->rank+s);
1724  int q=IDELEMS(B); // p x q
1725  for(int j=0;j<q;j++)
1726  {
1727  poly h=pp_Mult_qq(f,B->m[j],r);
1728  if (h!=NULL)
1729  {
1730  if (s>0) p_Shift(&h,s,r);
1731  res->m[j]=h;
1732  }
1733  }
1734  p_Delete(&f,r);
1735  return res;
1736 }
1737 // helper for mp_Tensor
1738 // updates res, destroyes contents of sm
1739 static void mp_AddSubMat(ideal res, ideal sm, int col, const ring r)
1740 {
1741  for(int i=0;i<IDELEMS(sm);i++)
1742  {
1743  res->m[col+i]=p_Add_q(res->m[col+i],sm->m[i],r);
1744  sm->m[i]=NULL;
1745  }
1746 }
1748 ideal mp_Tensor(ideal A, ideal B, const ring r)
1749 {
1750  // size of the result n*q x m*p
1751  int n=IDELEMS(A); // m x n
1752  int m=A->rank;
1753  int q=IDELEMS(B); // p x q
1754  int p=B->rank;
1755  ideal res=idInit(n*q,m*p);
1756  poly *a=(poly*)omAlloc(m*sizeof(poly));
1757  for(int i=0; i<n; i++)
1758  {
1759  memset(a,0,m*sizeof(poly));
1760  p_DecomposeComp(p_Copy(A->m[i],r),a,m,r);
1761  for(int j=0;j<m;j++)
1762  {
1763  if (a[j]!=NULL)
1764  {
1765  ideal sm=mp_MultAndShift(a[j], // A_i_j
1766  B,
1767  j*p, // shift j*p down
1768  r);
1769  mp_AddSubMat(res,sm,i*q,r); // add this columns to col i*q ff
1770  id_Delete(&sm,r); // delete the now empty ideal
1771  }
1772  }
1773  }
1774  omFreeSize(a,m*sizeof(poly));
1775  return res;
1776 }
1777 
dim
int dim(ideal I, ring r)
Definition: tropicalStrategy.cc:22
si_min
static int si_min(const int a, const int b)
Definition: auxiliary.h:139
FALSE
#define FALSE
Definition: auxiliary.h:94
sip_sideal_bin
omBin sip_sideal_bin
Definition: simpleideals.cc:28
omalloc.h
matrix
ip_smatrix * matrix
Definition: matpol.h:30
p_GetComp
#define p_GetComp(p, r)
Definition: monomials.h:68
ip_smatrix
Definition: matpol.h:13
StringAppendS
void StringAppendS(const char *st)
Definition: reporter.cc:106
row_col_weight::yn
int yn
Definition: matpol.cc:790
sparsmat.h
p_GetExp
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent @Note: the integer VarOffset encodes:
Definition: p_polys.h:457
j
int j
Definition: facHensel.cc:105
f
FILE * f
Definition: checklibs.c:9
mp_permmatrix::mp_permmatrix
mp_permmatrix()
Definition: matpol.cc:841
p_Normalize
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3706
mp_permmatrix::mpGetCol
int mpGetCol()
k
int k
Definition: cfEzgcd.cc:92
p_Write0
void p_Write0(poly p, ring lmRing, ring tailRing)
Definition: polys0.cc:193
mp_CoeffProc
matrix mp_CoeffProc(poly f, poly vars, const ring R)
Definition: matpol.cc:400
mp_PivBar
static int mp_PivBar(matrix a, int lr, int lc, const ring r)
Definition: matpol.cc:1233
x
Variable x
Definition: cfModGcd.cc:4023
MATELEM
#define MATELEM(mat, i, j)
Definition: matpol.h:28
mp_MultAndShift
static ideal mp_MultAndShift(poly f, ideal B, int s, const ring r)
Definition: matpol.cc:1719
result
return result
Definition: facAbsBiFact.cc:76
mp_permmatrix::mpToIntvec
void mpToIntvec(intvec *)
mp_permmatrix::mpPivotBareiss
int mpPivotBareiss(row_col_weight *)
Definition: matpol.cc:1058
mp_Coeffs
matrix mp_Coeffs(ideal I, int var, const ring R)
corresponds to Maple's coeffs: var has to be the number of a variable
Definition: matpol.cc:314
SM_MULT
#define SM_MULT
Definition: sparsmat.h:22
mp_Trace
poly mp_Trace(matrix a, const ring R)
Definition: matpol.cc:276
pEnlargeSet
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3630
ADDRESS
void * ADDRESS
Definition: auxiliary.h:133
mp_ElimBar
static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
Definition: matpol.cc:1348
p_Head
static poly p_Head(poly p, const ring r)
Definition: p_polys.h:810
mp_permmatrix::a_m
int a_m
Definition: matpol.cc:827
mp_permmatrix::mpColWeight
void mpColWeight(float *)
Definition: matpol.cc:916
p_Neg
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1028
mp_Select
static poly mp_Select(poly fro, poly what, const ring)
Definition: matpol.cc:669
mp_Exdiv
static poly mp_Exdiv(poly m, poly d, poly vars, const ring)
Definition: matpol.cc:482
mp_PrepareRow
static int mp_PrepareRow(matrix a, int lr, int lc, const ring R)
Definition: matpol.cc:1280
mp_permmatrix::mpColAdr
poly * mpColAdr(int c)
Definition: matpol.cc:834
mp_InitP
matrix mp_InitP(int r, int c, poly p, const ring R)
make it a p * unit matrix
Definition: matpol.cc:111
simpleideals.h
idGetNextChoise
void idGetNextChoise(int r, int end, BOOLEAN *endch, int *choise)
Definition: simpleideals.cc:854
mp_permmatrix::s_m
int s_m
Definition: matpol.cc:827
row_col_weight::~row_col_weight
~row_col_weight()
Definition: matpol.cc:805
sign
static int sign(int x)
Definition: ring.cc:3327
__p_GetComp
#define __p_GetComp(p, r)
Definition: monomials.h:67
mp_permmatrix::mpPivotRow
int mpPivotRow(row_col_weight *, int)
mp_RecMin
void mp_RecMin(int ar, ideal result, int &elems, matrix a, int lr, int lc, poly barDiv, ideal R, const ring r)
produces recursively the ideal of all arxar-minors of a
Definition: matpol.cc:1502
mp_PartClean
static void mp_PartClean(matrix a, int lr, int lc, const ring R)
Definition: matpol.cc:703
omAllocBin
#define omAllocBin(bin)
Definition: omAllocDecl.h:203
auxiliary.h
mp_permmatrix::mpRowWeight
void mpRowWeight(float *)
Definition: matpol.cc:935
Variable::next
Variable next() const
Definition: factory.h:137
reporter.h
mp_Tensor
ideal mp_Tensor(ideal A, ideal B, const ring r)
Definition: matpol.cc:1747
StringEndS
char * StringEndS()
Definition: reporter.cc:150
loop
#define loop
Definition: structs.h:77
mp_permmatrix::mpRowAdr
poly * mpRowAdr(int r)
Definition: matpol.cc:832
w
const CanonicalForm & w
Definition: facAbsFact.cc:55
mp_permmatrix::mpGetElem
poly mpGetElem(int, int)
Definition: matpol.cc:1135
b
CanonicalForm b
Definition: cfModGcd.cc:4044
iiWriteMatrix
void iiWriteMatrix(matrix im, const char *n, int dim, const ring r, int spaces)
set spaces to zero by default
Definition: matpol.cc:733
mp_permmatrix::qcol
int * qcol
Definition: matpol.cc:828
ap
Definition: ap.h:35
mpFinalClean
static void mpFinalClean(matrix a)
Definition: matpol.cc:1494
mp_permmatrix::mpInitMat
void mpInitMat()
Definition: matpol.cc:984
mp_permmatrix::mpColReorder
void mpColReorder()
Definition: matpol.cc:997
p_SetmComp
#define p_SetmComp
Definition: p_polys.h:235
mp_permmatrix::qrow
int * qrow
Definition: matpol.cc:828
p_SetExp
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent @Note: VarOffset encodes the position in p->exp
Definition: p_polys.h:476
pLength
static unsigned pLength(poly a)
Definition: p_polys.h:183
mp_permmatrix::mpElimBareiss
void mpElimBareiss(poly)
Definition: matpol.cc:1143
mp_MinorToResult
void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c, ideal R, const ring)
entries of a are minors and go to result (only if not in R)
Definition: matpol.cc:1406
row_col_weight
Definition: matpol.cc:787
mpSwapCol
static void mpSwapCol(matrix a, int pos, int lr, int lc)
Definition: matpol.cc:1318
mp_permmatrix::_R
ring _R
Definition: matpol.cc:830
for
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:65
p_Copy
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:797
mp_permmatrix::mpGetRdim
int mpGetRdim()
Definition: matpol.cc:847
mpReplace
static void mpReplace(int j, int n, int &sign, int *perm)
Definition: matpol.cc:1043
p_DecomposeComp
static void p_DecomposeComp(poly p, poly *a, int l, const ring r)
Definition: matpol.cc:1700
rVar
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:582
TRUE
#define TRUE
Definition: auxiliary.h:98
i
int i
Definition: cfEzgcd.cc:125
ip_smatrix::nrows
int nrows
Definition: matpol.h:20
rHasLocalOrMixedOrdering
BOOLEAN rHasLocalOrMixedOrdering(const ring r)
Definition: ring.h:751
id_Delete
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
Definition: simpleideals.cc:113
p_Sub
poly p_Sub(poly p1, poly p2, const ring r)
Definition: p_polys.cc:1934
res
CanonicalForm res
Definition: facAbsFact.cc:64
mp_Wedge
matrix mp_Wedge(matrix a, int ar, const ring R)
Definition: matpol.cc:1650
prCopy.h
matpol.h
SM_DIV
#define SM_DIV
Definition: sparsmat.h:23
M
#define M
Definition: sirandom.c:24
omFreeSize
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:258
mp_permmatrix::s_n
int s_n
Definition: matpol.cc:827
lc
CanonicalForm lc(const CanonicalForm &f)
Definition: canonicalform.h:297
BOOLEAN
int BOOLEAN
Definition: auxiliary.h:85
p_Cmp
static int p_Cmp(poly p1, poly p2, ring r)
Definition: p_polys.h:1642
omfreeSize
#define omfreeSize(addr, size)
Definition: omAllocDecl.h:234
mp_Transp
matrix mp_Transp(matrix a, const ring R)
Definition: matpol.cc:255
row_col_weight::wcol
float * wcol
Definition: matpol.cc:792
mp_MultI
matrix mp_MultI(matrix a, int f, const ring R)
c = f*a
Definition: matpol.cc:133
mp_permmatrix::~mp_permmatrix
~mp_permmatrix()
Definition: matpol.cc:897
sm_MultDiv
poly sm_MultDiv(poly a, poly b, const poly c, const ring R)
Definition: sparsmat.cc:1813
h
static Poly * h
Definition: janet.cc:972
mp_PreparePiv
static int mp_PreparePiv(matrix a, int lr, int lc, const ring r)
Definition: matpol.cc:1338
ip_smatrix::m
poly * m
Definition: matpol.h:18
p_LmDelete
static void p_LmDelete(poly p, const ring r)
Definition: p_polys.h:698
row_col_weight::row_col_weight
row_col_weight()
Definition: matpol.cc:793
p_String0
void p_String0(poly p, ring lmRing, ring tailRing)
print p according to ShortOut in lmRing & tailRing
Definition: polys0.cc:133
pp_Mult_qq
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1078
p_Write
void p_Write(poly p, ring lmRing, ring tailRing)
Definition: polys0.cc:203
intvec
Definition: intvec.h:16
pIter
#define pIter(p)
Definition: monomials.h:41
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:208
p_polys.h
mp_permmatrix
Definition: matpol.cc:824
row_col_weight::wrow
float * wrow
Definition: matpol.cc:792
pMultMp
matrix pMultMp(poly p, matrix a, const ring R)
Definition: matpol.cc:163
p_Compare
int p_Compare(const poly a, const poly b, const ring R)
Definition: p_polys.cc:4790
mp_permmatrix::mpDelElem
void mpDelElem(int, int)
sm_SpecialPolyDiv
void sm_SpecialPolyDiv(poly a, poly b, const ring R)
Definition: sparsmat.cc:1894
mp_Delete
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:779
pp
CanonicalForm pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:253
p_IsUnit
static BOOLEAN p_IsUnit(const poly p, const ring r)
Definition: p_polys.h:1928
mp_permmatrix::mpColSwap
void mpColSwap(int, int)
Definition: matpol.cc:970
intvec.h
mpNew
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:35
p_Shift
void p_Shift(poly *p, int i, const ring r)
shifts components of the vector p by i
Definition: p_polys.cc:4586
mpSwapRow
static void mpSwapRow(matrix a, int pos, int lr, int lc)
Definition: matpol.cc:1260
mp_Coef2
void mp_Coef2(poly v, poly mon, matrix *c, matrix *m, const ring R)
corresponds to Macauley's coef: the exponent vector of vars has to contain the variables,...
Definition: matpol.cc:502
mp_permmatrix::mpRowSwap
void mpRowSwap(int, int)
Definition: matpol.cc:955
mp_permmatrix::mpRowReorder
void mpRowReorder()
Definition: matpol.cc:1018
mp_InitI
matrix mp_InitI(int r, int c, int v, const ring R)
make it a v * unit matrix
Definition: matpol.cc:127
mp_permmatrix::mpGetRow
int mpGetRow()
mp_permmatrix::mpGetCdim
int mpGetCdim()
Definition: matpol.cc:848
mp_Mult
matrix mp_Mult(matrix a, matrix b, const ring R)
Definition: matpol.cc:211
div
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
ring.h
p_Delete
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:842
p_Add_q
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
p_One
poly p_One(const ring r)
Definition: p_polys.cc:1302
p_EqualPolys
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4392
StringSetS
void StringSetS(const char *st)
Definition: reporter.cc:127
mp_MultP
matrix mp_MultP(matrix a, poly p, const ring R)
multiply a matrix 'a' by a poly 'p', destroy the args
Definition: matpol.cc:146
si_max
static int si_max(const int a, const int b)
Definition: auxiliary.h:138
mp_permmatrix::sign
int sign
Definition: matpol.cc:827
Print
#define Print
Definition: emacs.cc:79
B
b *CanonicalForm B
Definition: facBivar.cc:52
mylimits.h
mp_IsDiagUnit
BOOLEAN mp_IsDiagUnit(matrix U, const ring R)
Definition: matpol.cc:715
Werror
void Werror(const char *fmt,...)
Definition: reporter.cc:188
idInitChoise
void idInitChoise(int r, int beg, int end, BOOLEAN *endch, int *choise)
Definition: simpleideals.cc:832
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:36
p_Vec2Polys
void p_Vec2Polys(poly v, poly **p, int *len, const ring r)
Definition: p_polys.cc:3524
mp_permmatrix::a_n
int a_n
Definition: matpol.cc:827
mp_permmatrix::mpSetSearch
void mpSetSearch(int s)
m
int m
Definition: cfEzgcd.cc:121
mp_Add
matrix mp_Add(matrix a, matrix b, const ring R)
Definition: matpol.cc:177
MATCOLS
#define MATCOLS(i)
Definition: matpol.h:27
mp_Monomials
void mp_Monomials(matrix c, int r, int var, matrix m, const ring R)
Definition: matpol.cc:363
assume
#define assume(x)
Definition: mod2.h:384
mp_Equal
BOOLEAN mp_Equal(matrix a, matrix b, const ring R)
Definition: matpol.cc:583
ip_smatrix::rank
long rank
Definition: matpol.h:19
NULL
#define NULL
Definition: omList.c:9
mp_permmatrix::piv_s
int piv_s
Definition: matpol.cc:827
p_SetComp
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:238
l
int l
Definition: cfEzgcd.cc:93
binom
int binom(int n, int r)
Definition: simpleideals.cc:912
R
#define R
Definition: sirandom.c:26
p_Setm
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:224
mp_DetBareiss
poly mp_DetBareiss(matrix a, const ring r)
returns the determinant of the matrix m; uses Bareiss algorithm
Definition: matpol.cc:1575
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
ip_smatrix::ncols
int ncols
Definition: matpol.h:21
p
int p
Definition: cfModGcd.cc:4019
p_IsConstant
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:1901
row_col_weight::ym
int ym
Definition: matpol.cc:790
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
count
int status int void size_t count
Definition: si_signals.h:58
mp_permmatrix::Xarray
poly * Xarray
Definition: matpol.cc:829
IDELEMS
#define IDELEMS(i)
Definition: simpleideals.h:24
mp_AddSubMat
static void mp_AddSubMat(ideal res, ideal sm, int col, const ring r)
Definition: matpol.cc:1738
mp_permmatrix::mpSaveArray
void mpSaveArray()
Definition: matpol.cc:851
p_ISet
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1286
p_Mult_q
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1035
mp_PolyWeight
static float mp_PolyWeight(poly p, const ring r)
Definition: matpol.cc:1201
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:48
mp_Compare
int mp_Compare(matrix a, matrix b, const ring R)
Definition: matpol.cc:564
prCopyR_NoSort
poly prCopyR_NoSort(poly p, ring src_r, ring dest_r)
Definition: prCopy.cc:77
n_Size
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:570
mp_Copy
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:62
id_Test
#define id_Test(A, lR)
Definition: simpleideals.h:79
MATROWS
#define MATROWS(i)
Definition: matpol.h:26
TraceOfProd
poly TraceOfProd(matrix a, matrix b, int n, const ring R)
Definition: matpol.cc:290
omFreeBin
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:257
p_Insert
static poly p_Insert(poly p1, poly p2, const ring R)
Definition: matpol.cc:611
A
#define A
Definition: sirandom.c:23
mp_Sub
matrix mp_Sub(matrix a, matrix b, const ring R)
Definition: matpol.cc:194
numbers.h
pNext
#define pNext(p)
Definition: monomials.h:40
pReverse
static poly pReverse(poly p)
Definition: p_polys.h:326
mp_PivRow
static int mp_PivRow(matrix a, int lr, int lc, const ring r)
Definition: matpol.cc:1293
omAlloc0
#define omAlloc0(size)
Definition: omAllocDecl.h:209
mp_permmatrix::mpGetSign
int mpGetSign()
Definition: matpol.cc:849
mp_permmatrix::mpSetElem
void mpSetElem(poly, int, int)
iiStringMatrix
char * iiStringMatrix(matrix im, int dim, const ring r, char ch)
Definition: matpol.cc:754