source: git/Singular/matpol.cc @ 31cb04b

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