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

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