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