source: git/Singular/matpol.cc @ 8f1473

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