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