source: git/libpolys/polys/matpol.cc @ 599813

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