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