# source:git/libpolys/polys/matpol.cc@fe115e

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