source: git/Singular/matpol.cc @ c616d1

spielwiese
Last change on this file since c616d1 was c616d1, checked in by Hans Schönemann <hannes@…>, 26 years ago
hannes/pohl: minor fix git-svn-id: file:///usr/local/Singular/svn/trunk@2529 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.20 1998-09-30 14:12:48 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 "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    b->m[i] = pCopy(t);
94  }
95  b->rank=a->rank;
96  return b;
97}
98
99/*2
100* make it a p * unit matrix
101*/
102matrix mpInitP(int r, int c, poly p)
103{
104  matrix rc = mpNew(r,c);
105  int i=min(r,c), n = c*(i-1)+i-1, inc = c+1;
106
107  pNormalize(p);
108  while (n>0)
109  {
110    rc->m[n] = pCopy(p);
111    n -= inc;
112  }
113  rc->m[0]=p;
114  return rc;
115}
116
117/*2
118* make it a v * unit matrix
119*/
120matrix mpInitI(int r, int c, int v)
121{
122  return mpInitP(r,c,pISet(v));
123}
124
125/*2
126* c = f*a
127*/
128matrix mpMultI(matrix a, int f)
129{
130  int k, n = a->nrows, m = a->ncols;
131  poly p = pISet(f);
132  matrix c = mpNew(n,m);
133
134  for (k=m*n-1; k>0; k--)
135    c->m[k] = pMult(pCopy(a->m[k]), pCopy(p));
136  c->m[0] = pMult(pCopy(a->m[0]), p);
137  return c;
138}
139
140/*2
141* multiply a matrix 'a' by a poly 'p', destroy the args
142*/
143matrix mpMultP(matrix a, poly p)
144{
145  int k, n = a->nrows, m = a->ncols;
146
147  pNormalize(p);
148  for (k=m*n-1; k>0; k--)
149    a->m[k] = pMult(a->m[k], pCopy(p));
150  a->m[0] = pMult(a->m[0], p);
151  return a;
152}
153
154matrix mpAdd(matrix a, matrix b)
155{
156  int k, n = a->nrows, m = a->ncols;
157  if ((n != b->nrows) || (m != b->ncols))
158  {
159/*
160*    Werror("cannot add %dx%d matrix and %dx%d matrix",
161*      m,n,b->cols(),b->rows());
162*/
163    return NULL;
164  }
165  matrix c = mpNew(n,m);
166  for (k=m*n-1; k>=0; k--)
167    c->m[k] = pAdd(pCopy(a->m[k]), pCopy(b->m[k]));
168  return c;
169}
170
171matrix mpSub(matrix a, matrix b)
172{
173  int k, n = a->nrows, m = a->ncols;
174  if ((n != b->nrows) || (m != b->ncols))
175  {
176/*
177*    Werror("cannot sub %dx%d matrix and %dx%d matrix",
178*      m,n,b->cols(),b->rows());
179*/
180    return NULL;
181  }
182  matrix c = mpNew(n,m);
183  for (k=m*n-1; k>=0; k--)
184    c->m[k] = pSub(pCopy(a->m[k]), pCopy(b->m[k]));
185  return c;
186}
187
188matrix mpMult(matrix a, matrix b)
189{
190  int i, j, k;
191  poly s, t, aik, bkj;
192  int m = MATROWS(a);
193  int p = MATCOLS(a);
194  int q = MATCOLS(b);
195
196  if (p!=MATROWS(b))
197  {
198/*
199*   Werror("cannot multiply %dx%d matrix and %dx%d matrix",
200*     m,p,b->rows(),q);
201*/
202    return NULL;
203  }
204  matrix c = mpNew(m,q);
205
206  for (i=1; i<=m; i++)
207  {
208    for (j=1; j<=q; j++)
209    {
210      t = NULL;
211      for (k=1; k<=p; k++)
212      {
213        aik = pCopy(MATELEM(a,i,k));
214        bkj = pCopy(MATELEM(b,k,j));
215        s = pMult(aik,bkj);
216        t = pAdd(t,s);
217      }
218      pNormalize(t);
219      MATELEM(c,i,j) = t;
220    }
221  }
222  return c;
223}
224
225matrix mpTransp(matrix a)
226{
227  int    i, j, r = MATROWS(a), c = MATCOLS(a);
228  poly *p;
229  matrix b =  mpNew(c,r);
230
231  p = b->m;
232  for (i=0; i<c; i++)
233  {
234    for (j=0; j<r; j++)
235    {
236      *p++ = pCopy(a->m[j*c+i]);
237    }
238  }
239  return b;
240}
241
242/*2
243*returns the trace of matrix a
244*/
245poly mpTrace ( matrix a)
246{
247  int i;
248  int n = (MATCOLS(a)<MATROWS(a)) ? MATCOLS(a) : MATROWS(a);
249  poly  t = NULL;
250
251  for (i=1; i<=n; i++)
252    t = pAdd(t, pCopy(MATELEM(a,i,i)));
253  return t;
254}
255
256/*2
257*returns the trace of the product of a and b
258*/
259poly TraceOfProd ( matrix a, matrix b, int n)
260{
261  int i, j;
262  poly  p, t = NULL;
263
264  for (i=1; i<=n; i++)
265  {
266    for (j=1; j<=n; j++)
267    {
268      p = pMult(pCopy(MATELEM(a,i,j)), pCopy(MATELEM(b,j,i)));
269      t = pAdd(t, p);
270    }
271  }
272  return t;
273}
274
275/*
276* C++ classes for Bareiss algorithm
277*/
278class row_col_weight
279{
280  private:
281  int ym, yn;
282  public:
283  float *wrow, *wcol;
284  row_col_weight() : ym(0) {}
285  row_col_weight(int, int);
286  ~row_col_weight();
287};
288
289/*2
290*  a submatrix M of a matrix X[m,n]:
291*    0 <= i < s_m <= a_m
292*    0 <= j < s_n <= a_n
293*    M = ( Xarray[qrow[i],qcol[j]] )
294*    if a_m = a_n and s_m = s_n
295*      det(X) = sign*div^(s_m-1)*det(M)
296*    resticted pivot for elimination
297*      0 <= j < piv_s
298*/
299class mp_permmatrix
300{
301  private:
302  int       a_m, a_n, s_m, s_n, sign, piv_s;
303  int       *qrow, *qcol;
304  poly      *Xarray;
305  void mpInitMat();
306  poly * mpRowAdr(int);
307  poly * mpColAdr(int);
308  void mpRowWeight(float *);
309  void mpColWeight(float *);
310  void mpRowSwap(int, int);
311  void mpColSwap(int, int);
312  public:
313  mp_permmatrix() : a_m(0) {}
314  mp_permmatrix(matrix);
315  mp_permmatrix(mp_permmatrix *);
316  ~mp_permmatrix();
317  int mpGetRow();
318  int mpGetCol();
319  int mpGetRdim();
320  int mpGetCdim();
321  int mpGetSign();
322  void mpSetSearch(int s);
323  void mpSaveArray();
324  poly mpGetElem(int, int);
325  void mpSetElem(poly, int, int);
326  void mpDelElem(int, int);
327  void mpElimBareiss(poly);
328  int mpPivotBareiss(row_col_weight *);
329  int mpPivotRow(row_col_weight *, int);
330  void mpToIntvec(intvec *);
331  void mpRowReorder();
332  void mpColReorder();
333};
334
335/*2
336*  caller of 'Bareiss' algorithm,
337*  return an list of a matrix and an intvec:
338*    the matrix is lower triangular and the result,
339*    the intvec is the performed permutation of columns.
340*/
341lists mpBareiss (matrix a, BOOLEAN sw)
342{
343  poly div;
344  matrix c = mpCopy(a);
345  mp_permmatrix *Bareiss = new mp_permmatrix(c);
346  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
347  intvec *v = new intvec(Bareiss->mpGetCdim());
348  lists res=(lists)Alloc(sizeof(slists));
349
350  if (sw) WarnS(feNotImplemented);
351  /* Bareiss */
352  div = NULL;
353  while(Bareiss->mpPivotBareiss(&w))
354  {
355    Bareiss->mpElimBareiss(div);
356    div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
357  }
358  Bareiss->mpToIntvec(v);
359  Bareiss->mpRowReorder();
360  Bareiss->mpColReorder();
361  Bareiss->mpSaveArray();
362  delete Bareiss;
363
364  res->Init(2);
365  res->m[0].rtyp=MATRIX_CMD;
366  res->m[0].data=(void *)c;
367  res->m[1].rtyp=INTVEC_CMD;
368  res->m[1].data=(void *)v;
369  return res;
370}
371
372/*2
373*  one step of 'Bareiss' algorithm
374*  for application in minor
375*  assume to have a full matrix
376*  if *H!=0, then divide by *H (the pivot from the step before)
377*  returns the choosen pivot *H=m[*r,*c]
378*  the result has the pivot at the right lower corner
379*/
380matrix mpOneStepBareiss (matrix a, poly *H, int *r, int *c)
381{
382  poly div=*H;
383  matrix re = mpCopy(a);
384  mp_permmatrix *Bareiss = new mp_permmatrix(re);
385  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
386  int row = *r;
387
388  /* step of Bareiss */
389  if(((row!=0) && Bareiss->mpPivotRow(&w,row-1)) || Bareiss->mpPivotBareiss(&w))
390  {
391    Bareiss->mpElimBareiss(div);
392    div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
393    pDelete(H);
394    *H = pCopy(div);
395    *c = Bareiss->mpGetCol()+1;
396    *r = Bareiss->mpGetRow()+1;
397    Bareiss->mpRowReorder();
398    Bareiss->mpColReorder();
399  }
400  else
401  {
402    pDelete(H);
403    *H = NULL;
404    *c = *r = 0;
405  }
406  Bareiss->mpSaveArray();
407  idTest((ideal)re);
408  delete Bareiss;
409  return re;
410}
411
412/*2
413*returns the determinant of the matrix m;
414*uses Bareiss algorithm
415*/
416poly mpDetBareiss (matrix a)
417{
418  int s;
419  poly div, res;
420  if (MATROWS(a) != MATCOLS(a))
421  {
422    Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
423    return NULL;
424  }
425  matrix c = mpCopy(a);
426  mp_permmatrix *Bareiss = new mp_permmatrix(c);
427  row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
428
429  /* Bareiss */
430  div = NULL;
431  while(Bareiss->mpPivotBareiss(&w))
432  {
433    Bareiss->mpElimBareiss(div);
434    div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
435  }
436  Bareiss->mpRowReorder();
437  Bareiss->mpColReorder();
438  Bareiss->mpSaveArray();
439  s = Bareiss->mpGetSign();
440  delete Bareiss;
441
442  /* result */
443  res = MATELEM(c,1,1);
444  MATELEM(c,1,1) = NULL;
445  idDelete((ideal *)&c);
446  if (s < 0)
447    res = pNeg(res);
448  return res;
449}
450
451/*2
452*returns the determinant of the matrix m;
453*uses Newtons formulea for symmetric functions
454*/
455poly mpDet (matrix m)
456{
457  int i,j,k,n;
458  poly p,q;
459  matrix a, s;
460  matrix ma[100];
461  number c=NULL, d=NULL, ONE=NULL;
462
463  n = MATROWS(m);
464  if (n != MATCOLS(m))
465  {
466    Werror("det of %d x %d matrix",n,MATCOLS(m));
467    return NULL;
468  }
469  k=currRing->ch;
470  if (k<0) k=-k;
471  else if (k==1) k=0;
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  matrix co;
784  if (f==NULL)
785  {
786    co = mpNew(2, 1);
787    MATELEM(co,1,1) = pOne();
788    MATELEM(co,2,1) = NULL;
789    return co;
790  }
791  sel = mpSelect(f, vars);
792  l = pLength(sel);
793  co = mpNew(2, l);
794  if (pOrdSgn==-1)
795  {
796    for (i=l; i>=1; i--)
797    {
798      h = pHead(sel);
799      MATELEM(co,1,i) = h;
800      MATELEM(co,2,i) = NULL;
801      sel = sel->next;
802    }
803  }
804  else
805  {
806    for (i=1; i<=l; i++)
807    {
808      h = pHead(sel);
809      MATELEM(co,1,i) = h;
810      MATELEM(co,2,i) = NULL;
811      sel = sel->next;
812    }
813  }
814  while (f!=NULL)
815  {
816    i = 1;
817    loop
818    {
819      h = mpExdiv(f, MATELEM(co,1,i));
820      if (h!=NULL)
821      {
822        MATELEM(co,2,i) = pAdd(MATELEM(co,2,i), h);
823        break;
824      }
825      if (i < l)
826        i++;
827      else
828        break;
829    }
830    pIter(f);
831  }
832  return co;
833}
834
835void mpCoef2(poly v, poly mon, matrix *c, matrix *m)
836{
837  polyset s;
838  poly p;
839  int sl,i,j;
840  int l=0;
841  poly sel=mpSelect(v,mon);
842
843  pVec2Polys(sel,&s,&sl);
844  for (i=0; i<sl; i++)
845    l=max(l,pLength(s[i]));
846  *c=mpNew(sl,l);
847  *m=mpNew(sl,l);
848  poly h;
849  int isConst;
850  for (j=1; j<=sl;j++)
851  {
852    p=s[j-1];
853    if (pIsConstant(p)) /*p != NULL */
854    {
855      isConst=-1;
856      i=l;
857    }
858    else
859    {
860      isConst=1;
861      i=1;
862    }
863    while(p!=NULL)
864    {
865      h = pHead(p);
866      MATELEM(*m,j,i) = h;
867      i+=isConst;
868      p = p->next;
869    }
870  }
871  while (v!=NULL)
872  {
873    i = 1;
874    j = pGetComp(v);
875    loop
876    {
877      poly mp=MATELEM(*m,j,i);
878      if (mp!=NULL)
879      {
880        h = mpExdiv(v, mp /*MATELEM(*m,j,i)*/);
881        if (h!=NULL)
882        {
883          pSetComp(h,0);
884          MATELEM(*c,j,i) = pAdd(MATELEM(*c,j,i), h);
885          break;
886        }
887      }
888      if (i < l)
889        i++;
890      else
891        break;
892    }
893    v = v->next;
894  }
895}
896
897
898BOOLEAN mpEqual(matrix a, matrix b)
899{
900  if ((MATCOLS(a)!=MATCOLS(b)) || (MATROWS(a)!=MATROWS(b)))
901    return FALSE;
902  int i=MATCOLS(a)*MATROWS(b)-1;
903  while (i>=0)
904  {
905    if (a->m[i]==NULL)
906    {
907      if (b->m[i]!=NULL) return FALSE;
908    }
909    else
910      if (pComp(a->m[i],b->m[i])!=0) return FALSE;
911    i--;
912  }
913  i=MATCOLS(a)*MATROWS(b)-1;
914  while (i>=0)
915  {
916    if (pSub(pCopy(a->m[i]),pCopy(b->m[i]))!=NULL) return FALSE;
917    i--;
918  }
919  return TRUE;
920}
921
922/* --------------- internal stuff ------------------- */
923
924row_col_weight::row_col_weight(int i, int j)
925{
926  ym = i;
927  yn = j;
928  wrow = (float *)Alloc(i*sizeof(float));
929  wcol = (float *)Alloc(j*sizeof(float));
930}
931
932row_col_weight::~row_col_weight()
933{
934  if (ym!=0)
935  {
936    Free((ADDRESS)wcol, yn*sizeof(float));
937    Free((ADDRESS)wrow, ym*sizeof(float));
938  }
939}
940
941mp_permmatrix::mp_permmatrix(matrix A) : sign(1)
942{
943  a_m = A->nrows;
944  a_n = A->ncols;
945  this->mpInitMat();
946  Xarray = A->m;
947}
948
949mp_permmatrix::mp_permmatrix(mp_permmatrix *M)
950{
951  poly p, *athis, *aM;
952  int i, j;
953
954  a_m = M->s_m;
955  a_n = M->s_n;
956  sign = M->sign;
957  this->mpInitMat();
958  Xarray = (poly *)Alloc0(a_m*a_n*sizeof(poly));
959  for (i=a_m-1; i>=0; i--)
960  {
961    athis = this->mpRowAdr(i);
962    aM = M->mpRowAdr(i);
963    for (j=a_n-1; j>=0; j--)
964    {
965      p = aM[M->qcol[j]];
966      if (p)
967      {
968        athis[j] = pCopy(p);
969      }
970    }
971  }
972}
973
974mp_permmatrix::~mp_permmatrix()
975{
976  int k;
977
978  if (a_m != 0)
979  {
980    Free((ADDRESS)qrow,a_m*sizeof(int));
981    Free((ADDRESS)qcol,a_n*sizeof(int));
982    if (Xarray != NULL)
983    {
984      for (k=a_m*a_n-1; k>=0; k--)
985        pDelete(&Xarray[k]);
986      Free((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
987    }
988  }
989}
990
991int mp_permmatrix::mpGetRdim() { return s_m; }
992
993int mp_permmatrix::mpGetCdim() { return s_n; }
994
995int mp_permmatrix::mpGetSign() { return sign; }
996
997void mp_permmatrix::mpSetSearch(int s) { piv_s = s; }
998
999void mp_permmatrix::mpSaveArray() { Xarray = NULL; }
1000
1001poly mp_permmatrix::mpGetElem(int r, int c)
1002{
1003  return Xarray[a_n*qrow[r]+qcol[c]];
1004}
1005
1006void mp_permmatrix::mpSetElem(poly p, int r, int c)
1007{
1008  Xarray[a_n*qrow[r]+qcol[c]] = p;
1009}
1010
1011void mp_permmatrix::mpDelElem(int r, int c)
1012{
1013  pDelete(&Xarray[a_n*qrow[r]+qcol[c]]);
1014}
1015
1016/*
1017* the Bareiss-type elimination with division by div (div != NULL)
1018*/
1019void mp_permmatrix::mpElimBareiss(poly div)
1020{
1021  poly piv, elim, q1, q2, *ap, *a;
1022  int i, j, jj;
1023
1024  ap = this->mpRowAdr(s_m);
1025  piv = ap[qcol[s_n]];
1026  for(i=s_m-1; i>=0; i--)
1027  {
1028    a = this->mpRowAdr(i);
1029    elim = a[qcol[s_n]];
1030    if (elim != NULL)
1031    {
1032      for (j=s_n-1; j>=0; j--)
1033      {
1034        q2 = NULL;
1035        jj = qcol[j];
1036        q1 = a[jj];
1037        if (ap[jj] != NULL)
1038        {
1039          q2 = pNeg(pCopy(ap[jj]));
1040          q2 = pMult(q2, pCopy(elim));
1041          if (q1 != NULL)
1042          {
1043            q1 = pMult(q1,pCopy(piv));
1044            q2 = pAdd(q2, q1);
1045          }
1046        }
1047        else if (q1 != NULL)
1048        {
1049          q2 = pMult(q1, pCopy(piv));
1050        }
1051        if ((q2!=NULL) && div)
1052          q2 = mpDivide(q2, div);
1053        a[jj] = q2;
1054      }
1055      pDelete(&a[qcol[s_n]]);
1056    }
1057    else
1058    {
1059      for (j=s_n-1; j>=0; j--)
1060      {
1061        jj = qcol[j];
1062        q1 = a[jj];
1063        if (q1 != NULL)
1064        {
1065          q1 = pMult(q1, pCopy(piv));
1066          if (div)
1067            q1 = mpDivide(q1, div);
1068          a[jj] = q1;
1069        }
1070      }
1071    }
1072  }
1073}
1074
1075/*2
1076* pivot strategy for Bareiss algorithm
1077*/
1078int mp_permmatrix::mpPivotBareiss(row_col_weight *C)
1079{
1080  poly p, *a;
1081  int i, j, iopt, jopt;
1082  float sum, f1, f2, fo, r, ro, lp;
1083  float *dr = C->wrow, *dc = C->wcol;
1084
1085  fo = 1.0e20;
1086  ro = 0.0;
1087  iopt = jopt = -1;
1088
1089  s_n--;
1090  s_m--;
1091  if (s_m == 0)
1092    return 0;
1093  if (s_n == 0)
1094  {
1095    for(i=s_m; i>=0; i--)
1096    {
1097      p = this->mpRowAdr(i)[qcol[0]];
1098      if (p)
1099      {
1100        f1 = mpPolyWeight(p);
1101        if (f1 < fo)
1102        {
1103          fo = f1;
1104          if (iopt >= 0)
1105            pDelete(&(this->mpRowAdr(iopt)[qcol[0]]));
1106          iopt = i;
1107        }
1108        else
1109          pDelete(&(this->mpRowAdr(i)[qcol[0]]));
1110      }
1111    }
1112    if (iopt >= 0)
1113      mpReplace(iopt, s_m, sign, qrow);
1114    return 0;
1115  }
1116  this->mpRowWeight(dr);
1117  this->mpColWeight(dc);
1118  sum = 0.0;
1119  for(i=s_m; i>=0; i--)
1120    sum += dr[i];
1121  for(i=s_m; i>=0; i--)
1122  {
1123    r = dr[i];
1124    a = this->mpRowAdr(i);
1125    for(j=s_n; j>=0; j--)
1126    {
1127      p = a[qcol[j]];
1128      if (p)
1129      {
1130        lp = mpPolyWeight(p);
1131        ro = r - lp;
1132        f1 = ro * (dc[j]-lp);
1133        if (f1 != 0.0)
1134        {
1135          f2 = lp * (sum - ro - dc[j]);
1136          f2 += f1;
1137        }
1138        else
1139          f2 = lp-r-dc[j];
1140        if (f2 < fo)
1141        {
1142          fo = f2;
1143          iopt = i;
1144          jopt = j;
1145        }
1146      }
1147    }
1148  }
1149  if (iopt < 0)
1150    return 0;
1151  mpReplace(iopt, s_m, sign, qrow);
1152  mpReplace(jopt, s_n, sign, qcol);
1153  return 1;
1154}
1155
1156/*2
1157* pivot strategy for Bareiss algorithm with defined row
1158*/
1159int mp_permmatrix::mpPivotRow(row_col_weight *C, int row)
1160{
1161  poly p, *a;
1162  int j, iopt, jopt;
1163  float sum, f1, f2, fo, r, ro, lp;
1164  float *dr = C->wrow, *dc = C->wcol;
1165
1166  fo = 1.0e20;
1167  ro = 0.0;
1168  iopt = jopt = -1;
1169
1170  s_n--;
1171  s_m--;
1172  if (s_m == 0)
1173    return 0;
1174  if (s_n == 0)
1175  {
1176    p = this->mpRowAdr(row)[qcol[0]];
1177    if (p)
1178    {
1179      f1 = mpPolyWeight(p);
1180      if (f1 < fo)
1181      {
1182        fo = f1;
1183        if (iopt >= 0)
1184        pDelete(&(this->mpRowAdr(iopt)[qcol[0]]));
1185        iopt = row;
1186      }
1187      else
1188        pDelete(&(this->mpRowAdr(row)[qcol[0]]));
1189    }
1190    if (iopt >= 0)
1191      mpReplace(iopt, s_m, sign, qrow);
1192    return 0;
1193  }
1194  this->mpRowWeight(dr);
1195  this->mpColWeight(dc);
1196  sum = 0.0;
1197  for(j=s_m; j>=0; j--)
1198    sum += dr[j];
1199  r = dr[row];
1200  a = this->mpRowAdr(row);
1201  for(j=s_n; j>=0; j--)
1202  {
1203    p = a[qcol[j]];
1204    if (p)
1205    {
1206      lp = mpPolyWeight(p);
1207      ro = r - lp;
1208      f1 = ro * (dc[j]-lp);
1209      if (f1 != 0.0)
1210      {
1211        f2 = lp * (sum - ro - dc[j]);
1212        f2 += f1;
1213      }
1214      else
1215        f2 = lp-r-dc[j];
1216      if (f2 < fo)
1217      {
1218        fo = f2;
1219        iopt = row;
1220        jopt = j;
1221      }
1222    }
1223  }
1224  if (iopt < 0)
1225    return 0;
1226  mpReplace(iopt, s_m, sign, qrow);
1227  mpReplace(jopt, s_n, sign, qcol);
1228  return 1;
1229}
1230
1231void mp_permmatrix::mpToIntvec(intvec *v)
1232{
1233  int i;
1234
1235  for (i=v->rows()-1; i>=0; i--)
1236    (*v)[i] = qcol[i]+1;
1237}
1238
1239void mp_permmatrix::mpRowReorder()
1240{
1241  int k, i, i1, i2;
1242
1243  if (a_m > a_n)
1244    k = a_m - a_n;
1245  else
1246    k = 0;
1247  for (i=a_m-1; i>=k; i--)
1248  {
1249    i1 = qrow[i];
1250    if (i1 != i)
1251    {
1252      this->mpRowSwap(i1, i);
1253      i2 = 0;
1254      while (qrow[i2] != i) i2++;
1255      qrow[i2] = i1;
1256    }
1257  }
1258}
1259
1260void mp_permmatrix::mpColReorder()
1261{
1262  int k, j, j1, j2;
1263
1264  if (a_n > a_m)
1265    k = a_n - a_m;
1266  else
1267    k = 0;
1268  for (j=a_n-1; j>=k; j--)
1269  {
1270    j1 = qcol[j];
1271    if (j1 != j)
1272    {
1273      this->mpColSwap(j1, j);
1274      j2 = 0;
1275      while (qcol[j2] != j) j2++;
1276      qcol[j2] = j1;
1277    }
1278  }
1279}
1280
1281// private
1282void mp_permmatrix::mpInitMat()
1283{
1284  int k;
1285
1286  s_m = a_m;
1287  s_n = a_n;
1288  piv_s = 0;
1289  qrow = (int *)Alloc(a_m*sizeof(int));
1290  qcol = (int *)Alloc(a_n*sizeof(int));
1291  for (k=a_m-1; k>=0; k--) qrow[k] = k;
1292  for (k=a_n-1; k>=0; k--) qcol[k] = k;
1293}
1294
1295poly * mp_permmatrix::mpRowAdr(int r)
1296{
1297  return &(Xarray[a_n*qrow[r]]);
1298}
1299
1300poly * mp_permmatrix::mpColAdr(int c)
1301{
1302  return &(Xarray[qcol[c]]);
1303}
1304
1305void mp_permmatrix::mpRowWeight(float *wrow)
1306{
1307  poly p, *a;
1308  int i, j;
1309  float count;
1310
1311  for (i=s_m; i>=0; i--)
1312  {
1313    a = this->mpRowAdr(i);
1314    count = 0.0;
1315    for(j=s_n; j>=0; j--)
1316    {
1317      p = a[qcol[j]];
1318      if (p)
1319        count += mpPolyWeight(p);
1320    }
1321    wrow[i] = count;
1322  }
1323}
1324
1325void mp_permmatrix::mpColWeight(float *wcol)
1326{
1327  poly p, *a;
1328  int i, j;
1329  float count;
1330
1331  for (j=s_n; j>=0; j--)
1332  {
1333    a = this->mpColAdr(j);
1334    count = 0.0;
1335    for(i=s_m; i>=0; i--)
1336    {
1337      p = a[a_n*qrow[i]];
1338      if (p)
1339        count += mpPolyWeight(p);
1340    }
1341    wcol[j] = count;
1342  }
1343}
1344
1345void mp_permmatrix::mpRowSwap(int i1, int i2)
1346{
1347   poly p, *a1, *a2;
1348   int j;
1349
1350   a1 = &(Xarray[a_n*i1]);
1351   a2 = &(Xarray[a_n*i2]);
1352   for (j=a_n-1; j>= 0; j--)
1353   {
1354     p = a1[j];
1355     a1[j] = a2[j];
1356     a2[j] = p;
1357   }
1358}
1359
1360void mp_permmatrix::mpColSwap(int j1, int j2)
1361{
1362   poly p, *a1, *a2;
1363   int i, k = a_n*a_m;
1364
1365   a1 = &(Xarray[j1]);
1366   a2 = &(Xarray[j2]);
1367   for (i=0; i< k; i+=a_n)
1368   {
1369     p = a1[i];
1370     a1[i] = a2[i];
1371     a2[i] = p;
1372   }
1373}
1374
1375int mp_permmatrix::mpGetRow()
1376{
1377  return qrow[s_m];
1378}
1379
1380int mp_permmatrix::mpGetCol()
1381{
1382  return qcol[s_n];
1383}
1384
1385/*2
1386* exact division a/b, used in Bareiss algorithm
1387* a destroyed, b NOT destroyed
1388*/
1389static poly mpDivide(poly a, poly b)
1390{
1391  number x, y, z;
1392  poly r, tail, t, h, h0;
1393  int i, deg;
1394
1395  r = a;
1396  x = pGetCoeff(b);
1397  deg = pTotaldegree(b);
1398  tail = pNext(b);
1399  if (tail == NULL)
1400  {
1401    do
1402    {
1403      if (deg != 0)
1404      {
1405        for (i=pVariables; i; i--)
1406          pSubExp(r,i,  pGetExp(b,i));
1407        pSetm(r);
1408      }
1409      y = nDiv(pGetCoeff(r),x);
1410      nNormalize(y);
1411      pSetCoeff(r,y);
1412      pIter(r);
1413    } while (r != NULL);
1414    //pTest(a);
1415    return a;
1416  }
1417  h0 = pInit();
1418  do
1419  {
1420    if (deg != 0)
1421    {
1422      for (i=pVariables; i>0; i--)
1423        pSubExp(r,i,pGetExp(b,i));
1424      pSetm(r);
1425    }
1426    y = nDiv(pGetCoeff(r), x);
1427    nNormalize(y);
1428    pSetCoeff(r,y);
1429    t = tail;
1430    h = h0;
1431    do
1432    {
1433      h = pNext(h) = pInit();
1434      for (i=pVariables; i>0; i--)
1435        pSetExp(h,i, pGetExp(r,i)+pGetExp(t,i));
1436      pSetm(h);
1437      z = nMult(y, pGetCoeff(t));
1438      pSetCoeff0(h,nNeg(z));
1439      pIter(t);
1440    } while (t != NULL);
1441    pNext(h) = NULL;
1442    r = pNext(r) = pAdd(pNext(r),pNext(h0));
1443  } while (r!=NULL);
1444  pFree1(h0);
1445  //pTest(a);
1446  return a;
1447}
1448
1449/*
1450* perform replacement for pivot strategy in Bareiss algorithm
1451* change sign of determinant
1452*/
1453static void mpReplace(int j, int n, int &sign, int *perm)
1454{
1455  int k;
1456
1457  if (j != n)
1458  {
1459    k = perm[n];
1460    perm[n] = perm[j];
1461    perm[j] = k;
1462    sign = -sign;
1463  }
1464}
1465
1466/*
1467* weigth of a polynomial, for pivot strategy
1468* modify this for characteristic 0 !!!
1469*/
1470static float mpPolyWeight(poly p)
1471{
1472  int i;
1473  float res;
1474
1475  if (pNext(p) == NULL)
1476  {
1477    res = (float)nSize(pGetCoeff(p));
1478    if (pTotaldegree(p) != 0) res += 1.0;
1479  }
1480  else
1481  {
1482    i = 0;
1483    res = 0.0;
1484    do
1485    {
1486      i++;
1487      res += (float)nSize(pGetCoeff(p));
1488      pIter(p);
1489    }
1490    while (p);
1491    res += (float)i;
1492  }
1493  return res;
1494}
1495
1496static int mpNextperm(perm * z, int max)
1497{
1498  int s, i, k, t;
1499  s = max;
1500  do
1501  {
1502    s--;
1503  }
1504  while ((s > 0) && ((*z)[s] >= (*z)[s+1]));
1505  if (s==0)
1506    return 0;
1507  do
1508  {
1509    (*z)[s]++;
1510    k = 0;
1511    do
1512    {
1513      k++;
1514    }
1515    while (((*z)[k] != (*z)[s]) && (k!=s));
1516  }
1517  while (k < s);
1518  for (i=s+1; i <= max; i++)
1519  {
1520    (*z)[i]=0;
1521    do
1522    {
1523      (*z)[i]++;
1524      k=0;
1525      do
1526      {
1527        k++;
1528      }
1529      while (((*z)[k] != (*z)[i]) && (k != i));
1530    }
1531    while (k < i);
1532  }
1533  s = max+1;
1534  do
1535  {
1536    s--;
1537  }
1538  while ((s > 0) && ((*z)[s] > (*z)[s+1]));
1539  t = 1;
1540  for (i=1; i<max; i++)
1541    for (k=i+1; k<=max; k++)
1542      if ((*z)[k] < (*z)[i])
1543        t = -t;
1544  (*z)[0] = t;
1545  return s;
1546}
1547
1548static poly mpLeibnitz(matrix a)
1549{
1550  int i, e, n;
1551  poly p, d;
1552  perm z;
1553
1554  n = MATROWS(a);
1555  memset(&z,0,(n+2)*sizeof(int));
1556  p = pOne();
1557  for (i=1; i <= n; i++)
1558    p = pMult(p, pCopy(MATELEM(a, i, i)));
1559  d = p;
1560  for (i=1; i<= n; i++)
1561    z[i] = i;
1562  z[0]=1;
1563  e = 1;
1564  if (n!=1)
1565  {
1566    while (e)
1567    {
1568      e = mpNextperm((perm *)&z, n);
1569      p = pOne();
1570      for (i = 1; i <= n; i++)
1571        p = pMult(p, pCopy(MATELEM(a, i, z[i])));
1572      if (z[0] > 0)
1573        d = pAdd(d, p);
1574      else
1575        d = pSub(d, p);
1576    }
1577  }
1578  return d;
1579}
1580
1581static poly minuscopy (poly p)
1582{
1583  poly w;
1584  number  e;
1585  e = nInit(-1);
1586  w = pCopy(p);
1587  pMultN(w, e);
1588  nDelete(&e);
1589  return w;
1590}
1591
1592/*2
1593* insert a monomial into a list, avoid duplicates
1594* arguments are destroyed
1595*/
1596static poly pInsert(poly p1, poly p2)
1597{
1598  poly a1, p, a2, a;
1599  int c;
1600
1601  if (p1==NULL) return p2;
1602  if (p2==NULL) return p1;
1603  a1 = p1;
1604  a2 = p2;
1605  a = p  = pOne();
1606  loop
1607  {
1608    c = pComp(a1, a2);
1609    if (c == 1)
1610    {
1611      a = pNext(a) = a1;
1612      pIter(a1);
1613      if (a1==NULL)
1614      {
1615        pNext(a) = a2;
1616        break;
1617      }
1618    }
1619    else if (c == -1)
1620    {
1621      a = pNext(a) = a2;
1622      pIter(a2);
1623      if (a2==NULL)
1624      {
1625        pNext(a) = a1;
1626        break;
1627      }
1628    }
1629    else
1630    {
1631      pDelete1(&a2);
1632      a = pNext(a) = a1;
1633      pIter(a1);
1634      if (a1==NULL)
1635      {
1636        pNext(a) = a2;
1637        break;
1638      }
1639      else if (a2==NULL)
1640      {
1641        pNext(a) = a1;
1642        break;
1643      }
1644    }
1645  }
1646  pDelete1(&p);
1647  return p;
1648}
1649
1650/*2
1651*if what == xy the result is the list of all different power products
1652*    x^i*y^j (i, j >= 0) that appear in fro
1653*/
1654static poly mpSelect (poly fro, poly what)
1655{
1656  int i;
1657  poly h, res;
1658  res = NULL;
1659  while (fro!=NULL)
1660  {
1661    h = pOne();
1662    for (i=1; i<=pVariables; i++)
1663      pSetExp(h,i, pGetExp(fro,i) * pGetExp(what, i));
1664    pSetComp(h, pGetComp(fro));
1665    pSetm(h);
1666    res = pInsert(h, res);
1667    fro = fro->next;
1668  }
1669  return res;
1670}
1671
1672/*2
1673*exact divisor: let d  == x^i*y^j, m is thought to have only one term;
1674*    return m/d iff d divides m, and no x^k*y^l (k>i or l>j) divides m
1675*/
1676static poly mpExdiv ( poly m, poly d)
1677{
1678  int i;
1679  poly h = pHead(m);
1680  for (i=1; i<=pVariables; i++)
1681  {
1682    if (pGetExp(d,i) > 0)
1683    {
1684      if (pGetExp(d,i) != pGetExp(h,i))
1685      {
1686        pDelete(&h);
1687        return NULL;
1688      }
1689      pSetExp(h,i,0);
1690    }
1691  }
1692  pSetm(h);
1693  return h;
1694}
1695
Note: See TracBrowser for help on using the repository browser.