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