source: git/Singular/matpol.cc @ 194f5e5

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