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